diff --git a/std/internal/math/errorfunction.d b/std/internal/math/errorfunction.d index 235e0e820..1708148e4 100644 --- a/std/internal/math/errorfunction.d +++ b/std/internal/math/errorfunction.d @@ -1,465 +1,457 @@ -/** - * Error Functions and Normal Distribution. - * - * License: $(WEB boost.org/LICENSE_1_0.txt, Boost License 1.0). - * Copyright: Based on the CEPHES math library, which is - * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com). - * Authors: Stephen L. Moshier, ported to D by Don Clugston - */ -/** - * Macros: - * NAN = $(RED NAN) - * SUP = $0 - * GAMMA = Γ - * INTEGRAL = ∫ - * INTEGRATE = $(BIG ∫$(SMALL $1)$2) - * POWER = $1$2 - * BIGSUM = $(BIG Σ $2$(SMALL $1)) - * CHOOSE = $(BIG () $(SMALL $1)$(SMALL $2) $(BIG )) - * TABLE_SV = - * - * $0
Special Values
- * SVH = $(TR $(TH $1) $(TH $2)) - * SV = $(TR $(TD $1) $(TD $2)) - */ -module std.internal.math.errorfunction; -import std.math; - -private { -immutable real EXP_2 = 0.13533528323661269189L; /* exp(-2) */ -enum real SQRT2PI = 2.50662827463100050242E0L; // sqrt(2pi) - - -enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max) -enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min*real.epsilon) = log(smallest denormal) -} - -T rationalPoly(T)(T x, const(T) [] numerator, const(T) [] denominator) pure nothrow -{ - return poly(x, numerator)/poly(x, denominator); -} - - - -private { - -/* erfc(x) = exp(-x^2) P(1/x)/Q(1/x) - 1/8 <= 1/x <= 1 - Peak relative error 5.8e-21 */ -immutable real[10] P = [ -0x1.30dfa809b3cc6676p-17, 0x1.38637cd0913c0288p+18, - 0x1.2f015e047b4476bp+22, 0x1.24726f46aa9ab08p+25, 0x1.64b13c6395dc9c26p+27, - 0x1.294c93046ad55b5p+29, 0x1.5962a82f92576dap+30, 0x1.11a709299faba04ap+31, - 0x1.11028065b087be46p+31, 0x1.0d8ef40735b097ep+30 -]; - -immutable real[11] Q = [ 0x1.14d8e2a72dec49f4p+19, 0x1.0c880ff467626e1p+23, - 0x1.04417ef060b58996p+26, 0x1.404e61ba86df4ebap+28, 0x1.0f81887bc82b873ap+30, - 0x1.4552a5e39fb49322p+31, 0x1.11779a0ceb2a01cep+32, 0x1.3544dd691b5b1d5cp+32, - 0x1.a91781f12251f02ep+31, 0x1.0d8ef3da605a1c86p+30, 1.0 -]; - - -/* erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2) - 1/128 <= 1/x < 1/8 - Peak relative error 1.9e-21 */ -immutable real[5] R = [ 0x1.b9f6d8b78e22459ep-6, 0x1.1b84686b0a4ea43ap-1, - 0x1.b8f6aebe96000c2ap+1, 0x1.cb1dbedac27c8ec2p+2, 0x1.cf885f8f572a4c14p+1 -]; - -immutable real[6] S = [ - 0x1.87ae3cae5f65eb5ep-5, 0x1.01616f266f306d08p+0, 0x1.a4abe0411eed6c22p+2, - 0x1.eac9ce3da600abaap+3, 0x1.5752a9ac2faebbccp+3, 1.0 -]; - -/* erf(x) = x P(x^2)/Q(x^2) - 0 <= x <= 1 - Peak relative error 7.6e-23 */ -immutable real[7] T = [ 0x1.0da01654d757888cp+20, 0x1.2eb7497bc8b4f4acp+17, - 0x1.79078c19530f72a8p+15, 0x1.4eaf2126c0b2c23p+11, 0x1.1f2ea81c9d272a2ep+8, - 0x1.59ca6e2d866e625p+2, 0x1.c188e0b67435faf4p-4 -]; - -immutable real[7] U = [ 0x1.dde6025c395ae34ep+19, 0x1.c4bc8b6235df35aap+18, - 0x1.8465900e88b6903ap+16, 0x1.855877093959ffdp+13, 0x1.e5c44395625ee358p+9, - 0x1.6a0fed103f1c68a6p+5, 1.0 -]; - -} - -/** - * Complementary error function - * - * erfc(x) = 1 - erf(x), and has high relative accuracy for - * values of x far from zero. (For values near zero, use erf(x)). - * - * 1 - erf(x) = 2/ $(SQRT)(π) - * $(INTEGRAL x, $(INFINITY)) exp( - $(POWER t, 2)) dt - * - * - * For small x, erfc(x) = 1 - erf(x); otherwise rational - * approximations are computed. - * - * A special function expx2(x) is used to suppress error amplification - * in computing exp(-x^2). - */ -real erfc(real a) -{ - if (a == real.infinity) - return 0.0; - if (a == -real.infinity) - return 2.0; - - real x; - - if (a < 0.0L ) - x = -a; - else - x = a; - if (x < 1.0) - return 1.0 - erf(a); - - real z = -a * a; - - if (z < -MAXLOG){ -// mtherr( "erfcl", UNDERFLOW ); - if (a < 0) return 2.0; - else return 0.0; - } - - /* Compute z = exp(z). */ - z = expx2(a, -1); - real y = 1.0/x; - - real p, q; - - if( x < 8.0 ) y = z * rationalPoly(y, P, Q); - else y = z * y * rationalPoly(y * y, R, S); - - if (a < 0.0L) - y = 2.0L - y; - - if (y == 0.0) { -// mtherr( "erfcl", UNDERFLOW ); - if (a < 0) return 2.0; - else return 0.0; - } - - return y; -} - - -private { -/* Exponentially scaled erfc function - exp(x^2) erfc(x) - valid for x > 1. - Use with normalDistribution and expx2. */ - -real erfce(real x) -{ - real p, q; - - real y = 1.0/x; - - if (x < 8.0) { - return rationalPoly( y, P, Q); - } else { - return y * rationalPoly(y*y, R, S); - } -} - -} - -/** - * Error function - * - * The integral is - * - * erf(x) = 2/ $(SQRT)(π) - * $(INTEGRAL 0, x) exp( - $(POWER t, 2)) dt - * - * The magnitude of x is limited to about 106.56 for IEEE 80-bit - * arithmetic; 1 or -1 is returned outside this range. - * - * For 0 <= |x| < 1, a rational polynomials are used; otherwise - * erf(x) = 1 - erfc(x). - * - * ACCURACY: - * Relative error: - * arithmetic domain # trials peak rms - * IEEE 0,1 50000 2.0e-19 5.7e-20 - */ -real erf(real x) -{ - if (x == 0.0) - return x; // deal with negative zero - if (x == -real.infinity) - return -1.0; - if (x == real.infinity) - return 1.0; - if (abs(x) > 1.0L) - return 1.0L - erfc(x); - - real z = x * x; - return x * rationalPoly(z, T, U); -} - -unittest { - // High resolution test points. - enum real erfc0_250 = 0.723663330078125 + 1.0279753638067014931732235184287934646022E-5; - enum real erfc0_375 = 0.5958709716796875 + 1.2118885490201676174914080878232469565953E-5; - enum real erfc0_500 = 0.4794921875 + 7.9346869534623172533461080354712635484242E-6; - enum real erfc0_625 = 0.3767547607421875 + 4.3570693945275513594941232097252997287766E-6; - enum real erfc0_750 = 0.2888336181640625 + 1.0748182422368401062165408589222625794046E-5; - enum real erfc0_875 = 0.215911865234375 + 1.3073705765341685464282101150637224028267E-5; - enum real erfc1_000 = 0.15728759765625 + 1.1609394035130658779364917390740703933002E-5; - enum real erfc1_125 = 0.111602783203125 + 8.9850951672359304215530728365232161564636E-6; - - enum real erf0_875 = (1-0.215911865234375) - 1.3073705765341685464282101150637224028267E-5; - - - assert(feqrel(erfc(0.250L), erfc0_250 )>=real.mant_dig-1); - assert(feqrel(erfc(0.375L), erfc0_375 )>=real.mant_dig-0); - assert(feqrel(erfc(0.500L), erfc0_500 )>=real.mant_dig-1); - assert(feqrel(erfc(0.625L), erfc0_625 )>=real.mant_dig-1); - assert(feqrel(erfc(0.750L), erfc0_750 )>=real.mant_dig-1); - assert(feqrel(erfc(0.875L), erfc0_875 )>=real.mant_dig-4); - assert(feqrel(erfc(1.000L), erfc1_000 )>=real.mant_dig-0); - assert(feqrel(erfc(1.125L), erfc1_125 )>=real.mant_dig-2); - assert(feqrel(erf(0.875L), erf0_875 )>=real.mant_dig-1); - // The DMC implementation of erfc() fails this next test (just) - assert(feqrel(erfc(4.1L),0.67000276540848983727e-8L)>=real.mant_dig-4); - - assert(isIdentical(erf(0.0),0.0)); - assert(isIdentical(erf(-0.0),-0.0)); - assert(erf(real.infinity) == 1.0); - assert(erf(-real.infinity) == -1.0); - assert(isIdentical(erf(NaN(0xDEF)),NaN(0xDEF))); - assert(isIdentical(erfc(NaN(0xDEF)),NaN(0xDEF))); - assert(isIdentical(erfc(real.infinity),0.0)); - assert(erfc(-real.infinity) == 2.0); - assert(erfc(0) == 1.0); -} - -/* - * Exponential of squared argument - * - * Computes y = exp(x*x) while suppressing error amplification - * that would ordinarily arise from the inexactness of the - * exponential argument x*x. - * - * If sign < 0, the result is inverted; i.e., y = exp(-x*x) . - * - * ACCURACY: - * Relative error: - * arithmetic domain # trials peak rms - * IEEE -106.566, 106.566 10^5 1.6e-19 4.4e-20 - */ - -real expx2(real x, int sign) -{ - /* - Cephes Math Library Release 2.9: June, 2000 - Copyright 2000 by Stephen L. Moshier - */ - const real M = 32768.0; - const real MINV = 3.0517578125e-5L; - - x = abs(x); - if (sign < 0) - x = -x; - - /* Represent x as an exact multiple of M plus a residual. - M is a power of 2 chosen so that exp(m * m) does not overflow - or underflow and so that |x - m| is small. */ - real m = MINV * floor(M * x + 0.5L); - real f = x - m; - - /* x^2 = m^2 + 2mf + f^2 */ - real u = m * m; - real u1 = 2 * m * f + f * f; - - if (sign < 0) { - u = -u; - u1 = -u1; - } - - if ((u+u1) > MAXLOG) - return real.infinity; - - /* u is exact, u1 is small. */ - return exp(u) * exp(u1); -} - - -/* -Computes the normal distribution function. - -The normal (or Gaussian, or bell-shaped) distribution is -defined as: - -normalDist(x) = 1/$(SQRT) π $(INTEGRAL -$(INFINITY), x) exp( - $(POWER t, 2)/2) dt - = 0.5 + 0.5 * erf(x/sqrt(2)) - = 0.5 * erfc(- x/sqrt(2)) - -To maintain accuracy at high values of x, use -normalDistribution(x) = 1 - normalDistribution(-x). - -Accuracy: -Within a few bits of machine resolution over the entire -range. - -References: -$(LINK http://www.netlib.org/cephes/ldoubdoc.html), -G. Marsaglia, "Evaluating the Normal Distribution", -Journal of Statistical Software 11, (July 2004). -*/ -real normalDistributionImpl(real a) -{ - real x = a * SQRT1_2; - real z = abs(x); - - if( z < 1.0 ) - return 0.5L + 0.5L * erf(x); - else { - real y = 0.5L * erfce(z); - /* Multiply by exp(-x^2 / 2) */ - z = expx2(a, -1); - y = y * sqrt(z); - if( x > 0.0L ) - y = 1.0L - y; - return y; - } -} - -unittest { -assert(fabs(normalDistributionImpl(1L) - (0.841344746068543))< 0.0000000000000005); -assert(isIdentical(normalDistributionImpl(NaN(0x325)), NaN(0x325))); -} - -/* - * Inverse of Normal distribution function - * - * Returns the argument, x, for which the area under the - * Normal probability density function (integrated from - * minus infinity to x) is equal to p. - * - * For small arguments 0 < p < exp(-2), the program computes - * z = sqrt( -2 log(p) ); then the approximation is - * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z) . - * For larger arguments, x/sqrt(2 pi) = w + w^3 R(w^2)/S(w^2)) , - * where w = p - 0.5 . - */ -real normalDistributionInvImpl(real p) -in { - assert(p>=0.0L && p<=1.0L, "Domain error"); -} -body -{ -static immutable real[8] P0 = -[ -0x1.758f4d969484bfdcp-7, 0x1.53cee17a59259dd2p-3, - -0x1.ea01e4400a9427a2p-1, 0x1.61f7504a0105341ap+1, -0x1.09475a594d0399f6p+2, - 0x1.7c59e7a0df99e3e2p+1, -0x1.87a81da52edcdf14p-1, 0x1.1fb149fd3f83600cp-7 -]; - -static immutable real[8] Q0 = -[ -0x1.64b92ae791e64bb2p-7, 0x1.7585c7d597298286p-3, - -0x1.40011be4f7591ce6p+0, 0x1.1fc067d8430a425ep+2, -0x1.21008ffb1e7ccdf2p+3, - 0x1.3d1581cf9bc12fccp+3, -0x1.53723a89fd8f083cp+2, 1.0 -]; - -static immutable real[10] P1 = -[ 0x1.20ceea49ea142f12p-13, 0x1.cbe8a7267aea80bp-7, - 0x1.79fea765aa787c48p-2, 0x1.d1f59faa1f4c4864p+1, 0x1.1c22e426a013bb96p+4, - 0x1.a8675a0c51ef3202p+5, 0x1.75782c4f83614164p+6, 0x1.7a2f3d90948f1666p+6, - 0x1.5cd116ee4c088c3ap+5, 0x1.1361e3eb6e3cc20ap+2 -]; - -static immutable real[10] Q1 = -[ 0x1.3a4ce1406cea98fap-13, 0x1.f45332623335cda2p-7, - 0x1.98f28bbd4b98db1p-2, 0x1.ec3b24f9c698091cp+1, 0x1.1cc56ecda7cf58e4p+4, - 0x1.92c6f7376bf8c058p+5, 0x1.4154c25aa47519b4p+6, 0x1.1b321d3b927849eap+6, - 0x1.403a5f5a4ce7b202p+4, 1.0 -]; - -static immutable real[8] P2 = -[ 0x1.8c124a850116a6d8p-21, 0x1.534abda3c2fb90bap-13, - 0x1.29a055ec93a4718cp-7, 0x1.6468e98aad6dd474p-3, 0x1.3dab2ef4c67a601cp+0, - 0x1.e1fb3a1e70c67464p+1, 0x1.b6cce8035ff57b02p+2, 0x1.9f4c9e749ff35f62p+1 -]; - -static immutable real[8] Q2 = -[ 0x1.af03f4fc0655e006p-21, 0x1.713192048d11fb2p-13, - 0x1.4357e5bbf5fef536p-7, 0x1.7fdac8749985d43cp-3, 0x1.4a080c813a2d8e84p+0, - 0x1.c3a4b423cdb41bdap+1, 0x1.8160694e24b5557ap+2, 1.0 -]; - -static immutable real[8] P3 = -[ -0x1.55da447ae3806168p-34, -0x1.145635641f8778a6p-24, - -0x1.abf46d6b48040128p-17, -0x1.7da550945da790fcp-11, -0x1.aa0b2a31157775fap-8, - 0x1.b11d97522eed26bcp-3, 0x1.1106d22f9ae89238p+1, 0x1.029a358e1e630f64p+1 -]; - -static immutable real[8] Q3 = -[ -0x1.74022dd5523e6f84p-34, -0x1.2cb60d61e29ee836p-24, - -0x1.d19e6ec03a85e556p-17, -0x1.9ea2a7b4422f6502p-11, -0x1.c54b1e852f107162p-8, - 0x1.e05268dd3c07989ep-3, 0x1.239c6aff14afbf82p+1, 1.0 -]; - - if(p<=0.0L || p>=1.0L) - { - if (p == 0.0L) - return -real.infinity; - if( p == 1.0L ) - return real.infinity; - return real.nan; // domain error - } - int code = 1; - real y = p; - if( y > (1.0L - EXP_2) ) { - y = 1.0L - y; - code = 0; - } - - real x, z, y2, x0, x1; - - if ( y > EXP_2 ) { - y = y - 0.5L; - y2 = y * y; - x = y + y * (y2 * rationalPoly( y2, P0, Q0)); - return x * SQRT2PI; - } - - x = sqrt( -2.0L * log(y) ); - x0 = x - log(x)/x; - z = 1.0L/x; - if ( x < 8.0L ) { - x1 = z * rationalPoly( z, P1, Q1); - } else if( x < 32.0L ) { - x1 = z * rationalPoly( z, P2, Q2); - } else { - x1 = z * rationalPoly( z, P3, Q3); - } - x = x0 - x1; - if ( code != 0 ) { - x = -x; - } - return x; -} - - -unittest { - // TODO: Use verified test points. - // The values below are from Excel 2003. - assert(fabs(normalDistributionInvImpl(0.001) - (-3.09023230616779))< 0.00000000000005); - assert(fabs(normalDistributionInvImpl(1e-50) - (-14.9333375347885))< 0.00000000000005); - assert(feqrel(normalDistributionInvImpl(0.999), -normalDistributionInvImpl(0.001)) > real.mant_dig-6); - - // Excel 2003 gets all the following values wrong! - assert(normalDistributionInvImpl(0.0) == -real.infinity); - assert(normalDistributionInvImpl(1.0) == real.infinity); - assert(normalDistributionInvImpl(0.5) == 0); - // (Excel 2003 returns norminv(p) = -30 for all p < 1e-200). - // The value tested here is the one the function returned in Jan 2006. - real unknown1 = normalDistributionInvImpl(1e-250L); - assert( fabs(unknown1 -(-33.79958617269L) ) < 0.00000005); -} +/** + * Error Functions and Normal Distribution. + * + * License: $(WEB boost.org/LICENSE_1_0.txt, Boost License 1.0). + * Copyright: Based on the CEPHES math library, which is + * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com). + * Authors: Stephen L. Moshier, ported to D by Don Clugston + */ +/** + * Macros: + * NAN = $(RED NAN) + * SUP = $0 + * GAMMA = Γ + * INTEGRAL = ∫ + * INTEGRATE = $(BIG ∫$(SMALL $1)$2) + * POWER = $1$2 + * BIGSUM = $(BIG Σ $2$(SMALL $1)) + * CHOOSE = $(BIG () $(SMALL $1)$(SMALL $2) $(BIG )) + * TABLE_SV = + * + * $0
Special Values
+ * SVH = $(TR $(TH $1) $(TH $2)) + * SV = $(TR $(TD $1) $(TD $2)) + */ +module std.internal.math.errorfunction; +import std.math; + +private { +immutable real EXP_2 = 0.13533528323661269189L; /* exp(-2) */ +enum real SQRT2PI = 2.50662827463100050242E0L; // sqrt(2pi) + + +enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max) +enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min*real.epsilon) = log(smallest denormal) +} + +T rationalPoly(T)(T x, const(T) [] numerator, const(T) [] denominator) pure nothrow +{ + return poly(x, numerator)/poly(x, denominator); +} + + + +private { + +/* erfc(x) = exp(-x^2) P(1/x)/Q(1/x) + 1/8 <= 1/x <= 1 + Peak relative error 5.8e-21 */ +immutable real [] P = [ -0x1.30dfa809b3cc6676p-17, 0x1.38637cd0913c0288p+18, + 0x1.2f015e047b4476bp+22, 0x1.24726f46aa9ab08p+25, 0x1.64b13c6395dc9c26p+27, + 0x1.294c93046ad55b5p+29, 0x1.5962a82f92576dap+30, 0x1.11a709299faba04ap+31, + 0x1.11028065b087be46p+31, 0x1.0d8ef40735b097ep+30 +]; + +immutable real [] Q = [ 0x1.14d8e2a72dec49f4p+19, 0x1.0c880ff467626e1p+23, + 0x1.04417ef060b58996p+26, 0x1.404e61ba86df4ebap+28, 0x1.0f81887bc82b873ap+30, + 0x1.4552a5e39fb49322p+31, 0x1.11779a0ceb2a01cep+32, 0x1.3544dd691b5b1d5cp+32, + 0x1.a91781f12251f02ep+31, 0x1.0d8ef3da605a1c86p+30, 1.0 +]; + + +/* erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2) + 1/128 <= 1/x < 1/8 + Peak relative error 1.9e-21 */ +immutable real [] R = [ 0x1.b9f6d8b78e22459ep-6, 0x1.1b84686b0a4ea43ap-1, + 0x1.b8f6aebe96000c2ap+1, 0x1.cb1dbedac27c8ec2p+2, 0x1.cf885f8f572a4c14p+1 +]; + +immutable real [] S = [ + 0x1.87ae3cae5f65eb5ep-5, 0x1.01616f266f306d08p+0, 0x1.a4abe0411eed6c22p+2, + 0x1.eac9ce3da600abaap+3, 0x1.5752a9ac2faebbccp+3, 1.0 +]; + +/* erf(x) = x P(x^2)/Q(x^2) + 0 <= x <= 1 + Peak relative error 7.6e-23 */ +immutable real [] T = [ 0x1.0da01654d757888cp+20, 0x1.2eb7497bc8b4f4acp+17, + 0x1.79078c19530f72a8p+15, 0x1.4eaf2126c0b2c23p+11, 0x1.1f2ea81c9d272a2ep+8, + 0x1.59ca6e2d866e625p+2, 0x1.c188e0b67435faf4p-4 +]; + +immutable real [] U = [ 0x1.dde6025c395ae34ep+19, 0x1.c4bc8b6235df35aap+18, + 0x1.8465900e88b6903ap+16, 0x1.855877093959ffdp+13, 0x1.e5c44395625ee358p+9, + 0x1.6a0fed103f1c68a6p+5, 1.0 +]; + +} + +/** + * Complementary error function + * + * erfc(x) = 1 - erf(x), and has high relative accuracy for + * values of x far from zero. (For values near zero, use erf(x)). + * + * 1 - erf(x) = 2/ $(SQRT)(π) + * $(INTEGRAL x, $(INFINITY)) exp( - $(POWER t, 2)) dt + * + * + * For small x, erfc(x) = 1 - erf(x); otherwise rational + * approximations are computed. + * + * A special function expx2(x) is used to suppress error amplification + * in computing exp(-x^2). + */ +real erfc(real a) +{ + if (a == real.infinity) + return 0.0; + if (a == -real.infinity) + return 2.0; + + real x; + + if (a < 0.0L ) + x = -a; + else + x = a; + if (x < 1.0) + return 1.0 - erf(a); + + real z = -a * a; + + if (z < -MAXLOG){ +// mtherr( "erfcl", UNDERFLOW ); + if (a < 0) return 2.0; + else return 0.0; + } + + /* Compute z = exp(z). */ + z = expx2(a, -1); + real y = 1.0/x; + + real p, q; + + if( x < 8.0 ) y = z * rationalPoly(y, P, Q); + else y = z * y * rationalPoly(y * y, R, S); + + if (a < 0.0L) + y = 2.0L - y; + + if (y == 0.0) { +// mtherr( "erfcl", UNDERFLOW ); + if (a < 0) return 2.0; + else return 0.0; + } + + return y; +} + + +private { +/* Exponentially scaled erfc function + exp(x^2) erfc(x) + valid for x > 1. + Use with normalDistribution and expx2. */ + +real erfce(real x) +{ + real p, q; + + real y = 1.0/x; + + if (x < 8.0) { + return rationalPoly( y, P, Q); + } else { + return y * rationalPoly(y*y, R, S); + } +} + +} + +/** + * Error function + * + * The integral is + * + * erf(x) = 2/ $(SQRT)(π) + * $(INTEGRAL 0, x) exp( - $(POWER t, 2)) dt + * + * The magnitude of x is limited to about 106.56 for IEEE 80-bit + * arithmetic; 1 or -1 is returned outside this range. + * + * For 0 <= |x| < 1, a rational polynomials are used; otherwise + * erf(x) = 1 - erfc(x). + * + * ACCURACY: + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,1 50000 2.0e-19 5.7e-20 + */ +real erf(real x) +{ + if (x == 0.0) + return x; // deal with negative zero + if (x == -real.infinity) + return -1.0; + if (x == real.infinity) + return 1.0; + if (abs(x) > 1.0L) + return 1.0L - erfc(x); + + real z = x * x; + return x * rationalPoly(z, T, U); +} + +unittest { + // High resolution test points. + enum real erfc0_250 = 0.723663330078125 + 1.0279753638067014931732235184287934646022E-5; + enum real erfc0_375 = 0.5958709716796875 + 1.2118885490201676174914080878232469565953E-5; + enum real erfc0_500 = 0.4794921875 + 7.9346869534623172533461080354712635484242E-6; + enum real erfc0_625 = 0.3767547607421875 + 4.3570693945275513594941232097252997287766E-6; + enum real erfc0_750 = 0.2888336181640625 + 1.0748182422368401062165408589222625794046E-5; + enum real erfc0_875 = 0.215911865234375 + 1.3073705765341685464282101150637224028267E-5; + enum real erfc1_000 = 0.15728759765625 + 1.1609394035130658779364917390740703933002E-5; + enum real erfc1_125 = 0.111602783203125 + 8.9850951672359304215530728365232161564636E-6; + + enum real erf0_875 = (1-0.215911865234375) - 1.3073705765341685464282101150637224028267E-5; + + + assert(feqrel(erfc(0.250L), erfc0_250 )>=real.mant_dig-1); + assert(feqrel(erfc(0.375L), erfc0_375 )>=real.mant_dig-0); + assert(feqrel(erfc(0.500L), erfc0_500 )>=real.mant_dig-1); + assert(feqrel(erfc(0.625L), erfc0_625 )>=real.mant_dig-1); + assert(feqrel(erfc(0.750L), erfc0_750 )>=real.mant_dig-1); + assert(feqrel(erfc(0.875L), erfc0_875 )>=real.mant_dig-4); + assert(feqrel(erfc(1.000L), erfc1_000 )>=real.mant_dig-0); + assert(feqrel(erfc(1.125L), erfc1_125 )>=real.mant_dig-2); + assert(feqrel(erf(0.875L), erf0_875 )>=real.mant_dig-1); + // The DMC implementation of erfc() fails this next test (just) + assert(feqrel(erfc(4.1L),0.67000276540848983727e-8L)>=real.mant_dig-4); + + assert(isIdentical(erf(0.0),0.0)); + assert(isIdentical(erf(-0.0),-0.0)); + assert(erf(real.infinity) == 1.0); + assert(erf(-real.infinity) == -1.0); + assert(isIdentical(erf(NaN(0xDEF)),NaN(0xDEF))); + assert(isIdentical(erfc(NaN(0xDEF)),NaN(0xDEF))); + assert(isIdentical(erfc(real.infinity),0.0)); + assert(erfc(-real.infinity) == 2.0); + assert(erfc(0) == 1.0); +} + +/* + * Exponential of squared argument + * + * Computes y = exp(x*x) while suppressing error amplification + * that would ordinarily arise from the inexactness of the + * exponential argument x*x. + * + * If sign < 0, the result is inverted; i.e., y = exp(-x*x) . + * + * ACCURACY: + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -106.566, 106.566 10^5 1.6e-19 4.4e-20 + */ + +real expx2(real x, int sign) +{ + /* + Cephes Math Library Release 2.9: June, 2000 + Copyright 2000 by Stephen L. Moshier + */ + const real M = 32768.0; + const real MINV = 3.0517578125e-5L; + + x = abs(x); + if (sign < 0) + x = -x; + + /* Represent x as an exact multiple of M plus a residual. + M is a power of 2 chosen so that exp(m * m) does not overflow + or underflow and so that |x - m| is small. */ + real m = MINV * floor(M * x + 0.5L); + real f = x - m; + + /* x^2 = m^2 + 2mf + f^2 */ + real u = m * m; + real u1 = 2 * m * f + f * f; + + if (sign < 0) { + u = -u; + u1 = -u1; + } + + if ((u+u1) > MAXLOG) + return real.infinity; + + /* u is exact, u1 is small. */ + return exp(u) * exp(u1); +} + + +/* +Computes the normal distribution function. + +The normal (or Gaussian, or bell-shaped) distribution is +defined as: + +normalDist(x) = 1/$(SQRT) π $(INTEGRAL -$(INFINITY), x) exp( - $(POWER t, 2)/2) dt + = 0.5 + 0.5 * erf(x/sqrt(2)) + = 0.5 * erfc(- x/sqrt(2)) + +To maintain accuracy at high values of x, use +normalDistribution(x) = 1 - normalDistribution(-x). + +Accuracy: +Within a few bits of machine resolution over the entire +range. + +References: +$(LINK http://www.netlib.org/cephes/ldoubdoc.html), +G. Marsaglia, "Evaluating the Normal Distribution", +Journal of Statistical Software 11, (July 2004). +*/ +real normalDistributionImpl(real a) +{ + real x = a * SQRT1_2; + real z = abs(x); + + if( z < 1.0 ) + return 0.5L + 0.5L * erf(x); + else { + real y = 0.5L * erfce(z); + /* Multiply by exp(-x^2 / 2) */ + z = expx2(a, -1); + y = y * sqrt(z); + if( x > 0.0L ) + y = 1.0L - y; + return y; + } +} + +unittest { +assert(fabs(normalDistributionImpl(1L) - (0.841344746068543))< 0.0000000000000005); +assert(isIdentical(normalDistributionImpl(NaN(0x325)), NaN(0x325))); +} + +/* + * Inverse of Normal distribution function + * + * Returns the argument, x, for which the area under the + * Normal probability density function (integrated from + * minus infinity to x) is equal to p. + * + * For small arguments 0 < p < exp(-2), the program computes + * z = sqrt( -2 log(p) ); then the approximation is + * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z) . + * For larger arguments, x/sqrt(2 pi) = w + w^3 R(w^2)/S(w^2)) , + * where w = p - 0.5 . + */ +real normalDistributionInvImpl(real p) +in { + assert(p>=0.0L && p<=1.0L, "Domain error"); +} +body +{ +immutable real[] P0 = [ -0x1.758f4d969484bfdcp-7, 0x1.53cee17a59259dd2p-3, + -0x1.ea01e4400a9427a2p-1, 0x1.61f7504a0105341ap+1, -0x1.09475a594d0399f6p+2, + 0x1.7c59e7a0df99e3e2p+1, -0x1.87a81da52edcdf14p-1, 0x1.1fb149fd3f83600cp-7 +]; + +immutable real[] Q0 = [ -0x1.64b92ae791e64bb2p-7, 0x1.7585c7d597298286p-3, + -0x1.40011be4f7591ce6p+0, 0x1.1fc067d8430a425ep+2, -0x1.21008ffb1e7ccdf2p+3, + 0x1.3d1581cf9bc12fccp+3, -0x1.53723a89fd8f083cp+2, 1.0 +]; + +immutable real[] P1 = [ 0x1.20ceea49ea142f12p-13, 0x1.cbe8a7267aea80bp-7, + 0x1.79fea765aa787c48p-2, 0x1.d1f59faa1f4c4864p+1, 0x1.1c22e426a013bb96p+4, + 0x1.a8675a0c51ef3202p+5, 0x1.75782c4f83614164p+6, 0x1.7a2f3d90948f1666p+6, + 0x1.5cd116ee4c088c3ap+5, 0x1.1361e3eb6e3cc20ap+2 +]; + +immutable real[] Q1 = [ 0x1.3a4ce1406cea98fap-13, 0x1.f45332623335cda2p-7, + 0x1.98f28bbd4b98db1p-2, 0x1.ec3b24f9c698091cp+1, 0x1.1cc56ecda7cf58e4p+4, + 0x1.92c6f7376bf8c058p+5, 0x1.4154c25aa47519b4p+6, 0x1.1b321d3b927849eap+6, + 0x1.403a5f5a4ce7b202p+4, 1.0 +]; + +immutable real[] P2 = [ 0x1.8c124a850116a6d8p-21, 0x1.534abda3c2fb90bap-13, + 0x1.29a055ec93a4718cp-7, 0x1.6468e98aad6dd474p-3, 0x1.3dab2ef4c67a601cp+0, + 0x1.e1fb3a1e70c67464p+1, 0x1.b6cce8035ff57b02p+2, 0x1.9f4c9e749ff35f62p+1 +]; + +immutable real[] Q2 = [ 0x1.af03f4fc0655e006p-21, 0x1.713192048d11fb2p-13, + 0x1.4357e5bbf5fef536p-7, 0x1.7fdac8749985d43cp-3, 0x1.4a080c813a2d8e84p+0, + 0x1.c3a4b423cdb41bdap+1, 0x1.8160694e24b5557ap+2, 1.0 +]; + +immutable real[] P3 = [ -0x1.55da447ae3806168p-34, -0x1.145635641f8778a6p-24, + -0x1.abf46d6b48040128p-17, -0x1.7da550945da790fcp-11, -0x1.aa0b2a31157775fap-8, + 0x1.b11d97522eed26bcp-3, 0x1.1106d22f9ae89238p+1, 0x1.029a358e1e630f64p+1 +]; + +immutable real[] Q3 = [ -0x1.74022dd5523e6f84p-34, -0x1.2cb60d61e29ee836p-24, + -0x1.d19e6ec03a85e556p-17, -0x1.9ea2a7b4422f6502p-11, -0x1.c54b1e852f107162p-8, + 0x1.e05268dd3c07989ep-3, 0x1.239c6aff14afbf82p+1, 1.0 +]; + + if(p<=0.0L || p>=1.0L) + { + if (p == 0.0L) + return -real.infinity; + if( p == 1.0L ) + return real.infinity; + return real.nan; // domain error + } + int code = 1; + real y = p; + if( y > (1.0L - EXP_2) ) { + y = 1.0L - y; + code = 0; + } + + real x, z, y2, x0, x1; + + if ( y > EXP_2 ) { + y = y - 0.5L; + y2 = y * y; + x = y + y * (y2 * rationalPoly( y2, P0, Q0)); + return x * SQRT2PI; + } + + x = sqrt( -2.0L * log(y) ); + x0 = x - log(x)/x; + z = 1.0L/x; + if ( x < 8.0L ) { + x1 = z * rationalPoly( z, P1, Q1); + } else if( x < 32.0L ) { + x1 = z * rationalPoly( z, P2, Q2); + } else { + x1 = z * rationalPoly( z, P3, Q3); + } + x = x0 - x1; + if ( code != 0 ) { + x = -x; + } + return x; +} + + +unittest { + // TODO: Use verified test points. + // The values below are from Excel 2003. + assert(fabs(normalDistributionInvImpl(0.001) - (-3.09023230616779))< 0.00000000000005); + assert(fabs(normalDistributionInvImpl(1e-50) - (-14.9333375347885))< 0.00000000000005); + assert(feqrel(normalDistributionInvImpl(0.999), -normalDistributionInvImpl(0.001)) > real.mant_dig-6); + + // Excel 2003 gets all the following values wrong! + assert(normalDistributionInvImpl(0.0) == -real.infinity); + assert(normalDistributionInvImpl(1.0) == real.infinity); + assert(normalDistributionInvImpl(0.5) == 0); + // (Excel 2003 returns norminv(p) = -30 for all p < 1e-200). + // The value tested here is the one the function returned in Jan 2006. + real unknown1 = normalDistributionInvImpl(1e-250L); + assert( fabs(unknown1 -(-33.79958617269L) ) < 0.00000005); +} diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d index 160fd59a6..cce6eaa85 100644 --- a/std/internal/math/gammafunction.d +++ b/std/internal/math/gammafunction.d @@ -1,1404 +1,1404 @@ -/** - * Implementation of the gamma and beta functions, and their integrals. - * - * License: $(WEB boost.org/LICENSE_1_0.txt, Boost License 1.0). - * Copyright: Based on the CEPHES math library, which is - * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com). - * Authors: Stephen L. Moshier (original C code). Conversion to D by Don Clugston - * - * -Macros: - * TABLE_SV = - * - * $0
Special Values
- * SVH = $(TR $(TH $1) $(TH $2)) - * SV = $(TR $(TD $1) $(TD $2)) - * GAMMA = Γ - * INTEGRATE = $(BIG ∫$(SMALL $1)$2) - * POWER = $1$2 - * NAN = $(RED NAN) - */ -module std.internal.math.gammafunction; -import std.internal.math.errorfunction; -import std.math; - -private { - -enum real SQRT2PI = 2.50662827463100050242E0L; // sqrt(2pi) -immutable real EULERGAMMA = 0.57721_56649_01532_86060_65120_90082_40243_10421_59335_93992L; /** Euler-Mascheroni constant 0.57721566.. */ - -// Polynomial approximations for gamma and loggamma. - -immutable real[8] GammaNumeratorCoeffs = [ 1.0, - 0x1.acf42d903366539ep-1, 0x1.73a991c8475f1aeap-2, 0x1.c7e918751d6b2a92p-4, - 0x1.86d162cca32cfe86p-6, 0x1.0c378e2e6eaf7cd8p-8, 0x1.dc5c66b7d05feb54p-12, - 0x1.616457b47e448694p-15 -]; - -immutable real[9] GammaDenominatorCoeffs = [ 1.0, - 0x1.a8f9faae5d8fc8bp-2, -0x1.cb7895a6756eebdep-3, -0x1.7b9bab006d30652ap-5, - 0x1.c671af78f312082ep-6, -0x1.a11ebbfaf96252dcp-11, -0x1.447b4d2230a77ddap-10, - 0x1.ec1d45bb85e06696p-13,-0x1.d4ce24d05bd0a8e6p-17 -]; - -immutable real[9] GammaSmallCoeffs = [ 1.0, - 0x1.2788cfc6fb618f52p-1, -0x1.4fcf4026afa2f7ecp-1, -0x1.5815e8fa24d7e306p-5, - 0x1.5512320aea2ad71ap-3, -0x1.59af0fb9d82e216p-5, -0x1.3b4b61d3bfdf244ap-7, - 0x1.d9358e9d9d69fd34p-8, -0x1.38fc4bcbada775d6p-10 -]; - -immutable real[9] GammaSmallNegCoeffs = [ -1.0, - 0x1.2788cfc6fb618f54p-1, 0x1.4fcf4026afa2bc4cp-1, -0x1.5815e8fa2468fec8p-5, - -0x1.5512320baedaf4b6p-3, -0x1.59af0fa283baf07ep-5, 0x1.3b4a70de31e05942p-7, - 0x1.d9398be3bad13136p-8, 0x1.291b73ee05bcbba2p-10 -]; - -immutable real[7] logGammaStirlingCoeffs = [ - 0x1.5555555555553f98p-4, -0x1.6c16c16c07509b1p-9, 0x1.a01a012461cbf1e4p-11, - -0x1.3813089d3f9d164p-11, 0x1.b911a92555a277b8p-11, -0x1.ed0a7b4206087b22p-10, - 0x1.402523859811b308p-8 -]; - -immutable real[7] logGammaNumerator = [ - -0x1.0edd25913aaa40a2p+23, -0x1.31c6ce2e58842d1ep+24, -0x1.f015814039477c3p+23, - -0x1.74ffe40c4b184b34p+22, -0x1.0d9c6d08f9eab55p+20, -0x1.54c6b71935f1fc88p+16, - -0x1.0e761b42932b2aaep+11 -]; - -immutable real[8] logGammaDenominator = [ - -0x1.4055572d75d08c56p+24, -0x1.deeb6013998e4d76p+24, -0x1.106f7cded5dcc79ep+24, - -0x1.25e17184848c66d2p+22, -0x1.301303b99a614a0ap+19, -0x1.09e76ab41ae965p+15, - -0x1.00f95ced9e5f54eep+9, 1.0 -]; - -/* - * Helper function: Gamma function computed by Stirling's formula. - * - * Stirling's formula for the gamma function is: - * - * $(GAMMA)(x) = sqrt(2 π) xx-0.5 exp(-x) (1 + 1/x P(1/x)) - * - */ -real gammaStirling(real x) -{ - // CEPHES code Copyright 1994 by Stephen L. Moshier - - static immutable real[9] SmallStirlingCoeffs = [ - 0x1.55555555555543aap-4, 0x1.c71c71c720dd8792p-9, -0x1.5f7268f0b5907438p-9, - -0x1.e13cd410e0477de6p-13, 0x1.9b0f31643442616ep-11, 0x1.2527623a3472ae08p-14, - -0x1.37f6bc8ef8b374dep-11,-0x1.8c968886052b872ap-16, 0x1.76baa9c6d3eeddbcp-11 - ]; - - static immutable real[9] LargeStirlingCoeffs = [ 1.0L, - 8.33333333333333333333E-2L, 3.47222222222222222222E-3L, - -2.68132716049382716049E-3L, -2.29472093621399176955E-4L, - 7.84039221720066627474E-4L, 6.97281375836585777429E-5L - ]; - - real w = 1.0L/x; - real y = exp(x); - if ( x > 1024.0L ) { - // For large x, use rational coefficients from the analytical expansion. - w = poly(w, LargeStirlingCoeffs); - // Avoid overflow in pow() - real v = pow( x, 0.5L * x - 0.25L ); - y = v * (v / y); - } - else { - w = 1.0L + w * poly( w, SmallStirlingCoeffs); - y = pow( x, x - 0.5L ) / y; - } - y = SQRT2PI * y * w; - return y; -} - -} // private - -public: -/// The maximum value of x for which gamma(x) < real.infinity. -enum real MAXGAMMA = 1755.5483429L; - - -/***************************************************** - * The Gamma function, $(GAMMA)(x) - * - * $(GAMMA)(x) is a generalisation of the factorial function - * to real and complex numbers. - * Like x!, $(GAMMA)(x+1) = x*$(GAMMA)(x). - * - * Mathematically, if z.re > 0 then - * $(GAMMA)(z) = $(INTEGRATE 0, ∞) $(POWER t, z-1)$(POWER e, -t) dt - * - * $(TABLE_SV - * $(SVH x, $(GAMMA)(x) ) - * $(SV $(NAN), $(NAN) ) - * $(SV ±0.0, ±∞) - * $(SV integer > 0, (x-1)! ) - * $(SV integer < 0, $(NAN) ) - * $(SV +∞, +∞ ) - * $(SV -∞, $(NAN) ) - * ) - */ -real gamma(real x) -{ -/* Based on code from the CEPHES library. - * CEPHES code Copyright 1994 by Stephen L. Moshier - * - * Arguments |x| <= 13 are reduced by recurrence and the function - * approximated by a rational function of degree 7/8 in the - * interval (2,3). Large arguments are handled by Stirling's - * formula. Large negative arguments are made positive using - * a reflection formula. - */ - - real q, z; - if (isNaN(x)) return x; - if (x == -x.infinity) return real.nan; - if ( fabs(x) > MAXGAMMA ) return real.infinity; - if (x==0) return 1.0 / x; // +- infinity depending on sign of x, create an exception. - - q = fabs(x); - - if ( q > 13.0L ) { - // Large arguments are handled by Stirling's - // formula. Large negative arguments are made positive using - // the reflection formula. - - if ( x < 0.0L ) { - if (x < -1/real.epsilon) - { - // Large negatives lose all precision - return real.nan; - } - int sgngam = 1; // sign of gamma. - long intpart = cast(long)(q); - if (q == intpart) - return real.nan; // poles for all integers <0. - real p = intpart; - if ( (intpart & 1) == 0 ) - sgngam = -1; - z = q - p; - if ( z > 0.5L ) { - p += 1.0L; - z = q - p; - } - z = q * sin( PI * z ); - z = fabs(z) * gammaStirling(q); - if ( z <= PI/real.max ) return sgngam * real.infinity; - return sgngam * PI/z; - } else { - return gammaStirling(x); - } - } - - // Arguments |x| <= 13 are reduced by recurrence and the function - // approximated by a rational function of degree 7/8 in the - // interval (2,3). - - z = 1.0L; - while ( x >= 3.0L ) { - x -= 1.0L; - z *= x; - } - - while ( x < -0.03125L ) { - z /= x; - x += 1.0L; - } - - if ( x <= 0.03125L ) { - if ( x == 0.0L ) - return real.nan; - else { - if ( x < 0.0L ) { - x = -x; - return z / (x * poly( x, GammaSmallNegCoeffs )); - } else { - return z / (x * poly( x, GammaSmallCoeffs )); - } - } - } - - while ( x < 2.0L ) { - z /= x; - x += 1.0L; - } - if ( x == 2.0L ) return z; - - x -= 2.0L; - return z * poly( x, GammaNumeratorCoeffs ) / poly( x, GammaDenominatorCoeffs ); -} - -unittest { - // gamma(n) = factorial(n-1) if n is an integer. - real fact = 1.0L; - for (int i=1; fact real.mant_dig-15); - fact *= (i*1.0L); - } - assert(gamma(0.0) == real.infinity); - assert(gamma(-0.0) == -real.infinity); - assert(isNaN(gamma(-1.0))); - assert(isNaN(gamma(-15.0))); - assert(isIdentical(gamma(NaN(0xABC)), NaN(0xABC))); - assert(gamma(real.infinity) == real.infinity); - assert(gamma(real.max) == real.infinity); - assert(isNaN(gamma(-real.infinity))); - assert(gamma(real.min*real.epsilon) == real.infinity); - assert(gamma(MAXGAMMA)< real.infinity); - assert(gamma(MAXGAMMA*2) == real.infinity); - - // Test some high-precision values (50 decimal digits) - real SQRT_PI = 1.77245385090551602729816748334114518279754945612238L; - - assert(feqrel(gamma(0.5L), SQRT_PI) == real.mant_dig); - - assert(feqrel(gamma(1.0 / 3.L), 2.67893853470774763365569294097467764412868937795730L) >= real.mant_dig-2); - assert(feqrel(gamma(0.25L), - 3.62560990822190831193068515586767200299516768288006L) >= real.mant_dig-1); - assert(feqrel(gamma(1.0 / 5.0L), - 4.59084371199880305320475827592915200343410999829340L) >= real.mant_dig-1); -} - -/***************************************************** - * Natural logarithm of gamma function. - * - * Returns the base e (2.718...) logarithm of the absolute - * value of the gamma function of the argument. - * - * For reals, logGamma is equivalent to log(fabs(gamma(x))). - * - * $(TABLE_SV - * $(SVH x, logGamma(x) ) - * $(SV $(NAN), $(NAN) ) - * $(SV integer <= 0, +∞ ) - * $(SV ±∞, +∞ ) - * ) - */ -real logGamma(real x) -{ - /* Based on code from the CEPHES library. - * CEPHES code Copyright 1994 by Stephen L. Moshier - * - * For arguments greater than 33, the logarithm of the gamma - * function is approximated by the logarithmic version of - * Stirling's formula using a polynomial approximation of - * degree 4. Arguments between -33 and +33 are reduced by - * recurrence to the interval [2,3] of a rational approximation. - * The cosecant reflection formula is employed for arguments - * less than -33. - */ - real q, w, z, f, nx; - - if (isNaN(x)) return x; - if (fabs(x) == x.infinity) return x.infinity; - - if( x < -34.0L ) - { - q = -x; - w = logGamma(q); - real p = floor(q); - if ( p == q ) - return real.infinity; - int intpart = cast(int)(p); - real sgngam = 1; - if ( (intpart & 1) == 0 ) - sgngam = -1; - z = q - p; - if ( z > 0.5L ) { - p += 1.0L; - z = p - q; - } - z = q * sin( PI * z ); - if ( z == 0.0L ) - return sgngam * real.infinity; - /* z = LOGPI - logl( z ) - w; */ - z = log( PI/z ) - w; - return z; - } - - if( x < 13.0L ) - { - z = 1.0L; - nx = floor( x + 0.5L ); - f = x - nx; - while ( x >= 3.0L ) { - nx -= 1.0L; - x = nx + f; - z *= x; - } - while ( x < 2.0L ) { - if( fabs(x) <= 0.03125 ) - { - if ( x == 0.0L ) - return real.infinity; - if ( x < 0.0L ) - { - x = -x; - q = z / (x * poly( x, GammaSmallNegCoeffs)); - } else - q = z / (x * poly( x, GammaSmallCoeffs)); - return log( fabs(q) ); - } - z /= nx + f; - nx += 1.0L; - x = nx + f; - } - z = fabs(z); - if ( x == 2.0L ) - return log(z); - x = (nx - 2.0L) + f; - real p = x * rationalPoly( x, logGammaNumerator, logGammaDenominator); - return log(z) + p; - } - - // const real MAXLGM = 1.04848146839019521116e+4928L; - // if( x > MAXLGM ) return sgngaml * real.infinity; - - const real LOGSQRT2PI = 0.91893853320467274178L; // log( sqrt( 2*pi ) ) - - q = ( x - 0.5L ) * log(x) - x + LOGSQRT2PI; - if (x > 1.0e10L) return q; - real p = 1.0L / (x*x); - q += poly( p, logGammaStirlingCoeffs ) / x; - return q ; -} - -unittest { - assert(isIdentical(logGamma(NaN(0xDEF)), NaN(0xDEF))); - assert(logGamma(real.infinity) == real.infinity); - assert(logGamma(-1.0) == real.infinity); - assert(logGamma(0.0) == real.infinity); - assert(logGamma(-50.0) == real.infinity); - assert(isIdentical(0.0L, logGamma(1.0L))); - assert(isIdentical(0.0L, logGamma(2.0L))); - assert(logGamma(real.min*real.epsilon) == real.infinity); - assert(logGamma(-real.min*real.epsilon) == real.infinity); - - // x, correct loggamma(x), correct d/dx loggamma(x). - static real[] testpoints = [ - 8.0L, 8.525146484375L + 1.48766904143001655310E-5, 2.01564147795560999654E0L, - 8.99993896484375e-1L, 6.6375732421875e-2L + 5.11505711292524166220E-6L, -7.54938684259372234258E-1, - 7.31597900390625e-1L, 2.2369384765625e-1 + 5.21506341809849792422E-6L, -1.13355566660398608343E0L, - 2.31639862060546875e-1L, 1.3686676025390625L + 1.12609441752996145670E-5L, -4.56670961813812679012E0, - 1.73162841796875L, -8.88214111328125e-2L + 3.36207740803753034508E-6L, 2.33339034686200586920E-1L, - 1.23162841796875L, -9.3902587890625e-2L + 1.28765089229009648104E-5L, -2.49677345775751390414E-1L, - 7.3786976294838206464e19L, 3.301798506038663053312e21L - 1.656137564136932662487046269677E5L, - 4.57477139169563904215E1L, - 1.08420217248550443401E-19L, 4.36682586669921875e1L + 1.37082843669932230418E-5L, - -9.22337203685477580858E18L, - 1.0L, 0.0L, -5.77215664901532860607E-1L, - 2.0L, 0.0L, 4.22784335098467139393E-1L, - -0.5L, 1.2655029296875L + 9.19379714539648894580E-6L, 3.64899739785765205590E-2L, - -1.5L, 8.6004638671875e-1L + 6.28657731014510932682E-7L, 7.03156640645243187226E-1L, - -2.5L, -5.6243896484375E-2L + 1.79986700949327405470E-7, 1.10315664064524318723E0L, - -3.5L, -1.30902099609375L + 1.43111007079536392848E-5L, 1.38887092635952890151E0L - ]; - // TODO: test derivatives as well. - for (int i=0; i real.mant_dig-5); - if (testpoints[i] real.mant_dig-5); - } - } - 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); -} - - -private { -enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max) -enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min*real.epsilon) = log(smallest denormal) -enum real BETA_BIG = 9.223372036854775808e18L; -enum real BETA_BIGINV = 1.084202172485504434007e-19L; -} - -/** Incomplete beta integral - * - * Returns incomplete beta integral of the arguments, evaluated - * from zero to x. The regularized incomplete beta function is defined as - * - * betaIncomplete(a, b, x) = Γ(a+b)/(Γ(a) Γ(b)) * - * $(INTEGRATE 0, x) $(POWER t, a-1)$(POWER (1-t),b-1) dt - * - * and is the same as the the cumulative distribution function. - * - * The domain of definition is 0 <= x <= 1. In this - * implementation a and b are restricted to positive values. - * The integral from x to 1 may be obtained by the symmetry - * relation - * - * betaIncompleteCompl(a, b, x ) = betaIncomplete( b, a, 1-x ) - * - * The integral is evaluated by a continued fraction expansion - * or, when b*x is small, by a power series. - */ -real betaIncomplete(real aa, real bb, real xx ) -{ - if ( !(aa>0 && bb>0) ) - { - if ( isNaN(aa) ) return aa; - if ( isNaN(bb) ) return bb; - return real.nan; // domain error - } - if (!(xx>0 && xx<1.0)) { - if (isNaN(xx)) return xx; - if ( xx == 0.0L ) return 0.0; - if ( xx == 1.0L ) return 1.0; - return real.nan; // domain error - } - if ( (bb * xx) <= 1.0L && xx <= 0.95L) { - return betaDistPowerSeries(aa, bb, xx); - } - real x; - real xc; // = 1 - x - - real a, b; - int flag = 0; - - /* Reverse a and b if x is greater than the mean. */ - if( xx > (aa/(aa+bb)) ) { - // here x > aa/(aa+bb) and (bb*x>1 or x>0.95) - flag = 1; - a = bb; - b = aa; - xc = xx; - x = 1.0L - xx; - } else { - a = aa; - b = bb; - xc = 1.0L - xx; - x = xx; - } - - if( flag == 1 && (b * x) <= 1.0L && x <= 0.95L) { - // here xx > aa/(aa+bb) and ((bb*xx>1) or xx>0.95) and (aa*(1-xx)<=1) and xx > 0.05 - return 1.0 - betaDistPowerSeries(a, b, x); // note loss of precision - } - - real w; - // Choose expansion for optimal convergence - // One is for x * (a+b+2) < (a+1), - // the other is for x * (a+b+2) > (a+1). - real y = x * (a+b-2.0L) - (a-1.0L); - if( y < 0.0L ) { - w = betaDistExpansion1( a, b, x ); - } else { - w = betaDistExpansion2( a, b, x ) / xc; - } - - /* Multiply w by the factor - a b - x (1-x) Gamma(a+b) / ( a Gamma(a) Gamma(b) ) . */ - - y = a * log(x); - real t = b * log(xc); - if ( (a+b) < MAXGAMMA && fabs(y) < MAXLOG && fabs(t) < MAXLOG ) { - t = pow(xc,b); - t *= pow(x,a); - t /= a; - t *= w; - t *= gamma(a+b) / (gamma(a) * gamma(b)); - } else { - /* Resort to logarithms. */ - y += t + logGamma(a+b) - logGamma(a) - logGamma(b); - y += log(w/a); - - t = exp(y); -/+ - // There seems to be a bug in Cephes at this point. - // Problems occur for y > MAXLOG, not y < MINLOG. - if( y < MINLOG ) { - t = 0.0L; - } else { - t = exp(y); - } -+/ - } - if( flag == 1 ) { -/+ // CEPHES includes this code, but I think it is erroneous. - if( t <= real.epsilon ) { - t = 1.0L - real.epsilon; - } else -+/ - t = 1.0L - t; - } - return t; -} - -/** Inverse of incomplete beta integral - * - * Given y, the function finds x such that - * - * betaIncomplete(a, b, x) == y - * - * Newton iterations or interval halving is used. - */ -real betaIncompleteInv(real aa, real bb, real yy0 ) -{ - real a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt; - int i, rflg, dir, nflg; - - if (isNaN(yy0)) return yy0; - if (isNaN(aa)) return aa; - if (isNaN(bb)) return bb; - if( yy0 <= 0.0L ) - return 0.0L; - if( yy0 >= 1.0L ) - return 1.0L; - x0 = 0.0L; - yl = 0.0L; - x1 = 1.0L; - yh = 1.0L; - if( aa <= 1.0L || bb <= 1.0L ) { - dithresh = 1.0e-7L; - rflg = 0; - a = aa; - b = bb; - y0 = yy0; - x = a/(a+b); - y = betaIncomplete( a, b, x ); - nflg = 0; - goto ihalve; - } else { - nflg = 0; - dithresh = 1.0e-4L; - } - - // approximation to inverse function - - yp = -normalDistributionInvImpl( yy0 ); - - if( yy0 > 0.5L ) { - rflg = 1; - a = bb; - b = aa; - y0 = 1.0L - yy0; - yp = -yp; - } else { - rflg = 0; - a = aa; - b = bb; - y0 = yy0; - } - - lgm = (yp * yp - 3.0L)/6.0L; - x = 2.0L/( 1.0L/(2.0L * a-1.0L) + 1.0L/(2.0L * b - 1.0L) ); - d = yp * sqrt( x + lgm ) / x - - ( 1.0L/(2.0L * b - 1.0L) - 1.0L/(2.0L * a - 1.0L) ) - * (lgm + (5.0L/6.0L) - 2.0L/(3.0L * x)); - d = 2.0L * d; - if( d < MINLOG ) { - x = 1.0L; - goto under; - } - x = a/( a + b * exp(d) ); - y = betaIncomplete( a, b, x ); - yp = (y - y0)/y0; - if( fabs(yp) < 0.2 ) - goto newt; - - /* Resort to interval halving if not close enough. */ -ihalve: - - dir = 0; - di = 0.5L; - for( i=0; i<400; i++ ) { - if( i != 0 ) { - x = x0 + di * (x1 - x0); - if( x == 1.0L ) { - x = 1.0L - real.epsilon; - } - if( x == 0.0L ) { - di = 0.5; - x = x0 + di * (x1 - x0); - if( x == 0.0 ) - goto under; - } - y = betaIncomplete( a, b, x ); - yp = (x1 - x0)/(x1 + x0); - if( fabs(yp) < dithresh ) - goto newt; - yp = (y-y0)/y0; - if( fabs(yp) < dithresh ) - goto newt; - } - if( y < y0 ) { - x0 = x; - yl = y; - if( dir < 0 ) { - dir = 0; - di = 0.5L; - } else if( dir > 3 ) - di = 1.0L - (1.0L - di) * (1.0L - di); - else if( dir > 1 ) - di = 0.5L * di + 0.5L; - else - di = (y0 - y)/(yh - yl); - dir += 1; - if( x0 > 0.95L ) { - if( rflg == 1 ) { - rflg = 0; - a = aa; - b = bb; - y0 = yy0; - } else { - rflg = 1; - a = bb; - b = aa; - y0 = 1.0 - yy0; - } - x = 1.0L - x; - y = betaIncomplete( a, b, x ); - x0 = 0.0; - yl = 0.0; - x1 = 1.0; - yh = 1.0; - goto ihalve; - } - } else { - x1 = x; - if( rflg == 1 && x1 < real.epsilon ) { - x = 0.0L; - goto done; - } - yh = y; - if( dir > 0 ) { - dir = 0; - di = 0.5L; - } - else if( dir < -3 ) - di = di * di; - else if( dir < -1 ) - di = 0.5L * di; - else - di = (y - y0)/(yh - yl); - dir -= 1; - } - } - if( x0 >= 1.0L ) { - // partial loss of precision - x = 1.0L - real.epsilon; - goto done; - } - if( x <= 0.0L ) { -under: - // underflow has occurred - x = real.min * real.min; - goto done; - } - -newt: - - if ( nflg ) { - goto done; - } - nflg = 1; - lgm = logGamma(a+b) - logGamma(a) - logGamma(b); - - for( i=0; i<15; i++ ) { - /* Compute the function at this point. */ - if ( i != 0 ) - y = betaIncomplete(a,b,x); - if ( y < yl ) { - x = x0; - y = yl; - } else if( y > yh ) { - x = x1; - y = yh; - } else if( y < y0 ) { - x0 = x; - yl = y; - } else { - x1 = x; - yh = y; - } - if( x == 1.0L || x == 0.0L ) - break; - /* Compute the derivative of the function at this point. */ - d = (a - 1.0L) * log(x) + (b - 1.0L) * log(1.0L - x) + lgm; - if ( d < MINLOG ) { - goto done; - } - if ( d > MAXLOG ) { - break; - } - d = exp(d); - /* Compute the step to the next approximation of x. */ - d = (y - y0)/d; - xt = x - d; - if ( xt <= x0 ) { - y = (x - x0) / (x1 - x0); - xt = x0 + 0.5L * y * (x - x0); - if( xt <= 0.0L ) - break; - } - if ( xt >= x1 ) { - y = (x1 - x) / (x1 - x0); - xt = x1 - 0.5L * y * (x1 - x); - if ( xt >= 1.0L ) - break; - } - x = xt; - if ( fabs(d/x) < (128.0L * real.epsilon) ) - goto done; - } - /* Did not converge. */ - dithresh = 256.0L * real.epsilon; - goto ihalve; - -done: - if ( rflg ) { - if( x <= real.epsilon ) - x = 1.0L - real.epsilon; - else - x = 1.0L - x; - } - return x; -} - -unittest { // also tested by the normal distribution - // check NaN propagation - assert(isIdentical(betaIncomplete(NaN(0xABC),2,3), NaN(0xABC))); - assert(isIdentical(betaIncomplete(7,NaN(0xABC),3), NaN(0xABC))); - assert(isIdentical(betaIncomplete(7,15,NaN(0xABC)), NaN(0xABC))); - assert(isIdentical(betaIncompleteInv(NaN(0xABC),1,17), NaN(0xABC))); - assert(isIdentical(betaIncompleteInv(2,NaN(0xABC),8), NaN(0xABC))); - assert(isIdentical(betaIncompleteInv(2,3, NaN(0xABC)), NaN(0xABC))); - - assert(isNaN(betaIncomplete(-1, 2, 3))); - - assert(betaIncomplete(1, 2, 0)==0); - assert(betaIncomplete(1, 2, 1)==1); - assert(isNaN(betaIncomplete(1, 2, 3))); - assert(betaIncompleteInv(1, 1, 0)==0); - assert(betaIncompleteInv(1, 1, 1)==1); - - // Test some values against Microsoft Excel 2003. - - assert(fabs(betaIncomplete(8, 10, 0.2) - 0.010_934_315_236_957_2L) < 0.000_000_000_5); - assert(fabs(betaIncomplete(2, 2.5, 0.9) - 0.989_722_597_604_107L) < 0.000_000_000_000_5); - assert(fabs(betaIncomplete(1000, 800, 0.5) - 1.17914088832798E-06L) < 0.000_000_05e-6); - - assert(fabs(betaIncomplete(0.0001, 10000, 0.0001) - 0.999978059369989L) < 0.000_000_000_05); - - assert(fabs(betaIncompleteInv(5, 10, 0.2) - 0.229121208190918L) < 0.000_000_5L); - assert(fabs(betaIncompleteInv(4, 7, 0.8) - 0.483657360076904L) < 0.000_000_5L); - - // Coverage tests. I don't have correct values for these tests, but - // these values cover most of the code, so they are useful for - // regression testing. - // 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. - -// Excel 2003 gives clearly erroneous results (betadist>1) when a and x are tiny and b is huge. -// The correct results are for these next tests are unknown. - -// real testpoint1 = betaIncomplete(1e-10, 5e20, 8e-21); -// assert(testpoint1 == 0x1.ffff_ffff_c906_404cp-1L); - - assert(betaIncomplete(0.01, 327726.7, 0.545113) == 1.0); - assert(betaIncompleteInv(0.01, 8e-48, 5.45464e-20)==1-real.epsilon); - assert(betaIncompleteInv(0.01, 8e-48, 9e-26)==1-real.epsilon); - - assert(betaIncomplete(0.01, 498.437, 0.0121433) == 0x1.ffff_8f72_19197402p-1); - assert(1- betaIncomplete(0.01, 328222, 4.0375e-5) == 0x1.5f62926b4p-30); - - // 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); - assert(betaIncompleteInv(0x1.ff1275ae5b939bcap-41, 4.6713e18, 0.0813601)==0x1.f97749d90c7adba8p-63); - 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. - assert(betaIncomplete(1.16251e20, 2.18e39, 5.45e-20)==-real.infinity); -} - - -private { -// Implementation functions - -// Continued fraction expansion #1 for incomplete beta integral -// Use when x < (a+1)/(a+b+2) -real betaDistExpansion1(real a, real b, real x ) -{ - real xk, pk, pkm1, pkm2, qk, qkm1, qkm2; - real k1, k2, k3, k4, k5, k6, k7, k8; - real r, t, ans; - int n; - - k1 = a; - k2 = a + b; - k3 = a; - k4 = a + 1.0L; - k5 = 1.0L; - k6 = b - 1.0L; - k7 = k4; - k8 = a + 2.0L; - - pkm2 = 0.0L; - qkm2 = 1.0L; - pkm1 = 1.0L; - qkm1 = 1.0L; - ans = 1.0L; - r = 1.0L; - n = 0; - const real thresh = 3.0L * real.epsilon; - do { - xk = -( x * k1 * k2 )/( k3 * k4 ); - pk = pkm1 + pkm2 * xk; - qk = qkm1 + qkm2 * xk; - pkm2 = pkm1; - pkm1 = pk; - qkm2 = qkm1; - qkm1 = qk; - - xk = ( x * k5 * k6 )/( k7 * k8 ); - pk = pkm1 + pkm2 * xk; - qk = qkm1 + qkm2 * xk; - pkm2 = pkm1; - pkm1 = pk; - qkm2 = qkm1; - qkm1 = qk; - - if( qk != 0.0L ) - r = pk/qk; - if( r != 0.0L ) { - t = fabs( (ans - r)/r ); - ans = r; - } else { - t = 1.0L; - } - - if( t < thresh ) - return ans; - - k1 += 1.0L; - k2 += 1.0L; - k3 += 2.0L; - k4 += 2.0L; - k5 += 1.0L; - k6 -= 1.0L; - k7 += 2.0L; - k8 += 2.0L; - - if( (fabs(qk) + fabs(pk)) > BETA_BIG ) { - pkm2 *= BETA_BIGINV; - pkm1 *= BETA_BIGINV; - qkm2 *= BETA_BIGINV; - qkm1 *= BETA_BIGINV; - } - if( (fabs(qk) < BETA_BIGINV) || (fabs(pk) < BETA_BIGINV) ) { - pkm2 *= BETA_BIG; - pkm1 *= BETA_BIG; - qkm2 *= BETA_BIG; - qkm1 *= BETA_BIG; - } - } - while( ++n < 400 ); -// loss of precision has occurred -// mtherr( "incbetl", PLOSS ); - return ans; -} - -// Continued fraction expansion #2 for incomplete beta integral -// Use when x > (a+1)/(a+b+2) -real betaDistExpansion2(real a, real b, real x ) -{ - real xk, pk, pkm1, pkm2, qk, qkm1, qkm2; - real k1, k2, k3, k4, k5, k6, k7, k8; - real r, t, ans, z; - - k1 = a; - k2 = b - 1.0L; - k3 = a; - k4 = a + 1.0L; - k5 = 1.0L; - k6 = a + b; - k7 = a + 1.0L; - k8 = a + 2.0L; - - pkm2 = 0.0L; - qkm2 = 1.0L; - pkm1 = 1.0L; - qkm1 = 1.0L; - z = x / (1.0L-x); - ans = 1.0L; - r = 1.0L; - int n = 0; - const real thresh = 3.0L * real.epsilon; - do { - - xk = -( z * k1 * k2 )/( k3 * k4 ); - pk = pkm1 + pkm2 * xk; - qk = qkm1 + qkm2 * xk; - pkm2 = pkm1; - pkm1 = pk; - qkm2 = qkm1; - qkm1 = qk; - - xk = ( z * k5 * k6 )/( k7 * k8 ); - pk = pkm1 + pkm2 * xk; - qk = qkm1 + qkm2 * xk; - pkm2 = pkm1; - pkm1 = pk; - qkm2 = qkm1; - qkm1 = qk; - - if( qk != 0.0L ) - r = pk/qk; - if( r != 0.0L ) { - t = fabs( (ans - r)/r ); - ans = r; - } else - t = 1.0L; - - if( t < thresh ) - return ans; - k1 += 1.0L; - k2 -= 1.0L; - k3 += 2.0L; - k4 += 2.0L; - k5 += 1.0L; - k6 += 1.0L; - k7 += 2.0L; - k8 += 2.0L; - - if( (fabs(qk) + fabs(pk)) > BETA_BIG ) { - pkm2 *= BETA_BIGINV; - pkm1 *= BETA_BIGINV; - qkm2 *= BETA_BIGINV; - qkm1 *= BETA_BIGINV; - } - if( (fabs(qk) < BETA_BIGINV) || (fabs(pk) < BETA_BIGINV) ) { - pkm2 *= BETA_BIG; - pkm1 *= BETA_BIG; - qkm2 *= BETA_BIG; - qkm1 *= BETA_BIG; - } - } while( ++n < 400 ); -// loss of precision has occurred -//mtherr( "incbetl", PLOSS ); - return ans; -} - -/* Power series for incomplete gamma integral. - Use when b*x is small. */ -real betaDistPowerSeries(real a, real b, real x ) -{ - real ai = 1.0L / a; - real u = (1.0L - b) * x; - real v = u / (a + 1.0L); - real t1 = v; - real t = u; - real n = 2.0L; - real s = 0.0L; - real z = real.epsilon * ai; - while( fabs(v) > z ) { - u = (n - b) * x / n; - t *= u; - v = t / (a + n); - s += v; - n += 1.0L; - } - s += t1; - s += ai; - - u = a * log(x); - if ( (a+b) < MAXGAMMA && fabs(u) < MAXLOG ) { - t = gamma(a+b)/(gamma(a)*gamma(b)); - s = s * t * pow(x,a); - } else { - t = logGamma(a+b) - logGamma(a) - logGamma(b) + u + log(s); - - if( t < MINLOG ) { - s = 0.0L; - } else - s = exp(t); - } - return s; -} - -} - -/*************************************** - * Incomplete gamma integral and its complement - * - * These functions are defined by - * - * gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a) - * - * gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x) - * = ($(INTEGRATE x, ∞) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a) - * - * In this implementation both arguments must be positive. - * The integral is evaluated by either a power series or - * continued fraction expansion, depending on the relative - * values of a and x. - */ -real gammaIncomplete(real a, real x ) -in { - assert(x >= 0); - assert(a > 0); -} -body { - /* left tail of incomplete gamma function: - * - * inf. k - * a -x - x - * x e > ---------- - * - - - * k=0 | (a+k+1) - * - */ - if (x==0) - return 0.0L; - - if ( (x > 1.0L) && (x > a ) ) - return 1.0L - gammaIncompleteCompl(a,x); - - real ax = a * log(x) - x - logGamma(a); -/+ - if( ax < MINLOGL ) return 0; // underflow - // { mtherr( "igaml", UNDERFLOW ); return( 0.0L ); } -+/ - ax = exp(ax); - - /* power series */ - real r = a; - real c = 1.0L; - real ans = 1.0L; - - do { - r += 1.0L; - c *= x/r; - ans += c; - } while( c/ans > real.epsilon ); - - return ans * ax/a; -} - -/** ditto */ -real gammaIncompleteCompl(real a, real x ) -in { - assert(x >= 0); - assert(a > 0); -} -body { - if (x==0) - return 1.0L; - if ( (x < 1.0L) || (x < a) ) - return 1.0L - gammaIncomplete(a,x); - - // DAC (Cephes bug fix): This is necessary to avoid - // spurious nans, eg - // log(x)-x = NaN when x = real.infinity - const real MAXLOGL = 1.1356523406294143949492E4L; - if (x > MAXLOGL) return 0; // underflow - - real ax = a * log(x) - x - logGamma(a); -//const real MINLOGL = -1.1355137111933024058873E4L; -// if ( ax < MINLOGL ) return 0; // underflow; - ax = exp(ax); - - - /* continued fraction */ - real y = 1.0L - a; - real z = x + y + 1.0L; - real c = 0.0L; - - real pk, qk, t; - - real pkm2 = 1.0L; - real qkm2 = x; - real pkm1 = x + 1.0L; - real qkm1 = z * x; - real ans = pkm1/qkm1; - - do { - c += 1.0L; - y += 1.0L; - z += 2.0L; - real yc = y * c; - pk = pkm1 * z - pkm2 * yc; - qk = qkm1 * z - qkm2 * yc; - if( qk != 0.0L ) { - real r = pk/qk; - t = fabs( (ans - r)/r ); - ans = r; - } else { - t = 1.0L; - } - pkm2 = pkm1; - pkm1 = pk; - qkm2 = qkm1; - qkm1 = qk; - - const real BIG = 9.223372036854775808e18L; - - if ( fabs(pk) > BIG ) { - pkm2 /= BIG; - pkm1 /= BIG; - qkm2 /= BIG; - qkm1 /= BIG; - } - } while ( t > real.epsilon ); - - return ans * ax; -} - -/** Inverse of complemented incomplete gamma integral - * - * Given a and y, the function finds x such that - * - * gammaIncompleteCompl( a, x ) = p. - * - * Starting with the approximate value x = a $(POWER t, 3), where - * t = 1 - d - normalDistributionInv(p) sqrt(d), - * and d = 1/9a, - * the routine performs up to 10 Newton iterations to find the - * root of incompleteGammaCompl(a,x) - p = 0. - */ -real gammaIncompleteComplInv(real a, real p) -in { - assert(p>=0 && p<= 1); - assert(a>0); -} -body { - if (p==0) return real.infinity; - - real y0 = p; - const real MAXLOGL = 1.1356523406294143949492E4L; - real x0, x1, x, yl, yh, y, d, lgm, dithresh; - int i, dir; - - /* bound the solution */ - x0 = real.max; - yl = 0.0L; - x1 = 0.0L; - yh = 1.0L; - dithresh = 4.0 * real.epsilon; - - /* approximation to inverse function */ - d = 1.0L/(9.0L*a); - y = 1.0L - d - normalDistributionInvImpl(y0) * sqrt(d); - x = a * y * y * y; - - lgm = logGamma(a); - - for( i=0; i<10; i++ ) { - if( x > x0 || x < x1 ) - goto ihalve; - y = gammaIncompleteCompl(a,x); - if ( y < yl || y > yh ) - goto ihalve; - if ( y < y0 ) { - x0 = x; - yl = y; - } else { - x1 = x; - yh = y; - } - /* compute the derivative of the function at this point */ - d = (a - 1.0L) * log(x0) - x0 - lgm; - if ( d < -MAXLOGL ) - goto ihalve; - d = -exp(d); - /* compute the step to the next approximation of x */ - d = (y - y0)/d; - x = x - d; - if ( i < 3 ) continue; - if ( fabs(d/x) < dithresh ) return x; - } - - /* Resort to interval halving if Newton iteration did not converge. */ -ihalve: - d = 0.0625L; - if ( x0 == real.max ) { - if( x <= 0.0L ) - x = 1.0L; - while( x0 == real.max ) { - x = (1.0L + d) * x; - y = gammaIncompleteCompl( a, x ); - if ( y < y0 ) { - x0 = x; - yl = y; - break; - } - d = d + d; - } - } - d = 0.5L; - dir = 0; - - for( i=0; i<400; i++ ) { - x = x1 + d * (x0 - x1); - y = gammaIncompleteCompl( a, x ); - lgm = (x0 - x1)/(x1 + x0); - if ( fabs(lgm) < dithresh ) - break; - lgm = (y - y0)/y0; - if ( fabs(lgm) < dithresh ) - break; - if ( x <= 0.0L ) - break; - if ( y > y0 ) { - x1 = x; - yh = y; - if ( dir < 0 ) { - dir = 0; - d = 0.5L; - } else if ( dir > 1 ) - d = 0.5L * d + 0.5L; - else - d = (y0 - yl)/(yh - yl); - dir += 1; - } else { - x0 = x; - yl = y; - if ( dir > 0 ) { - dir = 0; - d = 0.5L; - } else if ( dir < -1 ) - d = 0.5L * d; - else - d = (y0 - yl)/(yh - yl); - dir -= 1; - } - } - /+ - if( x == 0.0L ) - mtherr( "igamil", UNDERFLOW ); - +/ - return x; -} - -unittest { -//Values from Excel's GammaInv(1-p, x, 1) -assert(fabs(gammaIncompleteComplInv(1, 0.5) - 0.693147188044814) < 0.00000005); -assert(fabs(gammaIncompleteComplInv(12, 0.99) - 5.42818075054289) < 0.00000005); -assert(fabs(gammaIncompleteComplInv(100, 0.8) - 91.5013985848288L) < 0.000005); - -assert(gammaIncomplete(1, 0)==0); -assert(gammaIncompleteCompl(1, 0)==1); -assert(gammaIncomplete(4545, real.infinity)==1); - -// Values from Excel's (1-GammaDist(x, alpha, 1, TRUE)) - -assert(fabs(1.0L-gammaIncompleteCompl(0.5, 2) - 0.954499729507309L) < 0.00000005); -assert(fabs(gammaIncomplete(0.5, 2) - 0.954499729507309L) < 0.00000005); -// Fixed Cephes bug: -assert(gammaIncompleteCompl(384, real.infinity)==0); -assert(gammaIncompleteComplInv(3, 0)==real.infinity); -} - - -/** Digamma function -* -* The digamma function is the logarithmic derivative of the gamma function. -* -* digamma(x) = d/dx logGamma(x) -* -*/ -real digamma(real x) -{ - // Based on CEPHES, Stephen L. Moshier. - - // DAC: These values are Bn / n for n=2,4,6,8,10,12,14. - const real [] Bn_n = [ - 1.0L/(6*2), -1.0L/(30*4), 1.0L/(42*6), -1.0L/(30*8), - 5.0L/(66*10), -691.0L/(2730*12), 7.0L/(6*14) ]; - - real p, q, nz, s, w, y, z; - long i, n; - int negative; - - negative = 0; - nz = 0.0; - - if ( x <= 0.0 ) { - negative = 1; - q = x; - p = floor(q); - if( p == q ) { - return real.nan; // singularity. - } - /* Remove the zeros of tan(PI x) - * by subtracting the nearest integer from x - */ - nz = q - p; - if ( nz != 0.5 ) { - if ( nz > 0.5 ) { - p += 1.0; - nz = q - p; - } - nz = PI/tan(PI*nz); - } else { - nz = 0.0; - } - x = 1.0 - x; - } - - // check for small positive integer - if ((x <= 13.0) && (x == floor(x)) ) { - y = 0.0; - n = lrint(x); - // DAC: CEPHES bugfix. Cephes did this in reverse order, which - // created a larger roundoff error. - for (i=n-1; i>0; --i) { - y+=1.0L/i; - } - y -= EULERGAMMA; - goto done; - } - - s = x; - w = 0.0; - while ( s < 10.0 ) { - w += 1.0/s; - s += 1.0; - } - - if ( s < 1.0e17 ) { - z = 1.0/(s * s); - y = z * poly(z, Bn_n); - } else - y = 0.0; - - y = log(s) - 0.5L/s - y - w; - -done: - if ( negative ) { - y -= nz; - } - return y; -} - -version(unittest) import std.c.stdio; -unittest { - // Exact values - assert(digamma(1)== -EULERGAMMA); - assert(feqrel(digamma(0.25), -PI/2 - 3* LN2 - EULERGAMMA)>=real.mant_dig-7); - assert(feqrel(digamma(1.0L/6), -PI/2 *sqrt(3.0L) - 2* LN2 -1.5*log(3.0L) - EULERGAMMA)>=real.mant_dig-7); - assert(digamma(-5)!<>0); - assert(feqrel(digamma(2.5), -EULERGAMMA - 2*LN2 + 2.0 + 2.0L/3)>=real.mant_dig-9); - assert(isIdentical(digamma(NaN(0xABC)), NaN(0xABC))); - - for (int k=1; k<40; ++k) { - real y=0; - for (int u=k; u>=1; --u) { - y+= 1.0L/u; - } - assert(feqrel(digamma(k+1),-EULERGAMMA + y) >=real.mant_dig-2); - } -} +/** + * Implementation of the gamma and beta functions, and their integrals. + * + * License: $(WEB boost.org/LICENSE_1_0.txt, Boost License 1.0). + * Copyright: Based on the CEPHES math library, which is + * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com). + * Authors: Stephen L. Moshier (original C code). Conversion to D by Don Clugston + * + * +Macros: + * TABLE_SV = + * + * $0
Special Values
+ * SVH = $(TR $(TH $1) $(TH $2)) + * SV = $(TR $(TD $1) $(TD $2)) + * GAMMA = Γ + * INTEGRATE = $(BIG ∫$(SMALL $1)$2) + * POWER = $1$2 + * NAN = $(RED NAN) + */ +module std.internal.math.gammafunction; +import std.internal.math.errorfunction; +import std.math; + +private { + +enum real SQRT2PI = 2.50662827463100050242E0L; // sqrt(2pi) +immutable real EULERGAMMA = 0.57721_56649_01532_86060_65120_90082_40243_10421_59335_93992L; /** Euler-Mascheroni constant 0.57721566.. */ + +// Polynomial approximations for gamma and loggamma. + +immutable real[] GammaNumeratorCoeffs = [ 1.0, + 0x1.acf42d903366539ep-1, 0x1.73a991c8475f1aeap-2, 0x1.c7e918751d6b2a92p-4, + 0x1.86d162cca32cfe86p-6, 0x1.0c378e2e6eaf7cd8p-8, 0x1.dc5c66b7d05feb54p-12, + 0x1.616457b47e448694p-15 +]; + +immutable real[] GammaDenominatorCoeffs = [ 1.0, + 0x1.a8f9faae5d8fc8bp-2, -0x1.cb7895a6756eebdep-3, -0x1.7b9bab006d30652ap-5, + 0x1.c671af78f312082ep-6, -0x1.a11ebbfaf96252dcp-11, -0x1.447b4d2230a77ddap-10, + 0x1.ec1d45bb85e06696p-13,-0x1.d4ce24d05bd0a8e6p-17 +]; + +immutable real[] GammaSmallCoeffs = [ 1.0, + 0x1.2788cfc6fb618f52p-1, -0x1.4fcf4026afa2f7ecp-1, -0x1.5815e8fa24d7e306p-5, + 0x1.5512320aea2ad71ap-3, -0x1.59af0fb9d82e216p-5, -0x1.3b4b61d3bfdf244ap-7, + 0x1.d9358e9d9d69fd34p-8, -0x1.38fc4bcbada775d6p-10 +]; + +immutable real[] GammaSmallNegCoeffs = [ -1.0, + 0x1.2788cfc6fb618f54p-1, 0x1.4fcf4026afa2bc4cp-1, -0x1.5815e8fa2468fec8p-5, + -0x1.5512320baedaf4b6p-3, -0x1.59af0fa283baf07ep-5, 0x1.3b4a70de31e05942p-7, + 0x1.d9398be3bad13136p-8, 0x1.291b73ee05bcbba2p-10 +]; + +immutable real[] logGammaStirlingCoeffs = [ + 0x1.5555555555553f98p-4, -0x1.6c16c16c07509b1p-9, 0x1.a01a012461cbf1e4p-11, + -0x1.3813089d3f9d164p-11, 0x1.b911a92555a277b8p-11, -0x1.ed0a7b4206087b22p-10, + 0x1.402523859811b308p-8 +]; + +immutable real[] logGammaNumerator = [ + -0x1.0edd25913aaa40a2p+23, -0x1.31c6ce2e58842d1ep+24, -0x1.f015814039477c3p+23, + -0x1.74ffe40c4b184b34p+22, -0x1.0d9c6d08f9eab55p+20, -0x1.54c6b71935f1fc88p+16, + -0x1.0e761b42932b2aaep+11 +]; + +immutable real[] logGammaDenominator = [ + -0x1.4055572d75d08c56p+24, -0x1.deeb6013998e4d76p+24, -0x1.106f7cded5dcc79ep+24, + -0x1.25e17184848c66d2p+22, -0x1.301303b99a614a0ap+19, -0x1.09e76ab41ae965p+15, + -0x1.00f95ced9e5f54eep+9, 1.0 +]; + +/* + * Helper function: Gamma function computed by Stirling's formula. + * + * Stirling's formula for the gamma function is: + * + * $(GAMMA)(x) = sqrt(2 π) xx-0.5 exp(-x) (1 + 1/x P(1/x)) + * + */ +real gammaStirling(real x) +{ + // CEPHES code Copyright 1994 by Stephen L. Moshier + + const real SmallStirlingCoeffs[] = [ + 0x1.55555555555543aap-4, 0x1.c71c71c720dd8792p-9, -0x1.5f7268f0b5907438p-9, + -0x1.e13cd410e0477de6p-13, 0x1.9b0f31643442616ep-11, 0x1.2527623a3472ae08p-14, + -0x1.37f6bc8ef8b374dep-11,-0x1.8c968886052b872ap-16, 0x1.76baa9c6d3eeddbcp-11 + ]; + + const real LargeStirlingCoeffs[] = [ 1.0L, + 8.33333333333333333333E-2L, 3.47222222222222222222E-3L, + -2.68132716049382716049E-3L, -2.29472093621399176955E-4L, + 7.84039221720066627474E-4L, 6.97281375836585777429E-5L + ]; + + real w = 1.0L/x; + real y = exp(x); + if ( x > 1024.0L ) { + // For large x, use rational coefficients from the analytical expansion. + w = poly(w, LargeStirlingCoeffs); + // Avoid overflow in pow() + real v = pow( x, 0.5L * x - 0.25L ); + y = v * (v / y); + } + else { + w = 1.0L + w * poly( w, SmallStirlingCoeffs); + y = pow( x, x - 0.5L ) / y; + } + y = SQRT2PI * y * w; + return y; +} + +} // private + +public: +/// The maximum value of x for which gamma(x) < real.infinity. +enum real MAXGAMMA = 1755.5483429L; + + +/***************************************************** + * The Gamma function, $(GAMMA)(x) + * + * $(GAMMA)(x) is a generalisation of the factorial function + * to real and complex numbers. + * Like x!, $(GAMMA)(x+1) = x*$(GAMMA)(x). + * + * Mathematically, if z.re > 0 then + * $(GAMMA)(z) = $(INTEGRATE 0, ∞) $(POWER t, z-1)$(POWER e, -t) dt + * + * $(TABLE_SV + * $(SVH x, $(GAMMA)(x) ) + * $(SV $(NAN), $(NAN) ) + * $(SV ±0.0, ±∞) + * $(SV integer > 0, (x-1)! ) + * $(SV integer < 0, $(NAN) ) + * $(SV +∞, +∞ ) + * $(SV -∞, $(NAN) ) + * ) + */ +real gamma(real x) +{ +/* Based on code from the CEPHES library. + * CEPHES code Copyright 1994 by Stephen L. Moshier + * + * Arguments |x| <= 13 are reduced by recurrence and the function + * approximated by a rational function of degree 7/8 in the + * interval (2,3). Large arguments are handled by Stirling's + * formula. Large negative arguments are made positive using + * a reflection formula. + */ + + real q, z; + if (isNaN(x)) return x; + if (x == -x.infinity) return real.nan; + if ( fabs(x) > MAXGAMMA ) return real.infinity; + if (x==0) return 1.0 / x; // +- infinity depending on sign of x, create an exception. + + q = fabs(x); + + if ( q > 13.0L ) { + // Large arguments are handled by Stirling's + // formula. Large negative arguments are made positive using + // the reflection formula. + + if ( x < 0.0L ) { + if (x < -1/real.epsilon) + { + // Large negatives lose all precision + return real.nan; + } + int sgngam = 1; // sign of gamma. + long intpart = cast(long)(q); + if (q == intpart) + return real.nan; // poles for all integers <0. + real p = intpart; + if ( (intpart & 1) == 0 ) + sgngam = -1; + z = q - p; + if ( z > 0.5L ) { + p += 1.0L; + z = q - p; + } + z = q * sin( PI * z ); + z = fabs(z) * gammaStirling(q); + if ( z <= PI/real.max ) return sgngam * real.infinity; + return sgngam * PI/z; + } else { + return gammaStirling(x); + } + } + + // Arguments |x| <= 13 are reduced by recurrence and the function + // approximated by a rational function of degree 7/8 in the + // interval (2,3). + + z = 1.0L; + while ( x >= 3.0L ) { + x -= 1.0L; + z *= x; + } + + while ( x < -0.03125L ) { + z /= x; + x += 1.0L; + } + + if ( x <= 0.03125L ) { + if ( x == 0.0L ) + return real.nan; + else { + if ( x < 0.0L ) { + x = -x; + return z / (x * poly( x, GammaSmallNegCoeffs )); + } else { + return z / (x * poly( x, GammaSmallCoeffs )); + } + } + } + + while ( x < 2.0L ) { + z /= x; + x += 1.0L; + } + if ( x == 2.0L ) return z; + + x -= 2.0L; + return z * poly( x, GammaNumeratorCoeffs ) / poly( x, GammaDenominatorCoeffs ); +} + +unittest { + // gamma(n) = factorial(n-1) if n is an integer. + real fact = 1.0L; + for (int i=1; fact real.mant_dig-15); + fact *= (i*1.0L); + } + assert(gamma(0.0) == real.infinity); + assert(gamma(-0.0) == -real.infinity); + assert(isNaN(gamma(-1.0))); + assert(isNaN(gamma(-15.0))); + assert(isIdentical(gamma(NaN(0xABC)), NaN(0xABC))); + assert(gamma(real.infinity) == real.infinity); + assert(gamma(real.max) == real.infinity); + assert(isNaN(gamma(-real.infinity))); + assert(gamma(real.min*real.epsilon) == real.infinity); + assert(gamma(MAXGAMMA)< real.infinity); + assert(gamma(MAXGAMMA*2) == real.infinity); + + // Test some high-precision values (50 decimal digits) + real SQRT_PI = 1.77245385090551602729816748334114518279754945612238L; + + assert(feqrel(gamma(0.5L), SQRT_PI) == real.mant_dig); + + assert(feqrel(gamma(1.0 / 3.L), 2.67893853470774763365569294097467764412868937795730L) >= real.mant_dig-2); + assert(feqrel(gamma(0.25L), + 3.62560990822190831193068515586767200299516768288006L) >= real.mant_dig-1); + assert(feqrel(gamma(1.0 / 5.0L), + 4.59084371199880305320475827592915200343410999829340L) >= real.mant_dig-1); +} + +/***************************************************** + * Natural logarithm of gamma function. + * + * Returns the base e (2.718...) logarithm of the absolute + * value of the gamma function of the argument. + * + * For reals, logGamma is equivalent to log(fabs(gamma(x))). + * + * $(TABLE_SV + * $(SVH x, logGamma(x) ) + * $(SV $(NAN), $(NAN) ) + * $(SV integer <= 0, +∞ ) + * $(SV ±∞, +∞ ) + * ) + */ +real logGamma(real x) +{ + /* Based on code from the CEPHES library. + * CEPHES code Copyright 1994 by Stephen L. Moshier + * + * For arguments greater than 33, the logarithm of the gamma + * function is approximated by the logarithmic version of + * Stirling's formula using a polynomial approximation of + * degree 4. Arguments between -33 and +33 are reduced by + * recurrence to the interval [2,3] of a rational approximation. + * The cosecant reflection formula is employed for arguments + * less than -33. + */ + real q, w, z, f, nx; + + if (isNaN(x)) return x; + if (fabs(x) == x.infinity) return x.infinity; + + if( x < -34.0L ) + { + q = -x; + w = logGamma(q); + real p = floor(q); + if ( p == q ) + return real.infinity; + int intpart = cast(int)(p); + real sgngam = 1; + if ( (intpart & 1) == 0 ) + sgngam = -1; + z = q - p; + if ( z > 0.5L ) { + p += 1.0L; + z = p - q; + } + z = q * sin( PI * z ); + if ( z == 0.0L ) + return sgngam * real.infinity; + /* z = LOGPI - logl( z ) - w; */ + z = log( PI/z ) - w; + return z; + } + + if( x < 13.0L ) + { + z = 1.0L; + nx = floor( x + 0.5L ); + f = x - nx; + while ( x >= 3.0L ) { + nx -= 1.0L; + x = nx + f; + z *= x; + } + while ( x < 2.0L ) { + if( fabs(x) <= 0.03125 ) + { + if ( x == 0.0L ) + return real.infinity; + if ( x < 0.0L ) + { + x = -x; + q = z / (x * poly( x, GammaSmallNegCoeffs)); + } else + q = z / (x * poly( x, GammaSmallCoeffs)); + return log( fabs(q) ); + } + z /= nx + f; + nx += 1.0L; + x = nx + f; + } + z = fabs(z); + if ( x == 2.0L ) + return log(z); + x = (nx - 2.0L) + f; + real p = x * rationalPoly( x, logGammaNumerator, logGammaDenominator); + return log(z) + p; + } + + // const real MAXLGM = 1.04848146839019521116e+4928L; + // if( x > MAXLGM ) return sgngaml * real.infinity; + + const real LOGSQRT2PI = 0.91893853320467274178L; // log( sqrt( 2*pi ) ) + + q = ( x - 0.5L ) * log(x) - x + LOGSQRT2PI; + if (x > 1.0e10L) return q; + real p = 1.0L / (x*x); + q += poly( p, logGammaStirlingCoeffs ) / x; + return q ; +} + +unittest { + assert(isIdentical(logGamma(NaN(0xDEF)), NaN(0xDEF))); + assert(logGamma(real.infinity) == real.infinity); + assert(logGamma(-1.0) == real.infinity); + assert(logGamma(0.0) == real.infinity); + assert(logGamma(-50.0) == real.infinity); + assert(isIdentical(0.0L, logGamma(1.0L))); + assert(isIdentical(0.0L, logGamma(2.0L))); + assert(logGamma(real.min*real.epsilon) == real.infinity); + assert(logGamma(-real.min*real.epsilon) == real.infinity); + + // x, correct loggamma(x), correct d/dx loggamma(x). + static real[] testpoints = [ + 8.0L, 8.525146484375L + 1.48766904143001655310E-5, 2.01564147795560999654E0L, + 8.99993896484375e-1L, 6.6375732421875e-2L + 5.11505711292524166220E-6L, -7.54938684259372234258E-1, + 7.31597900390625e-1L, 2.2369384765625e-1 + 5.21506341809849792422E-6L, -1.13355566660398608343E0L, + 2.31639862060546875e-1L, 1.3686676025390625L + 1.12609441752996145670E-5L, -4.56670961813812679012E0, + 1.73162841796875L, -8.88214111328125e-2L + 3.36207740803753034508E-6L, 2.33339034686200586920E-1L, + 1.23162841796875L, -9.3902587890625e-2L + 1.28765089229009648104E-5L, -2.49677345775751390414E-1L, + 7.3786976294838206464e19L, 3.301798506038663053312e21L - 1.656137564136932662487046269677E5L, + 4.57477139169563904215E1L, + 1.08420217248550443401E-19L, 4.36682586669921875e1L + 1.37082843669932230418E-5L, + -9.22337203685477580858E18L, + 1.0L, 0.0L, -5.77215664901532860607E-1L, + 2.0L, 0.0L, 4.22784335098467139393E-1L, + -0.5L, 1.2655029296875L + 9.19379714539648894580E-6L, 3.64899739785765205590E-2L, + -1.5L, 8.6004638671875e-1L + 6.28657731014510932682E-7L, 7.03156640645243187226E-1L, + -2.5L, -5.6243896484375E-2L + 1.79986700949327405470E-7, 1.10315664064524318723E0L, + -3.5L, -1.30902099609375L + 1.43111007079536392848E-5L, 1.38887092635952890151E0L + ]; + // TODO: test derivatives as well. + for (int i=0; i real.mant_dig-5); + if (testpoints[i] real.mant_dig-5); + } + } + 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); +} + + +private { +enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max) +enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min*real.epsilon) = log(smallest denormal) +enum real BETA_BIG = 9.223372036854775808e18L; +enum real BETA_BIGINV = 1.084202172485504434007e-19L; +} + +/** Incomplete beta integral + * + * Returns incomplete beta integral of the arguments, evaluated + * from zero to x. The regularized incomplete beta function is defined as + * + * betaIncomplete(a, b, x) = Γ(a+b)/(Γ(a) Γ(b)) * + * $(INTEGRATE 0, x) $(POWER t, a-1)$(POWER (1-t),b-1) dt + * + * and is the same as the the cumulative distribution function. + * + * The domain of definition is 0 <= x <= 1. In this + * implementation a and b are restricted to positive values. + * The integral from x to 1 may be obtained by the symmetry + * relation + * + * betaIncompleteCompl(a, b, x ) = betaIncomplete( b, a, 1-x ) + * + * The integral is evaluated by a continued fraction expansion + * or, when b*x is small, by a power series. + */ +real betaIncomplete(real aa, real bb, real xx ) +{ + if ( !(aa>0 && bb>0) ) + { + if ( isNaN(aa) ) return aa; + if ( isNaN(bb) ) return bb; + return real.nan; // domain error + } + if (!(xx>0 && xx<1.0)) { + if (isNaN(xx)) return xx; + if ( xx == 0.0L ) return 0.0; + if ( xx == 1.0L ) return 1.0; + return real.nan; // domain error + } + if ( (bb * xx) <= 1.0L && xx <= 0.95L) { + return betaDistPowerSeries(aa, bb, xx); + } + real x; + real xc; // = 1 - x + + real a, b; + int flag = 0; + + /* Reverse a and b if x is greater than the mean. */ + if( xx > (aa/(aa+bb)) ) { + // here x > aa/(aa+bb) and (bb*x>1 or x>0.95) + flag = 1; + a = bb; + b = aa; + xc = xx; + x = 1.0L - xx; + } else { + a = aa; + b = bb; + xc = 1.0L - xx; + x = xx; + } + + if( flag == 1 && (b * x) <= 1.0L && x <= 0.95L) { + // here xx > aa/(aa+bb) and ((bb*xx>1) or xx>0.95) and (aa*(1-xx)<=1) and xx > 0.05 + return 1.0 - betaDistPowerSeries(a, b, x); // note loss of precision + } + + real w; + // Choose expansion for optimal convergence + // One is for x * (a+b+2) < (a+1), + // the other is for x * (a+b+2) > (a+1). + real y = x * (a+b-2.0L) - (a-1.0L); + if( y < 0.0L ) { + w = betaDistExpansion1( a, b, x ); + } else { + w = betaDistExpansion2( a, b, x ) / xc; + } + + /* Multiply w by the factor + a b + x (1-x) Gamma(a+b) / ( a Gamma(a) Gamma(b) ) . */ + + y = a * log(x); + real t = b * log(xc); + if ( (a+b) < MAXGAMMA && fabs(y) < MAXLOG && fabs(t) < MAXLOG ) { + t = pow(xc,b); + t *= pow(x,a); + t /= a; + t *= w; + t *= gamma(a+b) / (gamma(a) * gamma(b)); + } else { + /* Resort to logarithms. */ + y += t + logGamma(a+b) - logGamma(a) - logGamma(b); + y += log(w/a); + + t = exp(y); +/+ + // There seems to be a bug in Cephes at this point. + // Problems occur for y > MAXLOG, not y < MINLOG. + if( y < MINLOG ) { + t = 0.0L; + } else { + t = exp(y); + } ++/ + } + if( flag == 1 ) { +/+ // CEPHES includes this code, but I think it is erroneous. + if( t <= real.epsilon ) { + t = 1.0L - real.epsilon; + } else ++/ + t = 1.0L - t; + } + return t; +} + +/** Inverse of incomplete beta integral + * + * Given y, the function finds x such that + * + * betaIncomplete(a, b, x) == y + * + * Newton iterations or interval halving is used. + */ +real betaIncompleteInv(real aa, real bb, real yy0 ) +{ + real a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt; + int i, rflg, dir, nflg; + + if (isNaN(yy0)) return yy0; + if (isNaN(aa)) return aa; + if (isNaN(bb)) return bb; + if( yy0 <= 0.0L ) + return 0.0L; + if( yy0 >= 1.0L ) + return 1.0L; + x0 = 0.0L; + yl = 0.0L; + x1 = 1.0L; + yh = 1.0L; + if( aa <= 1.0L || bb <= 1.0L ) { + dithresh = 1.0e-7L; + rflg = 0; + a = aa; + b = bb; + y0 = yy0; + x = a/(a+b); + y = betaIncomplete( a, b, x ); + nflg = 0; + goto ihalve; + } else { + nflg = 0; + dithresh = 1.0e-4L; + } + + // approximation to inverse function + + yp = -normalDistributionInvImpl( yy0 ); + + if( yy0 > 0.5L ) { + rflg = 1; + a = bb; + b = aa; + y0 = 1.0L - yy0; + yp = -yp; + } else { + rflg = 0; + a = aa; + b = bb; + y0 = yy0; + } + + lgm = (yp * yp - 3.0L)/6.0L; + x = 2.0L/( 1.0L/(2.0L * a-1.0L) + 1.0L/(2.0L * b - 1.0L) ); + d = yp * sqrt( x + lgm ) / x + - ( 1.0L/(2.0L * b - 1.0L) - 1.0L/(2.0L * a - 1.0L) ) + * (lgm + (5.0L/6.0L) - 2.0L/(3.0L * x)); + d = 2.0L * d; + if( d < MINLOG ) { + x = 1.0L; + goto under; + } + x = a/( a + b * exp(d) ); + y = betaIncomplete( a, b, x ); + yp = (y - y0)/y0; + if( fabs(yp) < 0.2 ) + goto newt; + + /* Resort to interval halving if not close enough. */ +ihalve: + + dir = 0; + di = 0.5L; + for( i=0; i<400; i++ ) { + if( i != 0 ) { + x = x0 + di * (x1 - x0); + if( x == 1.0L ) { + x = 1.0L - real.epsilon; + } + if( x == 0.0L ) { + di = 0.5; + x = x0 + di * (x1 - x0); + if( x == 0.0 ) + goto under; + } + y = betaIncomplete( a, b, x ); + yp = (x1 - x0)/(x1 + x0); + if( fabs(yp) < dithresh ) + goto newt; + yp = (y-y0)/y0; + if( fabs(yp) < dithresh ) + goto newt; + } + if( y < y0 ) { + x0 = x; + yl = y; + if( dir < 0 ) { + dir = 0; + di = 0.5L; + } else if( dir > 3 ) + di = 1.0L - (1.0L - di) * (1.0L - di); + else if( dir > 1 ) + di = 0.5L * di + 0.5L; + else + di = (y0 - y)/(yh - yl); + dir += 1; + if( x0 > 0.95L ) { + if( rflg == 1 ) { + rflg = 0; + a = aa; + b = bb; + y0 = yy0; + } else { + rflg = 1; + a = bb; + b = aa; + y0 = 1.0 - yy0; + } + x = 1.0L - x; + y = betaIncomplete( a, b, x ); + x0 = 0.0; + yl = 0.0; + x1 = 1.0; + yh = 1.0; + goto ihalve; + } + } else { + x1 = x; + if( rflg == 1 && x1 < real.epsilon ) { + x = 0.0L; + goto done; + } + yh = y; + if( dir > 0 ) { + dir = 0; + di = 0.5L; + } + else if( dir < -3 ) + di = di * di; + else if( dir < -1 ) + di = 0.5L * di; + else + di = (y - y0)/(yh - yl); + dir -= 1; + } + } + if( x0 >= 1.0L ) { + // partial loss of precision + x = 1.0L - real.epsilon; + goto done; + } + if( x <= 0.0L ) { +under: + // underflow has occurred + x = real.min * real.min; + goto done; + } + +newt: + + if ( nflg ) { + goto done; + } + nflg = 1; + lgm = logGamma(a+b) - logGamma(a) - logGamma(b); + + for( i=0; i<15; i++ ) { + /* Compute the function at this point. */ + if ( i != 0 ) + y = betaIncomplete(a,b,x); + if ( y < yl ) { + x = x0; + y = yl; + } else if( y > yh ) { + x = x1; + y = yh; + } else if( y < y0 ) { + x0 = x; + yl = y; + } else { + x1 = x; + yh = y; + } + if( x == 1.0L || x == 0.0L ) + break; + /* Compute the derivative of the function at this point. */ + d = (a - 1.0L) * log(x) + (b - 1.0L) * log(1.0L - x) + lgm; + if ( d < MINLOG ) { + goto done; + } + if ( d > MAXLOG ) { + break; + } + d = exp(d); + /* Compute the step to the next approximation of x. */ + d = (y - y0)/d; + xt = x - d; + if ( xt <= x0 ) { + y = (x - x0) / (x1 - x0); + xt = x0 + 0.5L * y * (x - x0); + if( xt <= 0.0L ) + break; + } + if ( xt >= x1 ) { + y = (x1 - x) / (x1 - x0); + xt = x1 - 0.5L * y * (x1 - x); + if ( xt >= 1.0L ) + break; + } + x = xt; + if ( fabs(d/x) < (128.0L * real.epsilon) ) + goto done; + } + /* Did not converge. */ + dithresh = 256.0L * real.epsilon; + goto ihalve; + +done: + if ( rflg ) { + if( x <= real.epsilon ) + x = 1.0L - real.epsilon; + else + x = 1.0L - x; + } + return x; +} + +unittest { // also tested by the normal distribution + // check NaN propagation + assert(isIdentical(betaIncomplete(NaN(0xABC),2,3), NaN(0xABC))); + assert(isIdentical(betaIncomplete(7,NaN(0xABC),3), NaN(0xABC))); + assert(isIdentical(betaIncomplete(7,15,NaN(0xABC)), NaN(0xABC))); + assert(isIdentical(betaIncompleteInv(NaN(0xABC),1,17), NaN(0xABC))); + assert(isIdentical(betaIncompleteInv(2,NaN(0xABC),8), NaN(0xABC))); + assert(isIdentical(betaIncompleteInv(2,3, NaN(0xABC)), NaN(0xABC))); + + assert(isNaN(betaIncomplete(-1, 2, 3))); + + assert(betaIncomplete(1, 2, 0)==0); + assert(betaIncomplete(1, 2, 1)==1); + assert(isNaN(betaIncomplete(1, 2, 3))); + assert(betaIncompleteInv(1, 1, 0)==0); + assert(betaIncompleteInv(1, 1, 1)==1); + + // Test some values against Microsoft Excel 2003. + + assert(fabs(betaIncomplete(8, 10, 0.2) - 0.010_934_315_236_957_2L) < 0.000_000_000_5); + assert(fabs(betaIncomplete(2, 2.5, 0.9) - 0.989_722_597_604_107L) < 0.000_000_000_000_5); + assert(fabs(betaIncomplete(1000, 800, 0.5) - 1.17914088832798E-06L) < 0.000_000_05e-6); + + assert(fabs(betaIncomplete(0.0001, 10000, 0.0001) - 0.999978059369989L) < 0.000_000_000_05); + + assert(fabs(betaIncompleteInv(5, 10, 0.2) - 0.229121208190918L) < 0.000_000_5L); + assert(fabs(betaIncompleteInv(4, 7, 0.8) - 0.483657360076904L) < 0.000_000_5L); + + // Coverage tests. I don't have correct values for these tests, but + // these values cover most of the code, so they are useful for + // regression testing. + // 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. + +// Excel 2003 gives clearly erroneous results (betadist>1) when a and x are tiny and b is huge. +// The correct results are for these next tests are unknown. + +// real testpoint1 = betaIncomplete(1e-10, 5e20, 8e-21); +// assert(testpoint1 == 0x1.ffff_ffff_c906_404cp-1L); + + assert(betaIncomplete(0.01, 327726.7, 0.545113) == 1.0); + assert(betaIncompleteInv(0.01, 8e-48, 5.45464e-20)==1-real.epsilon); + assert(betaIncompleteInv(0.01, 8e-48, 9e-26)==1-real.epsilon); + + assert(betaIncomplete(0.01, 498.437, 0.0121433) == 0x1.ffff_8f72_19197402p-1); + assert(1- betaIncomplete(0.01, 328222, 4.0375e-5) == 0x1.5f62926b4p-30); + + // 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); + assert(betaIncompleteInv(0x1.ff1275ae5b939bcap-41, 4.6713e18, 0.0813601)==0x1.f97749d90c7adba8p-63); + 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. + assert(betaIncomplete(1.16251e20, 2.18e39, 5.45e-20)==-real.infinity); +} + + +private { +// Implementation functions + +// Continued fraction expansion #1 for incomplete beta integral +// Use when x < (a+1)/(a+b+2) +real betaDistExpansion1(real a, real b, real x ) +{ + real xk, pk, pkm1, pkm2, qk, qkm1, qkm2; + real k1, k2, k3, k4, k5, k6, k7, k8; + real r, t, ans; + int n; + + k1 = a; + k2 = a + b; + k3 = a; + k4 = a + 1.0L; + k5 = 1.0L; + k6 = b - 1.0L; + k7 = k4; + k8 = a + 2.0L; + + pkm2 = 0.0L; + qkm2 = 1.0L; + pkm1 = 1.0L; + qkm1 = 1.0L; + ans = 1.0L; + r = 1.0L; + n = 0; + const real thresh = 3.0L * real.epsilon; + do { + xk = -( x * k1 * k2 )/( k3 * k4 ); + pk = pkm1 + pkm2 * xk; + qk = qkm1 + qkm2 * xk; + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + + xk = ( x * k5 * k6 )/( k7 * k8 ); + pk = pkm1 + pkm2 * xk; + qk = qkm1 + qkm2 * xk; + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + + if( qk != 0.0L ) + r = pk/qk; + if( r != 0.0L ) { + t = fabs( (ans - r)/r ); + ans = r; + } else { + t = 1.0L; + } + + if( t < thresh ) + return ans; + + k1 += 1.0L; + k2 += 1.0L; + k3 += 2.0L; + k4 += 2.0L; + k5 += 1.0L; + k6 -= 1.0L; + k7 += 2.0L; + k8 += 2.0L; + + if( (fabs(qk) + fabs(pk)) > BETA_BIG ) { + pkm2 *= BETA_BIGINV; + pkm1 *= BETA_BIGINV; + qkm2 *= BETA_BIGINV; + qkm1 *= BETA_BIGINV; + } + if( (fabs(qk) < BETA_BIGINV) || (fabs(pk) < BETA_BIGINV) ) { + pkm2 *= BETA_BIG; + pkm1 *= BETA_BIG; + qkm2 *= BETA_BIG; + qkm1 *= BETA_BIG; + } + } + while( ++n < 400 ); +// loss of precision has occurred +// mtherr( "incbetl", PLOSS ); + return ans; +} + +// Continued fraction expansion #2 for incomplete beta integral +// Use when x > (a+1)/(a+b+2) +real betaDistExpansion2(real a, real b, real x ) +{ + real xk, pk, pkm1, pkm2, qk, qkm1, qkm2; + real k1, k2, k3, k4, k5, k6, k7, k8; + real r, t, ans, z; + + k1 = a; + k2 = b - 1.0L; + k3 = a; + k4 = a + 1.0L; + k5 = 1.0L; + k6 = a + b; + k7 = a + 1.0L; + k8 = a + 2.0L; + + pkm2 = 0.0L; + qkm2 = 1.0L; + pkm1 = 1.0L; + qkm1 = 1.0L; + z = x / (1.0L-x); + ans = 1.0L; + r = 1.0L; + int n = 0; + const real thresh = 3.0L * real.epsilon; + do { + + xk = -( z * k1 * k2 )/( k3 * k4 ); + pk = pkm1 + pkm2 * xk; + qk = qkm1 + qkm2 * xk; + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + + xk = ( z * k5 * k6 )/( k7 * k8 ); + pk = pkm1 + pkm2 * xk; + qk = qkm1 + qkm2 * xk; + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + + if( qk != 0.0L ) + r = pk/qk; + if( r != 0.0L ) { + t = fabs( (ans - r)/r ); + ans = r; + } else + t = 1.0L; + + if( t < thresh ) + return ans; + k1 += 1.0L; + k2 -= 1.0L; + k3 += 2.0L; + k4 += 2.0L; + k5 += 1.0L; + k6 += 1.0L; + k7 += 2.0L; + k8 += 2.0L; + + if( (fabs(qk) + fabs(pk)) > BETA_BIG ) { + pkm2 *= BETA_BIGINV; + pkm1 *= BETA_BIGINV; + qkm2 *= BETA_BIGINV; + qkm1 *= BETA_BIGINV; + } + if( (fabs(qk) < BETA_BIGINV) || (fabs(pk) < BETA_BIGINV) ) { + pkm2 *= BETA_BIG; + pkm1 *= BETA_BIG; + qkm2 *= BETA_BIG; + qkm1 *= BETA_BIG; + } + } while( ++n < 400 ); +// loss of precision has occurred +//mtherr( "incbetl", PLOSS ); + return ans; +} + +/* Power series for incomplete gamma integral. + Use when b*x is small. */ +real betaDistPowerSeries(real a, real b, real x ) +{ + real ai = 1.0L / a; + real u = (1.0L - b) * x; + real v = u / (a + 1.0L); + real t1 = v; + real t = u; + real n = 2.0L; + real s = 0.0L; + real z = real.epsilon * ai; + while( fabs(v) > z ) { + u = (n - b) * x / n; + t *= u; + v = t / (a + n); + s += v; + n += 1.0L; + } + s += t1; + s += ai; + + u = a * log(x); + if ( (a+b) < MAXGAMMA && fabs(u) < MAXLOG ) { + t = gamma(a+b)/(gamma(a)*gamma(b)); + s = s * t * pow(x,a); + } else { + t = logGamma(a+b) - logGamma(a) - logGamma(b) + u + log(s); + + if( t < MINLOG ) { + s = 0.0L; + } else + s = exp(t); + } + return s; +} + +} + +/*************************************** + * Incomplete gamma integral and its complement + * + * These functions are defined by + * + * gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a) + * + * gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x) + * = ($(INTEGRATE x, ∞) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a) + * + * In this implementation both arguments must be positive. + * The integral is evaluated by either a power series or + * continued fraction expansion, depending on the relative + * values of a and x. + */ +real gammaIncomplete(real a, real x ) +in { + assert(x >= 0); + assert(a > 0); +} +body { + /* left tail of incomplete gamma function: + * + * inf. k + * a -x - x + * x e > ---------- + * - - + * k=0 | (a+k+1) + * + */ + if (x==0) + return 0.0L; + + if ( (x > 1.0L) && (x > a ) ) + return 1.0L - gammaIncompleteCompl(a,x); + + real ax = a * log(x) - x - logGamma(a); +/+ + if( ax < MINLOGL ) return 0; // underflow + // { mtherr( "igaml", UNDERFLOW ); return( 0.0L ); } ++/ + ax = exp(ax); + + /* power series */ + real r = a; + real c = 1.0L; + real ans = 1.0L; + + do { + r += 1.0L; + c *= x/r; + ans += c; + } while( c/ans > real.epsilon ); + + return ans * ax/a; +} + +/** ditto */ +real gammaIncompleteCompl(real a, real x ) +in { + assert(x >= 0); + assert(a > 0); +} +body { + if (x==0) + return 1.0L; + if ( (x < 1.0L) || (x < a) ) + return 1.0L - gammaIncomplete(a,x); + + // DAC (Cephes bug fix): This is necessary to avoid + // spurious nans, eg + // log(x)-x = NaN when x = real.infinity + const real MAXLOGL = 1.1356523406294143949492E4L; + if (x > MAXLOGL) return 0; // underflow + + real ax = a * log(x) - x - logGamma(a); +//const real MINLOGL = -1.1355137111933024058873E4L; +// if ( ax < MINLOGL ) return 0; // underflow; + ax = exp(ax); + + + /* continued fraction */ + real y = 1.0L - a; + real z = x + y + 1.0L; + real c = 0.0L; + + real pk, qk, t; + + real pkm2 = 1.0L; + real qkm2 = x; + real pkm1 = x + 1.0L; + real qkm1 = z * x; + real ans = pkm1/qkm1; + + do { + c += 1.0L; + y += 1.0L; + z += 2.0L; + real yc = y * c; + pk = pkm1 * z - pkm2 * yc; + qk = qkm1 * z - qkm2 * yc; + if( qk != 0.0L ) { + real r = pk/qk; + t = fabs( (ans - r)/r ); + ans = r; + } else { + t = 1.0L; + } + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + + const real BIG = 9.223372036854775808e18L; + + if ( fabs(pk) > BIG ) { + pkm2 /= BIG; + pkm1 /= BIG; + qkm2 /= BIG; + qkm1 /= BIG; + } + } while ( t > real.epsilon ); + + return ans * ax; +} + +/** Inverse of complemented incomplete gamma integral + * + * Given a and y, the function finds x such that + * + * gammaIncompleteCompl( a, x ) = p. + * + * Starting with the approximate value x = a $(POWER t, 3), where + * t = 1 - d - normalDistributionInv(p) sqrt(d), + * and d = 1/9a, + * the routine performs up to 10 Newton iterations to find the + * root of incompleteGammaCompl(a,x) - p = 0. + */ +real gammaIncompleteComplInv(real a, real p) +in { + assert(p>=0 && p<= 1); + assert(a>0); +} +body { + if (p==0) return real.infinity; + + real y0 = p; + const real MAXLOGL = 1.1356523406294143949492E4L; + real x0, x1, x, yl, yh, y, d, lgm, dithresh; + int i, dir; + + /* bound the solution */ + x0 = real.max; + yl = 0.0L; + x1 = 0.0L; + yh = 1.0L; + dithresh = 4.0 * real.epsilon; + + /* approximation to inverse function */ + d = 1.0L/(9.0L*a); + y = 1.0L - d - normalDistributionInvImpl(y0) * sqrt(d); + x = a * y * y * y; + + lgm = logGamma(a); + + for( i=0; i<10; i++ ) { + if( x > x0 || x < x1 ) + goto ihalve; + y = gammaIncompleteCompl(a,x); + if ( y < yl || y > yh ) + goto ihalve; + if ( y < y0 ) { + x0 = x; + yl = y; + } else { + x1 = x; + yh = y; + } + /* compute the derivative of the function at this point */ + d = (a - 1.0L) * log(x0) - x0 - lgm; + if ( d < -MAXLOGL ) + goto ihalve; + d = -exp(d); + /* compute the step to the next approximation of x */ + d = (y - y0)/d; + x = x - d; + if ( i < 3 ) continue; + if ( fabs(d/x) < dithresh ) return x; + } + + /* Resort to interval halving if Newton iteration did not converge. */ +ihalve: + d = 0.0625L; + if ( x0 == real.max ) { + if( x <= 0.0L ) + x = 1.0L; + while( x0 == real.max ) { + x = (1.0L + d) * x; + y = gammaIncompleteCompl( a, x ); + if ( y < y0 ) { + x0 = x; + yl = y; + break; + } + d = d + d; + } + } + d = 0.5L; + dir = 0; + + for( i=0; i<400; i++ ) { + x = x1 + d * (x0 - x1); + y = gammaIncompleteCompl( a, x ); + lgm = (x0 - x1)/(x1 + x0); + if ( fabs(lgm) < dithresh ) + break; + lgm = (y - y0)/y0; + if ( fabs(lgm) < dithresh ) + break; + if ( x <= 0.0L ) + break; + if ( y > y0 ) { + x1 = x; + yh = y; + if ( dir < 0 ) { + dir = 0; + d = 0.5L; + } else if ( dir > 1 ) + d = 0.5L * d + 0.5L; + else + d = (y0 - yl)/(yh - yl); + dir += 1; + } else { + x0 = x; + yl = y; + if ( dir > 0 ) { + dir = 0; + d = 0.5L; + } else if ( dir < -1 ) + d = 0.5L * d; + else + d = (y0 - yl)/(yh - yl); + dir -= 1; + } + } + /+ + if( x == 0.0L ) + mtherr( "igamil", UNDERFLOW ); + +/ + return x; +} + +unittest { +//Values from Excel's GammaInv(1-p, x, 1) +assert(fabs(gammaIncompleteComplInv(1, 0.5) - 0.693147188044814) < 0.00000005); +assert(fabs(gammaIncompleteComplInv(12, 0.99) - 5.42818075054289) < 0.00000005); +assert(fabs(gammaIncompleteComplInv(100, 0.8) - 91.5013985848288L) < 0.000005); + +assert(gammaIncomplete(1, 0)==0); +assert(gammaIncompleteCompl(1, 0)==1); +assert(gammaIncomplete(4545, real.infinity)==1); + +// Values from Excel's (1-GammaDist(x, alpha, 1, TRUE)) + +assert(fabs(1.0L-gammaIncompleteCompl(0.5, 2) - 0.954499729507309L) < 0.00000005); +assert(fabs(gammaIncomplete(0.5, 2) - 0.954499729507309L) < 0.00000005); +// Fixed Cephes bug: +assert(gammaIncompleteCompl(384, real.infinity)==0); +assert(gammaIncompleteComplInv(3, 0)==real.infinity); +} + + +/** Digamma function +* +* The digamma function is the logarithmic derivative of the gamma function. +* +* digamma(x) = d/dx logGamma(x) +* +*/ +real digamma(real x) +{ + // Based on CEPHES, Stephen L. Moshier. + + // DAC: These values are Bn / n for n=2,4,6,8,10,12,14. + const real [] Bn_n = [ + 1.0L/(6*2), -1.0L/(30*4), 1.0L/(42*6), -1.0L/(30*8), + 5.0L/(66*10), -691.0L/(2730*12), 7.0L/(6*14) ]; + + real p, q, nz, s, w, y, z; + long i, n; + int negative; + + negative = 0; + nz = 0.0; + + if ( x <= 0.0 ) { + negative = 1; + q = x; + p = floor(q); + if( p == q ) { + return real.nan; // singularity. + } + /* Remove the zeros of tan(PI x) + * by subtracting the nearest integer from x + */ + nz = q - p; + if ( nz != 0.5 ) { + if ( nz > 0.5 ) { + p += 1.0; + nz = q - p; + } + nz = PI/tan(PI*nz); + } else { + nz = 0.0; + } + x = 1.0 - x; + } + + // check for small positive integer + if ((x <= 13.0) && (x == floor(x)) ) { + y = 0.0; + n = lrint(x); + // DAC: CEPHES bugfix. Cephes did this in reverse order, which + // created a larger roundoff error. + for (i=n-1; i>0; --i) { + y+=1.0L/i; + } + y -= EULERGAMMA; + goto done; + } + + s = x; + w = 0.0; + while ( s < 10.0 ) { + w += 1.0/s; + s += 1.0; + } + + if ( s < 1.0e17 ) { + z = 1.0/(s * s); + y = z * poly(z, Bn_n); + } else + y = 0.0; + + y = log(s) - 0.5L/s - y - w; + +done: + if ( negative ) { + y -= nz; + } + return y; +} + +version(unittest) import std.c.stdio; +unittest { + // Exact values + assert(digamma(1)== -EULERGAMMA); + assert(feqrel(digamma(0.25), -PI/2 - 3* LN2 - EULERGAMMA)>=real.mant_dig-7); + assert(feqrel(digamma(1.0L/6), -PI/2 *sqrt(3.0L) - 2* LN2 -1.5*log(3.0L) - EULERGAMMA)>=real.mant_dig-7); + assert(digamma(-5)!<>0); + assert(feqrel(digamma(2.5), -EULERGAMMA - 2*LN2 + 2.0 + 2.0L/3)>=real.mant_dig-9); + assert(isIdentical(digamma(NaN(0xABC)), NaN(0xABC))); + + for (int k=1; k<40; ++k) { + real y=0; + for (int u=k; u>=1; --u) { + y+= 1.0L/u; + } + assert(feqrel(digamma(k+1),-EULERGAMMA + y) >=real.mant_dig-2); + } +}