From 468a7ed080e36762cd80a042c4f841ff9cc3af60 Mon Sep 17 00:00:00 2001 From: Martin Date: Thu, 14 May 2015 20:20:51 +0200 Subject: [PATCH] std.internal.math.gammafunction: loosen some tests for 64-bit reals. --- std/internal/math/gammafunction.d | 71 ++++++++++++++++++------------- 1 file changed, 41 insertions(+), 30 deletions(-) diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index 3f1b997e8d1..9515212bd83 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -523,7 +523,10 @@ unittest { assert(logGamma(-50.2) == log(fabs(gamma(-50.2)))); assert(logGamma(-0.008) == log(fabs(gamma(-0.008)))); assert(feqrel(logGamma(-38.8),log(fabs(gamma(-38.8)))) > real.mant_dig-4); - assert(feqrel(logGamma(1500.0L),log(gamma(1500.0L))) > real.mant_dig-2); + static if (real.mant_dig >= 64) // incl. 80-bit reals + assert(feqrel(logGamma(1500.0L),log(gamma(1500.0L))) > real.mant_dig-2); + else static if (real.mant_dig >= 53) // incl. 64-bit reals + assert(feqrel(logGamma(150.0L),log(gamma(150.0L))) > real.mant_dig-2); } @@ -898,12 +901,15 @@ unittest { // also tested by the normal distribution // Test against Mathematica betaRegularized[z,a,b] // These arbitrary points are chosen to give good code coverage. assert(feqrel(betaIncomplete(8, 10, 0.2), 0.010_934_315_234_099_2L) >= real.mant_dig - 5); - assert(feqrel(betaIncomplete(2, 2.5, 0.9),0.989_722_597_604_452_767_171_003_59L) >= real.mant_dig - 1 ); - assert(feqrel(betaIncomplete(1000, 800, 0.5), 1.179140859734704555102808541457164E-06L) >= real.mant_dig - 13 ); - assert(feqrel(betaIncomplete(0.0001, 10000, 0.0001),0.999978059362107134278786L) >= real.mant_dig - 18 ); + assert(feqrel(betaIncomplete(2, 2.5, 0.9), 0.989_722_597_604_452_767_171_003_59L) >= real.mant_dig - 1); + static if (real.mant_dig >= 64) // incl. 80-bit reals + assert(feqrel(betaIncomplete(1000, 800, 0.5), 1.179140859734704555102808541457164E-06L) >= real.mant_dig - 13); + else + assert(feqrel(betaIncomplete(1000, 800, 0.5), 1.179140859734704555102808541457164E-06L) >= real.mant_dig - 14); + assert(feqrel(betaIncomplete(0.0001, 10000, 0.0001), 0.999978059362107134278786L) >= real.mant_dig - 18); assert(betaIncomplete(0.01, 327726.7, 0.545113) == 1.0); assert(feqrel(betaIncompleteInv(8, 10, 0.010_934_315_234_099_2L), 0.2L) >= real.mant_dig - 2); - assert(feqrel(betaIncomplete(0.01, 498.437, 0.0121433),0.99999664562033077636065L) >= real.mant_dig - 1); + assert(feqrel(betaIncomplete(0.01, 498.437, 0.0121433), 0.99999664562033077636065L) >= real.mant_dig - 1); assert(feqrel(betaIncompleteInv(5, 10, 0.2000002972865658842), 0.229121208190918L) >= real.mant_dig - 3); assert(feqrel(betaIncompleteInv(4, 7, 0.8000002209179505L), 0.483657360076904L) >= real.mant_dig - 3); @@ -913,30 +919,32 @@ unittest { // also tested by the normal distribution // Extensive testing failed to increase the coverage. It seems likely that about // half the code in this function is unnecessary; there is potential for // significant improvement over the original CEPHES code. - - assert(betaIncompleteInv(0.01, 8e-48, 5.45464e-20)==1-real.epsilon); - assert(betaIncompleteInv(0.01, 8e-48, 9e-26)==1-real.epsilon); - - // Beware: a one-bit change in pow() changes almost all digits in the result! - assert(feqrel(betaIncompleteInv(0x1.b3d151fbba0eb18p+1, 1.2265e-19, 2.44859e-18),0x1.c0110c8531d0952cp-1L) > 10); - // This next case uncovered a one-bit difference in the FYL2X instruction - // between Intel and AMD processors. This difference gets magnified by 2^^38. - // WolframAlpha crashes attempting to calculate this. - assert(feqrel(betaIncompleteInv(0x1.ff1275ae5b939bcap-41, 4.6713e18, 0.0813601), - 0x1.f97749d90c7adba8p-63L) >= real.mant_dig - 39); - real a1 = 3.40483; - assert(betaIncompleteInv(a1, 4.0640301659679627772e19L, 0.545113)== 0x1.ba8c08108aaf5d14p-109); - real b1 = 2.82847e-25; - assert(feqrel(betaIncompleteInv(0.01, b1, 9e-26), 0x1.549696104490aa9p-830L) >= real.mant_dig-10); - - // --- Problematic cases --- - // This is a situation where the series expansion fails to converge - assert( isNaN(betaIncompleteInv(0.12167, 4.0640301659679627772e19L, 0.0813601))); - // This next result is almost certainly erroneous. - // Mathematica states: "(cannot be determined by current methods)" - assert(betaIncomplete(1.16251e20, 2.18e39, 5.45e-20)==-real.infinity); - // WolframAlpha gives no result for this, though indicates that it approximately 1.0 - 1.3e-9 - assert(1- betaIncomplete(0.01, 328222, 4.0375e-5) == 0x1.5f62926b4p-30); + static if (real.mant_dig == 64) // 80-bit reals + { + assert(betaIncompleteInv(0.01, 8e-48, 5.45464e-20) == 1-real.epsilon); + assert(betaIncompleteInv(0.01, 8e-48, 9e-26) == 1-real.epsilon); + + // Beware: a one-bit change in pow() changes almost all digits in the result! + assert(feqrel(betaIncompleteInv(0x1.b3d151fbba0eb18p+1, 1.2265e-19, 2.44859e-18), 0x1.c0110c8531d0952cp-1L) > 10); + // This next case uncovered a one-bit difference in the FYL2X instruction + // between Intel and AMD processors. This difference gets magnified by 2^^38. + // WolframAlpha crashes attempting to calculate this. + assert(feqrel(betaIncompleteInv(0x1.ff1275ae5b939bcap-41, 4.6713e18, 0.0813601), + 0x1.f97749d90c7adba8p-63L) >= real.mant_dig - 39); + real a1 = 3.40483; + assert(betaIncompleteInv(a1, 4.0640301659679627772e19L, 0.545113) == 0x1.ba8c08108aaf5d14p-109); + real b1 = 2.82847e-25; + assert(feqrel(betaIncompleteInv(0.01, b1, 9e-26), 0x1.549696104490aa9p-830L) >= real.mant_dig-10); + + // --- Problematic cases --- + // This is a situation where the series expansion fails to converge + assert( isNaN(betaIncompleteInv(0.12167, 4.0640301659679627772e19L, 0.0813601))); + // This next result is almost certainly erroneous. + // Mathematica states: "(cannot be determined by current methods)" + assert(betaIncomplete(1.16251e20, 2.18e39, 5.45e-20) == -real.infinity); + // WolframAlpha gives no result for this, though indicates that it approximately 1.0 - 1.3e-9 + assert(1 - betaIncomplete(0.01, 328222, 4.0375e-5) == 0x1.5f62926b4p-30); + } } @@ -1418,7 +1426,10 @@ assert(gammaIncompleteComplInv(3, 0)==real.infinity); // Fixed a bug that caused gammaIncompleteCompl to return a wrong value when // x was larger than a, but not by much, and both were large: // The value is from WolframAlpha (Gamma[100000, 100001, inf] / Gamma[100000]) -assert(fabs(gammaIncompleteCompl(100000, 100001) - 0.49831792109) < 0.000000000005); +static if (real.mant_dig >= 64) // incl. 80-bit reals + assert(fabs(gammaIncompleteCompl(100000, 100001) - 0.49831792109) < 0.000000000005); +else + assert(fabs(gammaIncompleteCompl(100000, 100001) - 0.49831792109) < 0.00000005); }