// from sal/inc/rtl/math.hxx #include #include #include #include #include using namespace std; class IllegalArgument {}; /** Test equality of two values with an accuracy of the magnitude of the given values scaled by 2^-48 (4 bits roundoff stripped). @ATTENTION approxEqual( value!=0.0, 0.0 ) _never_ yields true. */ inline bool approxEqual(double a, double b) { if ( a == b ) return true; double x = a - b; return (x < 0.0 ? -x : x) < ((a < 0.0 ? -a : a) * (1.0 / (16777216.0 * 16777216.0))); } /** floor() method taking approxEqual() into account. Use for expected integer values being calculated by double functions. @ATTENTION The threshhold value is 3.55271e-015 */ inline double approxFloor(double a) { double b = floor( a ); // The second approxEqual() is necessary for values that are near the limit // of numbers representable with 4 bits stripped off. (#i12446#) if ( approxEqual( a - 1.0, b ) && !approxEqual( a, b ) ) return b + 1.0; return b; } // from interpr3.cxx double BinomKoeff(double n, double k) { double nVal = 0.0; k = approxFloor(k); if (n < k) nVal = 0.0; else if (k == 0.0) nVal = 1.0; else { nVal = n/k; n--; k--; while (k > 0.0) { nVal *= n/k; k--; n--; } } return nVal; } double ScHypGeomDist( double N, double M, double n, double x ) { // N n_population // M successes // n n_sample // x X if( (x < 0.0) || (n < x) || (M < x) || (N < n) || (N < M) || (x < n - N + M) ) throw IllegalArgument(); double fFactor = BinomKoeff( n, x ) / BinomKoeff( N, M ) * BinomKoeff( N - n, M - x ); return fFactor; } void insertNumberIntoDenom( unsigned long k, multiset& setEnum, multiset& setDenom ) { multiset::iterator itr = setEnum.find( k ); if ( itr != setEnum.end() ) setEnum.erase( itr ); else setDenom.insert( k ); } void SetBinomElements( unsigned long n, unsigned long k, multiset& setEnum, multiset& setDenom ) { if ( n >= k ) { if ( n == 0 || k == 0 ) return; if ( n > 1 ) setEnum.insert( n ); if ( k > 1 ) insertNumberIntoDenom( k, setEnum, setDenom ); --n; --k; while ( k > 0 ) { if ( n > 1 ) setEnum.insert( n ); if ( k > 1 ) insertNumberIntoDenom( k, setEnum, setDenom ); --n; --k; } } } double getHypGeomDist( double N, double M, double n, double x ) { // N n_population // M successes // n n_sample // x X if( (x < 0.0) || (n < x) || (M < x) || (N < n) || (N < M) || (x < n - N + M) ) throw IllegalArgument(); typedef multiset< unsigned long > LongMultiSet; LongMultiSet setEnum, setDenom; SetBinomElements( static_cast(n), static_cast(x), setEnum, setDenom ); SetBinomElements( static_cast(N), static_cast(M), setDenom, setEnum ); SetBinomElements( static_cast(N-n), static_cast(M-x), setEnum, setDenom ); double fAns = 1.0; LongMultiSet::reverse_iterator it1 = setEnum.rbegin(), it1End = setEnum.rend(); LongMultiSet::reverse_iterator it2 = setDenom.rbegin(), it2End = setDenom.rend(); for ( ; it1 != it1End || it2 != it2End; ) { double fEnum = 1.0, fDenom = 1.0; if ( it1 != it1End ) fEnum = static_cast( *it1++ ); if ( it2 != it2End ) fDenom = static_cast( *it2++ ); fAns *= fEnum / fDenom; } return fAns; } int main( int argc, char** argv ) { double fPopUpper = 100.0; cout << setprecision( 16 ); if ( argc < 5 ) { // x >= 0.0 && x <= n && x <= M // n <= N // M <= N // x >= n - N + M >= 0 double fMaxDiff = 0.0; unsigned short nCount = 0; unsigned long nTotalCount = 0; unsigned long nNumDiff = 0; for ( double N = 0.0; N < fPopUpper; ++N ) for ( double M = 0.0; M <= N; ++M ) for ( double n = N - M; n <= N; ++n ) for ( double x = n - N + M; x <= n && x <= M; ++x ) { double fOld = ScHypGeomDist( N, M, n, x ); // old method double fNew = getHypGeomDist( N, M, n, x ); // new method double fAbsDiff = fabs( fOld - fNew ); if ( !approxEqual( fOld, fNew ) ) ++nNumDiff; fMaxDiff = fMaxDiff < fAbsDiff ? fAbsDiff : fMaxDiff; if ( nCount > 10000 ) { nCount = 0; cout << "max diff: " << fMaxDiff << endl; } ++nCount; ++nTotalCount; } cout << "maximum difference between old and new methods: " << fMaxDiff << endl; cout << "number of times comparisons returned a non-equal result: " << nNumDiff << endl; cout << "number of comparisons performed: " << nTotalCount << endl; } else { double N = atof( argv[4] ); // population double M = atof( argv[3] ); // # successes double n = atof( argv[2] ); // # sample double x = atof( argv[1] ); // i = x cout << "population (N) = " << N << endl; cout << " # success (M) = " << M << endl; cout << " # samples (n) = " << n << endl; cout << " x = " << x << endl; try { double fOld = ScHypGeomDist( N, M, n, x ); cout << "old: " << fOld << endl; double fNew = getHypGeomDist( N, M, n, x ); cout << "new: " << fNew << endl; } catch( const IllegalArgument& ex ) { cout << "illegal argument" << endl; } } return EXIT_SUCCESS; }