/std/internal/math/gammafunction.d
D | 1518 lines | 1096 code | 144 blank | 278 comment | 252 complexity | a20985dd8c78a99da3c02d40dcce0143 MD5 | raw file
1/** 2 * Implementation of the gamma and beta functions, and their integrals. 3 * 4 * License: $(WEB boost.org/LICENSE_1_0.txt, Boost License 1.0). 5 * Copyright: Based on the CEPHES math library, which is 6 * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com). 7 * Authors: Stephen L. Moshier (original C code). Conversion to D by Don Clugston 8 * 9 * 10Macros: 11 * TABLE_SV = <table border=1 cellpadding=4 cellspacing=0> 12 * <caption>Special Values</caption> 13 * $0</table> 14 * SVH = $(TR $(TH $1) $(TH $2)) 15 * SV = $(TR $(TD $1) $(TD $2)) 16 * GAMMA = Γ 17 * INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>) 18 * POWER = $1<sup>$2</sup> 19 * NAN = $(RED NAN) 20 */ 21module std.internal.math.gammafunction; 22import std.internal.math.errorfunction; 23import std.math; 24 25private { 26 27enum real SQRT2PI = 2.50662827463100050242E0L; // sqrt(2pi) 28immutable real EULERGAMMA = 0.57721_56649_01532_86060_65120_90082_40243_10421_59335_93992L; /** Euler-Mascheroni constant 0.57721566.. */ 29 30// Polynomial approximations for gamma and loggamma. 31 32immutable real[8] GammaNumeratorCoeffs = [ 1.0, 33 0x1.acf42d903366539ep-1, 0x1.73a991c8475f1aeap-2, 0x1.c7e918751d6b2a92p-4, 34 0x1.86d162cca32cfe86p-6, 0x1.0c378e2e6eaf7cd8p-8, 0x1.dc5c66b7d05feb54p-12, 35 0x1.616457b47e448694p-15 36]; 37 38immutable real[9] GammaDenominatorCoeffs = [ 1.0, 39 0x1.a8f9faae5d8fc8bp-2, -0x1.cb7895a6756eebdep-3, -0x1.7b9bab006d30652ap-5, 40 0x1.c671af78f312082ep-6, -0x1.a11ebbfaf96252dcp-11, -0x1.447b4d2230a77ddap-10, 41 0x1.ec1d45bb85e06696p-13,-0x1.d4ce24d05bd0a8e6p-17 42]; 43 44immutable real[9] GammaSmallCoeffs = [ 1.0, 45 0x1.2788cfc6fb618f52p-1, -0x1.4fcf4026afa2f7ecp-1, -0x1.5815e8fa24d7e306p-5, 46 0x1.5512320aea2ad71ap-3, -0x1.59af0fb9d82e216p-5, -0x1.3b4b61d3bfdf244ap-7, 47 0x1.d9358e9d9d69fd34p-8, -0x1.38fc4bcbada775d6p-10 48]; 49 50immutable real[9] GammaSmallNegCoeffs = [ -1.0, 51 0x1.2788cfc6fb618f54p-1, 0x1.4fcf4026afa2bc4cp-1, -0x1.5815e8fa2468fec8p-5, 52 -0x1.5512320baedaf4b6p-3, -0x1.59af0fa283baf07ep-5, 0x1.3b4a70de31e05942p-7, 53 0x1.d9398be3bad13136p-8, 0x1.291b73ee05bcbba2p-10 54]; 55 56immutable real[7] logGammaStirlingCoeffs = [ 57 0x1.5555555555553f98p-4, -0x1.6c16c16c07509b1p-9, 0x1.a01a012461cbf1e4p-11, 58 -0x1.3813089d3f9d164p-11, 0x1.b911a92555a277b8p-11, -0x1.ed0a7b4206087b22p-10, 59 0x1.402523859811b308p-8 60]; 61 62immutable real[7] logGammaNumerator = [ 63 -0x1.0edd25913aaa40a2p+23, -0x1.31c6ce2e58842d1ep+24, -0x1.f015814039477c3p+23, 64 -0x1.74ffe40c4b184b34p+22, -0x1.0d9c6d08f9eab55p+20, -0x1.54c6b71935f1fc88p+16, 65 -0x1.0e761b42932b2aaep+11 66]; 67 68immutable real[8] logGammaDenominator = [ 69 -0x1.4055572d75d08c56p+24, -0x1.deeb6013998e4d76p+24, -0x1.106f7cded5dcc79ep+24, 70 -0x1.25e17184848c66d2p+22, -0x1.301303b99a614a0ap+19, -0x1.09e76ab41ae965p+15, 71 -0x1.00f95ced9e5f54eep+9, 1.0 72]; 73 74/* 75 * Helper function: Gamma function computed by Stirling's formula. 76 * 77 * Stirling's formula for the gamma function is: 78 * 79 * $(GAMMA)(x) = sqrt(2 π) x<sup>x-0.5</sup> exp(-x) (1 + 1/x P(1/x)) 80 * 81 */ 82real gammaStirling(real x) 83{ 84 // CEPHES code Copyright 1994 by Stephen L. Moshier 85 86 static immutable real[9] SmallStirlingCoeffs = [ 87 0x1.55555555555543aap-4, 0x1.c71c71c720dd8792p-9, -0x1.5f7268f0b5907438p-9, 88 -0x1.e13cd410e0477de6p-13, 0x1.9b0f31643442616ep-11, 0x1.2527623a3472ae08p-14, 89 -0x1.37f6bc8ef8b374dep-11,-0x1.8c968886052b872ap-16, 0x1.76baa9c6d3eeddbcp-11 90 ]; 91 92 static immutable real[7] LargeStirlingCoeffs = [ 1.0L, 93 8.33333333333333333333E-2L, 3.47222222222222222222E-3L, 94 -2.68132716049382716049E-3L, -2.29472093621399176955E-4L, 95 7.84039221720066627474E-4L, 6.97281375836585777429E-5L 96 ]; 97 98 real w = 1.0L/x; 99 real y = exp(x); 100 if ( x > 1024.0L ) { 101 // For large x, use rational coefficients from the analytical expansion. 102 w = poly(w, LargeStirlingCoeffs); 103 // Avoid overflow in pow() 104 real v = pow( x, 0.5L * x - 0.25L ); 105 y = v * (v / y); 106 } 107 else { 108 w = 1.0L + w * poly( w, SmallStirlingCoeffs); 109 y = pow( x, x - 0.5L ) / y; 110 } 111 y = SQRT2PI * y * w; 112 return y; 113} 114 115/* 116 * Helper function: Incomplete gamma function computed by Temme's expansion. 117 * 118 * This is a port of igamma_temme_large from Boost. 119 * 120 */ 121real igammaTemmeLarge(real a, real x) 122{ 123 static immutable real[][13] coef = [ 124 [ -0.333333333333333333333, 0.0833333333333333333333, 125 -0.0148148148148148148148, 0.00115740740740740740741, 126 0.000352733686067019400353, -0.0001787551440329218107, 127 0.39192631785224377817e-4, -0.218544851067999216147e-5, 128 -0.18540622107151599607e-5, 0.829671134095308600502e-6, 129 -0.176659527368260793044e-6, 0.670785354340149858037e-8, 130 0.102618097842403080426e-7, -0.438203601845335318655e-8, 131 0.914769958223679023418e-9, -0.255141939949462497669e-10, 132 -0.583077213255042506746e-10, 0.243619480206674162437e-10, 133 -0.502766928011417558909e-11 ], 134 [ -0.00185185185185185185185, -0.00347222222222222222222, 135 0.00264550264550264550265, -0.000990226337448559670782, 136 0.000205761316872427983539, -0.40187757201646090535e-6, 137 -0.18098550334489977837e-4, 0.764916091608111008464e-5, 138 -0.161209008945634460038e-5, 0.464712780280743434226e-8, 139 0.137863344691572095931e-6, -0.575254560351770496402e-7, 140 0.119516285997781473243e-7, -0.175432417197476476238e-10, 141 -0.100915437106004126275e-8, 0.416279299184258263623e-9, 142 -0.856390702649298063807e-10 ], 143 [ 0.00413359788359788359788, -0.00268132716049382716049, 144 0.000771604938271604938272, 0.200938786008230452675e-5, 145 -0.000107366532263651605215, 0.529234488291201254164e-4, 146 -0.127606351886187277134e-4, 0.342357873409613807419e-7, 147 0.137219573090629332056e-5, -0.629899213838005502291e-6, 148 0.142806142060642417916e-6, -0.204770984219908660149e-9, 149 -0.140925299108675210533e-7, 0.622897408492202203356e-8, 150 -0.136704883966171134993e-8 ], 151 [ 0.000649434156378600823045, 0.000229472093621399176955, 152 -0.000469189494395255712128, 0.000267720632062838852962, 153 -0.756180167188397641073e-4, -0.239650511386729665193e-6, 154 0.110826541153473023615e-4, -0.56749528269915965675e-5, 155 0.142309007324358839146e-5, -0.278610802915281422406e-10, 156 -0.169584040919302772899e-6, 0.809946490538808236335e-7, 157 -0.191111684859736540607e-7 ], 158 [ -0.000861888290916711698605, 0.000784039221720066627474, 159 -0.000299072480303190179733, -0.146384525788434181781e-5, 160 0.664149821546512218666e-4, -0.396836504717943466443e-4, 161 0.113757269706784190981e-4, 0.250749722623753280165e-9, 162 -0.169541495365583060147e-5, 0.890750753220530968883e-6, 163 -0.229293483400080487057e-6], 164 [ -0.000336798553366358150309, -0.697281375836585777429e-4, 165 0.000277275324495939207873, -0.000199325705161888477003, 166 0.679778047793720783882e-4, 0.141906292064396701483e-6, 167 -0.135940481897686932785e-4, 0.801847025633420153972e-5, 168 -0.229148117650809517038e-5 ], 169 [ 0.000531307936463992223166, -0.000592166437353693882865, 170 0.000270878209671804482771, 0.790235323266032787212e-6, 171 -0.815396936756196875093e-4, 0.561168275310624965004e-4, 172 -0.183291165828433755673e-4, -0.307961345060330478256e-8, 173 0.346515536880360908674e-5, -0.20291327396058603727e-5, 174 0.57887928631490037089e-6 ], 175 [ 0.000344367606892377671254, 0.517179090826059219337e-4, 176 -0.000334931610811422363117, 0.000281269515476323702274, 177 -0.000109765822446847310235, -0.127410090954844853795e-6, 178 0.277444515115636441571e-4, -0.182634888057113326614e-4, 179 0.578769494973505239894e-5 ], 180 [ -0.000652623918595309418922, 0.000839498720672087279993, 181 -0.000438297098541721005061, -0.696909145842055197137e-6, 182 0.000166448466420675478374, -0.000127835176797692185853, 183 0.462995326369130429061e-4 ], 184 [ -0.000596761290192746250124, -0.720489541602001055909e-4, 185 0.000678230883766732836162, -0.0006401475260262758451, 186 0.000277501076343287044992 ], 187 [ 0.00133244544948006563713, -0.0019144384985654775265, 188 0.00110893691345966373396 ], 189 [ 0.00157972766073083495909, 0.000162516262783915816899, 190 -0.00206334210355432762645, 0.00213896861856890981541, 191 -0.00101085593912630031708 ], 192 [ -0.00407251211951401664727, 0.00640336283380806979482, 193 -0.00404101610816766177474 ] 194 ]; 195 196 // avoid nans when one of the arguments is inf: 197 if(x == real.infinity && a != real.infinity) 198 return 0; 199 200 if(x != real.infinity && a == real.infinity) 201 return 1; 202 203 real sigma = (x - a) / a; 204 real phi = sigma - log(sigma + 1); 205 206 real y = a * phi; 207 real z = sqrt(2 * phi); 208 if(x < a) 209 z = -z; 210 211 real[13] workspace; 212 foreach(i; 0 .. coef.length) 213 workspace[i] = poly(z, coef[i]); 214 215 real result = poly(1 / a, workspace); 216 result *= exp(-y) / sqrt(2 * PI * a); 217 if(x < a) 218 result = -result; 219 220 result += erfc(sqrt(y)) / 2; 221 222 return result; 223} 224 225} // private 226 227public: 228/// The maximum value of x for which gamma(x) < real.infinity. 229enum real MAXGAMMA = 1755.5483429L; 230 231 232/***************************************************** 233 * The Gamma function, $(GAMMA)(x) 234 * 235 * $(GAMMA)(x) is a generalisation of the factorial function 236 * to real and complex numbers. 237 * Like x!, $(GAMMA)(x+1) = x*$(GAMMA)(x). 238 * 239 * Mathematically, if z.re > 0 then 240 * $(GAMMA)(z) = $(INTEGRATE 0, ∞) $(POWER t, z-1)$(POWER e, -t) dt 241 * 242 * $(TABLE_SV 243 * $(SVH x, $(GAMMA)(x) ) 244 * $(SV $(NAN), $(NAN) ) 245 * $(SV ±0.0, ±∞) 246 * $(SV integer > 0, (x-1)! ) 247 * $(SV integer < 0, $(NAN) ) 248 * $(SV +∞, +∞ ) 249 * $(SV -∞, $(NAN) ) 250 * ) 251 */ 252real gamma(real x) 253{ 254/* Based on code from the CEPHES library. 255 * CEPHES code Copyright 1994 by Stephen L. Moshier 256 * 257 * Arguments |x| <= 13 are reduced by recurrence and the function 258 * approximated by a rational function of degree 7/8 in the 259 * interval (2,3). Large arguments are handled by Stirling's 260 * formula. Large negative arguments are made positive using 261 * a reflection formula. 262 */ 263 264 real q, z; 265 if (isNaN(x)) return x; 266 if (x == -x.infinity) return real.nan; 267 if ( fabs(x) > MAXGAMMA ) return real.infinity; 268 if (x==0) return 1.0 / x; // +- infinity depending on sign of x, create an exception. 269 270 q = fabs(x); 271 272 if ( q > 13.0L ) { 273 // Large arguments are handled by Stirling's 274 // formula. Large negative arguments are made positive using 275 // the reflection formula. 276 277 if ( x < 0.0L ) { 278 if (x < -1/real.epsilon) 279 { 280 // Large negatives lose all precision 281 return real.nan; 282 } 283 int sgngam = 1; // sign of gamma. 284 long intpart = cast(long)(q); 285 if (q == intpart) 286 return real.nan; // poles for all integers <0. 287 real p = intpart; 288 if ( (intpart & 1) == 0 ) 289 sgngam = -1; 290 z = q - p; 291 if ( z > 0.5L ) { 292 p += 1.0L; 293 z = q - p; 294 } 295 z = q * sin( PI * z ); 296 z = fabs(z) * gammaStirling(q); 297 if ( z <= PI/real.max ) return sgngam * real.infinity; 298 return sgngam * PI/z; 299 } else { 300 return gammaStirling(x); 301 } 302 } 303 304 // Arguments |x| <= 13 are reduced by recurrence and the function 305 // approximated by a rational function of degree 7/8 in the 306 // interval (2,3). 307 308 z = 1.0L; 309 while ( x >= 3.0L ) { 310 x -= 1.0L; 311 z *= x; 312 } 313 314 while ( x < -0.03125L ) { 315 z /= x; 316 x += 1.0L; 317 } 318 319 if ( x <= 0.03125L ) { 320 if ( x == 0.0L ) 321 return real.nan; 322 else { 323 if ( x < 0.0L ) { 324 x = -x; 325 return z / (x * poly( x, GammaSmallNegCoeffs )); 326 } else { 327 return z / (x * poly( x, GammaSmallCoeffs )); 328 } 329 } 330 } 331 332 while ( x < 2.0L ) { 333 z /= x; 334 x += 1.0L; 335 } 336 if ( x == 2.0L ) return z; 337 338 x -= 2.0L; 339 return z * poly( x, GammaNumeratorCoeffs ) / poly( x, GammaDenominatorCoeffs ); 340} 341 342unittest { 343 // gamma(n) = factorial(n-1) if n is an integer. 344 real fact = 1.0L; 345 for (int i=1; fact<real.max; ++i) { 346 // Require exact equality for small factorials 347 if (i<14) assert(gamma(i*1.0L) == fact); 348 assert(feqrel(gamma(i*1.0L), fact) >= real.mant_dig-15); 349 fact *= (i*1.0L); 350 } 351 assert(gamma(0.0) == real.infinity); 352 assert(gamma(-0.0) == -real.infinity); 353 assert(isNaN(gamma(-1.0))); 354 assert(isNaN(gamma(-15.0))); 355 assert(isIdentical(gamma(NaN(0xABC)), NaN(0xABC))); 356 assert(gamma(real.infinity) == real.infinity); 357 assert(gamma(real.max) == real.infinity); 358 assert(isNaN(gamma(-real.infinity))); 359 assert(gamma(real.min_normal*real.epsilon) == real.infinity); 360 assert(gamma(MAXGAMMA)< real.infinity); 361 assert(gamma(MAXGAMMA*2) == real.infinity); 362 363 // Test some high-precision values (50 decimal digits) 364 real SQRT_PI = 1.77245385090551602729816748334114518279754945612238L; 365 366 367 assert(feqrel(gamma(0.5L), SQRT_PI) >= real.mant_dig-1); 368 assert(feqrel(gamma(17.25L), 4.224986665692703551570937158682064589938e13L) >= real.mant_dig-4); 369 370 assert(feqrel(gamma(1.0 / 3.0L), 2.67893853470774763365569294097467764412868937795730L) >= real.mant_dig-2); 371 assert(feqrel(gamma(0.25L), 372 3.62560990822190831193068515586767200299516768288006L) >= real.mant_dig-1); 373 assert(feqrel(gamma(1.0 / 5.0L), 374 4.59084371199880305320475827592915200343410999829340L) >= real.mant_dig-1); 375} 376 377/***************************************************** 378 * Natural logarithm of gamma function. 379 * 380 * Returns the base e (2.718...) logarithm of the absolute 381 * value of the gamma function of the argument. 382 * 383 * For reals, logGamma is equivalent to log(fabs(gamma(x))). 384 * 385 * $(TABLE_SV 386 * $(SVH x, logGamma(x) ) 387 * $(SV $(NAN), $(NAN) ) 388 * $(SV integer <= 0, +∞ ) 389 * $(SV ±∞, +∞ ) 390 * ) 391 */ 392real logGamma(real x) 393{ 394 /* Based on code from the CEPHES library. 395 * CEPHES code Copyright 1994 by Stephen L. Moshier 396 * 397 * For arguments greater than 33, the logarithm of the gamma 398 * function is approximated by the logarithmic version of 399 * Stirling's formula using a polynomial approximation of 400 * degree 4. Arguments between -33 and +33 are reduced by 401 * recurrence to the interval [2,3] of a rational approximation. 402 * The cosecant reflection formula is employed for arguments 403 * less than -33. 404 */ 405 real q, w, z, f, nx; 406 407 if (isNaN(x)) return x; 408 if (fabs(x) == x.infinity) return x.infinity; 409 410 if( x < -34.0L ) 411 { 412 q = -x; 413 w = logGamma(q); 414 real p = floor(q); 415 if ( p == q ) 416 return real.infinity; 417 int intpart = cast(int)(p); 418 real sgngam = 1; 419 if ( (intpart & 1) == 0 ) 420 sgngam = -1; 421 z = q - p; 422 if ( z > 0.5L ) { 423 p += 1.0L; 424 z = p - q; 425 } 426 z = q * sin( PI * z ); 427 if ( z == 0.0L ) 428 return sgngam * real.infinity; 429 /* z = LOGPI - logl( z ) - w; */ 430 z = log( PI/z ) - w; 431 return z; 432 } 433 434 if( x < 13.0L ) 435 { 436 z = 1.0L; 437 nx = floor( x + 0.5L ); 438 f = x - nx; 439 while ( x >= 3.0L ) { 440 nx -= 1.0L; 441 x = nx + f; 442 z *= x; 443 } 444 while ( x < 2.0L ) { 445 if( fabs(x) <= 0.03125 ) 446 { 447 if ( x == 0.0L ) 448 return real.infinity; 449 if ( x < 0.0L ) 450 { 451 x = -x; 452 q = z / (x * poly( x, GammaSmallNegCoeffs)); 453 } else 454 q = z / (x * poly( x, GammaSmallCoeffs)); 455 return log( fabs(q) ); 456 } 457 z /= nx + f; 458 nx += 1.0L; 459 x = nx + f; 460 } 461 z = fabs(z); 462 if ( x == 2.0L ) 463 return log(z); 464 x = (nx - 2.0L) + f; 465 real p = x * rationalPoly( x, logGammaNumerator, logGammaDenominator); 466 return log(z) + p; 467 } 468 469 // const real MAXLGM = 1.04848146839019521116e+4928L; 470 // if( x > MAXLGM ) return sgngaml * real.infinity; 471 472 const real LOGSQRT2PI = 0.91893853320467274178L; // log( sqrt( 2*pi ) ) 473 474 q = ( x - 0.5L ) * log(x) - x + LOGSQRT2PI; 475 if (x > 1.0e10L) return q; 476 real p = 1.0L / (x*x); 477 q += poly( p, logGammaStirlingCoeffs ) / x; 478 return q ; 479} 480 481unittest { 482 assert(isIdentical(logGamma(NaN(0xDEF)), NaN(0xDEF))); 483 assert(logGamma(real.infinity) == real.infinity); 484 assert(logGamma(-1.0) == real.infinity); 485 assert(logGamma(0.0) == real.infinity); 486 assert(logGamma(-50.0) == real.infinity); 487 assert(isIdentical(0.0L, logGamma(1.0L))); 488 assert(isIdentical(0.0L, logGamma(2.0L))); 489 assert(logGamma(real.min_normal*real.epsilon) == real.infinity); 490 assert(logGamma(-real.min_normal*real.epsilon) == real.infinity); 491 492 // x, correct loggamma(x), correct d/dx loggamma(x). 493 static real[] testpoints = [ 494 8.0L, 8.525146484375L + 1.48766904143001655310E-5, 2.01564147795560999654E0L, 495 8.99993896484375e-1L, 6.6375732421875e-2L + 5.11505711292524166220E-6L, -7.54938684259372234258E-1, 496 7.31597900390625e-1L, 2.2369384765625e-1 + 5.21506341809849792422E-6L, -1.13355566660398608343E0L, 497 2.31639862060546875e-1L, 1.3686676025390625L + 1.12609441752996145670E-5L, -4.56670961813812679012E0, 498 1.73162841796875L, -8.88214111328125e-2L + 3.36207740803753034508E-6L, 2.33339034686200586920E-1L, 499 1.23162841796875L, -9.3902587890625e-2L + 1.28765089229009648104E-5L, -2.49677345775751390414E-1L, 500 7.3786976294838206464e19L, 3.301798506038663053312e21L - 1.656137564136932662487046269677E5L, 501 4.57477139169563904215E1L, 502 1.08420217248550443401E-19L, 4.36682586669921875e1L + 1.37082843669932230418E-5L, 503 -9.22337203685477580858E18L, 504 1.0L, 0.0L, -5.77215664901532860607E-1L, 505 2.0L, 0.0L, 4.22784335098467139393E-1L, 506 -0.5L, 1.2655029296875L + 9.19379714539648894580E-6L, 3.64899739785765205590E-2L, 507 -1.5L, 8.6004638671875e-1L + 6.28657731014510932682E-7L, 7.03156640645243187226E-1L, 508 -2.5L, -5.6243896484375E-2L + 1.79986700949327405470E-7, 1.10315664064524318723E0L, 509 -3.5L, -1.30902099609375L + 1.43111007079536392848E-5L, 1.38887092635952890151E0L 510 ]; 511 // TODO: test derivatives as well. 512 for (int i=0; i<testpoints.length; i+=3) { 513 assert( feqrel(logGamma(testpoints[i]), testpoints[i+1]) > real.mant_dig-5); 514 if (testpoints[i]<MAXGAMMA) { 515 assert( feqrel(log(fabs(gamma(testpoints[i]))), testpoints[i+1]) > real.mant_dig-5); 516 } 517 } 518 assert(logGamma(-50.2) == log(fabs(gamma(-50.2)))); 519 assert(logGamma(-0.008) == log(fabs(gamma(-0.008)))); 520 assert(feqrel(logGamma(-38.8),log(fabs(gamma(-38.8)))) > real.mant_dig-4); 521 assert(feqrel(logGamma(1500.0L),log(gamma(1500.0L))) > real.mant_dig-2); 522} 523 524 525private { 526enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max) 527enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min*real.epsilon) = log(smallest denormal) 528enum real BETA_BIG = 9.223372036854775808e18L; 529enum real BETA_BIGINV = 1.084202172485504434007e-19L; 530} 531 532/** Incomplete beta integral 533 * 534 * Returns incomplete beta integral of the arguments, evaluated 535 * from zero to x. The regularized incomplete beta function is defined as 536 * 537 * betaIncomplete(a, b, x) = Γ(a+b)/(Γ(a) Γ(b)) * 538 * $(INTEGRATE 0, x) $(POWER t, a-1)$(POWER (1-t),b-1) dt 539 * 540 * and is the same as the the cumulative distribution function. 541 * 542 * The domain of definition is 0 <= x <= 1. In this 543 * implementation a and b are restricted to positive values. 544 * The integral from x to 1 may be obtained by the symmetry 545 * relation 546 * 547 * betaIncompleteCompl(a, b, x ) = betaIncomplete( b, a, 1-x ) 548 * 549 * The integral is evaluated by a continued fraction expansion 550 * or, when b*x is small, by a power series. 551 */ 552real betaIncomplete(real aa, real bb, real xx ) 553{ 554 if ( !(aa>0 && bb>0) ) 555 { 556 if ( isNaN(aa) ) return aa; 557 if ( isNaN(bb) ) return bb; 558 return real.nan; // domain error 559 } 560 if (!(xx>0 && xx<1.0)) { 561 if (isNaN(xx)) return xx; 562 if ( xx == 0.0L ) return 0.0; 563 if ( xx == 1.0L ) return 1.0; 564 return real.nan; // domain error 565 } 566 if ( (bb * xx) <= 1.0L && xx <= 0.95L) { 567 return betaDistPowerSeries(aa, bb, xx); 568 } 569 real x; 570 real xc; // = 1 - x 571 572 real a, b; 573 int flag = 0; 574 575 /* Reverse a and b if x is greater than the mean. */ 576 if( xx > (aa/(aa+bb)) ) { 577 // here x > aa/(aa+bb) and (bb*x>1 or x>0.95) 578 flag = 1; 579 a = bb; 580 b = aa; 581 xc = xx; 582 x = 1.0L - xx; 583 } else { 584 a = aa; 585 b = bb; 586 xc = 1.0L - xx; 587 x = xx; 588 } 589 590 if( flag == 1 && (b * x) <= 1.0L && x <= 0.95L) { 591 // here xx > aa/(aa+bb) and ((bb*xx>1) or xx>0.95) and (aa*(1-xx)<=1) and xx > 0.05 592 return 1.0 - betaDistPowerSeries(a, b, x); // note loss of precision 593 } 594 595 real w; 596 // Choose expansion for optimal convergence 597 // One is for x * (a+b+2) < (a+1), 598 // the other is for x * (a+b+2) > (a+1). 599 real y = x * (a+b-2.0L) - (a-1.0L); 600 if( y < 0.0L ) { 601 w = betaDistExpansion1( a, b, x ); 602 } else { 603 w = betaDistExpansion2( a, b, x ) / xc; 604 } 605 606 /* Multiply w by the factor 607 a b 608 x (1-x) Gamma(a+b) / ( a Gamma(a) Gamma(b) ) . */ 609 610 y = a * log(x); 611 real t = b * log(xc); 612 if ( (a+b) < MAXGAMMA && fabs(y) < MAXLOG && fabs(t) < MAXLOG ) { 613 t = pow(xc,b); 614 t *= pow(x,a); 615 t /= a; 616 t *= w; 617 t *= gamma(a+b) / (gamma(a) * gamma(b)); 618 } else { 619 /* Resort to logarithms. */ 620 y += t + logGamma(a+b) - logGamma(a) - logGamma(b); 621 y += log(w/a); 622 623 t = exp(y); 624/+ 625 // There seems to be a bug in Cephes at this point. 626 // Problems occur for y > MAXLOG, not y < MINLOG. 627 if( y < MINLOG ) { 628 t = 0.0L; 629 } else { 630 t = exp(y); 631 } 632+/ 633 } 634 if( flag == 1 ) { 635/+ // CEPHES includes this code, but I think it is erroneous. 636 if( t <= real.epsilon ) { 637 t = 1.0L - real.epsilon; 638 } else 639+/ 640 t = 1.0L - t; 641 } 642 return t; 643} 644 645/** Inverse of incomplete beta integral 646 * 647 * Given y, the function finds x such that 648 * 649 * betaIncomplete(a, b, x) == y 650 * 651 * Newton iterations or interval halving is used. 652 */ 653real betaIncompleteInv(real aa, real bb, real yy0 ) 654{ 655 real a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt; 656 int i, rflg, dir, nflg; 657 658 if (isNaN(yy0)) return yy0; 659 if (isNaN(aa)) return aa; 660 if (isNaN(bb)) return bb; 661 if( yy0 <= 0.0L ) 662 return 0.0L; 663 if( yy0 >= 1.0L ) 664 return 1.0L; 665 x0 = 0.0L; 666 yl = 0.0L; 667 x1 = 1.0L; 668 yh = 1.0L; 669 if( aa <= 1.0L || bb <= 1.0L ) { 670 dithresh = 1.0e-7L; 671 rflg = 0; 672 a = aa; 673 b = bb; 674 y0 = yy0; 675 x = a/(a+b); 676 y = betaIncomplete( a, b, x ); 677 nflg = 0; 678 goto ihalve; 679 } else { 680 nflg = 0; 681 dithresh = 1.0e-4L; 682 } 683 684 // approximation to inverse function 685 686 yp = -normalDistributionInvImpl( yy0 ); 687 688 if( yy0 > 0.5L ) { 689 rflg = 1; 690 a = bb; 691 b = aa; 692 y0 = 1.0L - yy0; 693 yp = -yp; 694 } else { 695 rflg = 0; 696 a = aa; 697 b = bb; 698 y0 = yy0; 699 } 700 701 lgm = (yp * yp - 3.0L)/6.0L; 702 x = 2.0L/( 1.0L/(2.0L * a-1.0L) + 1.0L/(2.0L * b - 1.0L) ); 703 d = yp * sqrt( x + lgm ) / x 704 - ( 1.0L/(2.0L * b - 1.0L) - 1.0L/(2.0L * a - 1.0L) ) 705 * (lgm + (5.0L/6.0L) - 2.0L/(3.0L * x)); 706 d = 2.0L * d; 707 if( d < MINLOG ) { 708 x = 1.0L; 709 goto under; 710 } 711 x = a/( a + b * exp(d) ); 712 y = betaIncomplete( a, b, x ); 713 yp = (y - y0)/y0; 714 if( fabs(yp) < 0.2 ) 715 goto newt; 716 717 /* Resort to interval halving if not close enough. */ 718ihalve: 719 720 dir = 0; 721 di = 0.5L; 722 for( i=0; i<400; i++ ) { 723 if( i != 0 ) { 724 x = x0 + di * (x1 - x0); 725 if( x == 1.0L ) { 726 x = 1.0L - real.epsilon; 727 } 728 if( x == 0.0L ) { 729 di = 0.5; 730 x = x0 + di * (x1 - x0); 731 if( x == 0.0 ) 732 goto under; 733 } 734 y = betaIncomplete( a, b, x ); 735 yp = (x1 - x0)/(x1 + x0); 736 if( fabs(yp) < dithresh ) 737 goto newt; 738 yp = (y-y0)/y0; 739 if( fabs(yp) < dithresh ) 740 goto newt; 741 } 742 if( y < y0 ) { 743 x0 = x; 744 yl = y; 745 if( dir < 0 ) { 746 dir = 0; 747 di = 0.5L; 748 } else if( dir > 3 ) 749 di = 1.0L - (1.0L - di) * (1.0L - di); 750 else if( dir > 1 ) 751 di = 0.5L * di + 0.5L; 752 else 753 di = (y0 - y)/(yh - yl); 754 dir += 1; 755 if( x0 > 0.95L ) { 756 if( rflg == 1 ) { 757 rflg = 0; 758 a = aa; 759 b = bb; 760 y0 = yy0; 761 } else { 762 rflg = 1; 763 a = bb; 764 b = aa; 765 y0 = 1.0 - yy0; 766 } 767 x = 1.0L - x; 768 y = betaIncomplete( a, b, x ); 769 x0 = 0.0; 770 yl = 0.0; 771 x1 = 1.0; 772 yh = 1.0; 773 goto ihalve; 774 } 775 } else { 776 x1 = x; 777 if( rflg == 1 && x1 < real.epsilon ) { 778 x = 0.0L; 779 goto done; 780 } 781 yh = y; 782 if( dir > 0 ) { 783 dir = 0; 784 di = 0.5L; 785 } 786 else if( dir < -3 ) 787 di = di * di; 788 else if( dir < -1 ) 789 di = 0.5L * di; 790 else 791 di = (y - y0)/(yh - yl); 792 dir -= 1; 793 } 794 } 795 if( x0 >= 1.0L ) { 796 // partial loss of precision 797 x = 1.0L - real.epsilon; 798 goto done; 799 } 800 if( x <= 0.0L ) { 801under: 802 // underflow has occurred 803 x = real.min_normal * real.min_normal; 804 goto done; 805 } 806 807newt: 808 809 if ( nflg ) { 810 goto done; 811 } 812 nflg = 1; 813 lgm = logGamma(a+b) - logGamma(a) - logGamma(b); 814 815 for( i=0; i<15; i++ ) { 816 /* Compute the function at this point. */ 817 if ( i != 0 ) 818 y = betaIncomplete(a,b,x); 819 if ( y < yl ) { 820 x = x0; 821 y = yl; 822 } else if( y > yh ) { 823 x = x1; 824 y = yh; 825 } else if( y < y0 ) { 826 x0 = x; 827 yl = y; 828 } else { 829 x1 = x; 830 yh = y; 831 } 832 if( x == 1.0L || x == 0.0L ) 833 break; 834 /* Compute the derivative of the function at this point. */ 835 d = (a - 1.0L) * log(x) + (b - 1.0L) * log(1.0L - x) + lgm; 836 if ( d < MINLOG ) { 837 goto done; 838 } 839 if ( d > MAXLOG ) { 840 break; 841 } 842 d = exp(d); 843 /* Compute the step to the next approximation of x. */ 844 d = (y - y0)/d; 845 xt = x - d; 846 if ( xt <= x0 ) { 847 y = (x - x0) / (x1 - x0); 848 xt = x0 + 0.5L * y * (x - x0); 849 if( xt <= 0.0L ) 850 break; 851 } 852 if ( xt >= x1 ) { 853 y = (x1 - x) / (x1 - x0); 854 xt = x1 - 0.5L * y * (x1 - x); 855 if ( xt >= 1.0L ) 856 break; 857 } 858 x = xt; 859 if ( fabs(d/x) < (128.0L * real.epsilon) ) 860 goto done; 861 } 862 /* Did not converge. */ 863 dithresh = 256.0L * real.epsilon; 864 goto ihalve; 865 866done: 867 if ( rflg ) { 868 if( x <= real.epsilon ) 869 x = 1.0L - real.epsilon; 870 else 871 x = 1.0L - x; 872 } 873 return x; 874} 875 876unittest { // also tested by the normal distribution 877 // check NaN propagation 878 assert(isIdentical(betaIncomplete(NaN(0xABC),2,3), NaN(0xABC))); 879 assert(isIdentical(betaIncomplete(7,NaN(0xABC),3), NaN(0xABC))); 880 assert(isIdentical(betaIncomplete(7,15,NaN(0xABC)), NaN(0xABC))); 881 assert(isIdentical(betaIncompleteInv(NaN(0xABC),1,17), NaN(0xABC))); 882 assert(isIdentical(betaIncompleteInv(2,NaN(0xABC),8), NaN(0xABC))); 883 assert(isIdentical(betaIncompleteInv(2,3, NaN(0xABC)), NaN(0xABC))); 884 885 assert(isNaN(betaIncomplete(-1, 2, 3))); 886 887 assert(betaIncomplete(1, 2, 0)==0); 888 assert(betaIncomplete(1, 2, 1)==1); 889 assert(isNaN(betaIncomplete(1, 2, 3))); 890 assert(betaIncompleteInv(1, 1, 0)==0); 891 assert(betaIncompleteInv(1, 1, 1)==1); 892 893 // Test against Mathematica betaRegularized[z,a,b] 894 // These arbitrary points are chosen to give good code coverage. 895 assert(feqrel(betaIncomplete(8, 10, 0.2), 0.010_934_315_234_099_2L) >= real.mant_dig - 5); 896 assert(feqrel(betaIncomplete(2, 2.5, 0.9),0.989_722_597_604_452_767_171_003_59L) >= real.mant_dig - 1 ); 897 assert(feqrel(betaIncomplete(1000, 800, 0.5), 1.179140859734704555102808541457164E-06L) >= real.mant_dig - 13 ); 898 assert(feqrel(betaIncomplete(0.0001, 10000, 0.0001),0.999978059362107134278786L) >= real.mant_dig - 18 ); 899 assert(betaIncomplete(0.01, 327726.7, 0.545113) == 1.0); 900 assert(feqrel(betaIncompleteInv(8, 10, 0.010_934_315_234_099_2L), 0.2L) >= real.mant_dig - 2); 901 assert(feqrel(betaIncomplete(0.01, 498.437, 0.0121433),0.99999664562033077636065L) >= real.mant_dig - 1); 902 assert(feqrel(betaIncompleteInv(5, 10, 0.2000002972865658842), 0.229121208190918L) >= real.mant_dig - 3); 903 assert(feqrel(betaIncompleteInv(4, 7, 0.8000002209179505L), 0.483657360076904L) >= real.mant_dig - 3); 904 905 // Coverage tests. I don't have correct values for these tests, but 906 // these values cover most of the code, so they are useful for 907 // regression testing. 908 // Extensive testing failed to increase the coverage. It seems likely that about 909 // half the code in this function is unnecessary; there is potential for 910 // significant improvement over the original CEPHES code. 911 912 assert(betaIncompleteInv(0.01, 8e-48, 5.45464e-20)==1-real.epsilon); 913 assert(betaIncompleteInv(0.01, 8e-48, 9e-26)==1-real.epsilon); 914 915 // Beware: a one-bit change in pow() changes almost all digits in the result! 916 assert(feqrel(betaIncompleteInv(0x1.b3d151fbba0eb18p+1, 1.2265e-19, 2.44859e-18),0x1.c0110c8531d0952cp-1L) > 10); 917 // This next case uncovered a one-bit difference in the FYL2X instruction 918 // between Intel and AMD processors. This difference gets magnified by 2^^38. 919 // WolframAlpha crashes attempting to calculate this. 920 assert(feqrel(betaIncompleteInv(0x1.ff1275ae5b939bcap-41, 4.6713e18, 0.0813601), 921 0x1.f97749d90c7adba8p-63L) >= real.mant_dig - 39); 922 real a1 = 3.40483; 923 assert(betaIncompleteInv(a1, 4.0640301659679627772e19L, 0.545113)== 0x1.ba8c08108aaf5d14p-109); 924 real b1 = 2.82847e-25; 925 assert(feqrel(betaIncompleteInv(0.01, b1, 9e-26), 0x1.549696104490aa9p-830L) >= real.mant_dig-10); 926 927 // --- Problematic cases --- 928 // This is a situation where the series expansion fails to converge 929 assert( isNaN(betaIncompleteInv(0.12167, 4.0640301659679627772e19L, 0.0813601))); 930 // This next result is almost certainly erroneous. 931 // Mathematica states: "(cannot be determined by current methods)" 932 assert(betaIncomplete(1.16251e20, 2.18e39, 5.45e-20)==-real.infinity); 933 // WolframAlpha gives no result for this, though indicates that it approximately 1.0 - 1.3e-9 934 assert(1- betaIncomplete(0.01, 328222, 4.0375e-5) == 0x1.5f62926b4p-30); 935} 936 937 938private { 939// Implementation functions 940 941// Continued fraction expansion #1 for incomplete beta integral 942// Use when x < (a+1)/(a+b+2) 943real betaDistExpansion1(real a, real b, real x ) 944{ 945 real xk, pk, pkm1, pkm2, qk, qkm1, qkm2; 946 real k1, k2, k3, k4, k5, k6, k7, k8; 947 real r, t, ans; 948 int n; 949 950 k1 = a; 951 k2 = a + b; 952 k3 = a; 953 k4 = a + 1.0L; 954 k5 = 1.0L; 955 k6 = b - 1.0L; 956 k7 = k4; 957 k8 = a + 2.0L; 958 959 pkm2 = 0.0L; 960 qkm2 = 1.0L; 961 pkm1 = 1.0L; 962 qkm1 = 1.0L; 963 ans = 1.0L; 964 r = 1.0L; 965 n = 0; 966 const real thresh = 3.0L * real.epsilon; 967 do { 968 xk = -( x * k1 * k2 )/( k3 * k4 ); 969 pk = pkm1 + pkm2 * xk; 970 qk = qkm1 + qkm2 * xk; 971 pkm2 = pkm1; 972 pkm1 = pk; 973 qkm2 = qkm1; 974 qkm1 = qk; 975 976 xk = ( x * k5 * k6 )/( k7 * k8 ); 977 pk = pkm1 + pkm2 * xk; 978 qk = qkm1 + qkm2 * xk; 979 pkm2 = pkm1; 980 pkm1 = pk; 981 qkm2 = qkm1; 982 qkm1 = qk; 983 984 if( qk != 0.0L ) 985 r = pk/qk; 986 if( r != 0.0L ) { 987 t = fabs( (ans - r)/r ); 988 ans = r; 989 } else { 990 t = 1.0L; 991 } 992 993 if( t < thresh ) 994 return ans; 995 996 k1 += 1.0L; 997 k2 += 1.0L; 998 k3 += 2.0L; 999 k4 += 2.0L; 1000 k5 += 1.0L; 1001 k6 -= 1.0L; 1002 k7 += 2.0L; 1003 k8 += 2.0L; 1004 1005 if( (fabs(qk) + fabs(pk)) > BETA_BIG ) { 1006 pkm2 *= BETA_BIGINV; 1007 pkm1 *= BETA_BIGINV; 1008 qkm2 *= BETA_BIGINV; 1009 qkm1 *= BETA_BIGINV; 1010 } 1011 if( (fabs(qk) < BETA_BIGINV) || (fabs(pk) < BETA_BIGINV) ) { 1012 pkm2 *= BETA_BIG; 1013 pkm1 *= BETA_BIG; 1014 qkm2 *= BETA_BIG; 1015 qkm1 *= BETA_BIG; 1016 } 1017 } 1018 while( ++n < 400 ); 1019// loss of precision has occurred 1020// mtherr( "incbetl", PLOSS ); 1021 return ans; 1022} 1023 1024// Continued fraction expansion #2 for incomplete beta integral 1025// Use when x > (a+1)/(a+b+2) 1026real betaDistExpansion2(real a, real b, real x ) 1027{ 1028 real xk, pk, pkm1, pkm2, qk, qkm1, qkm2; 1029 real k1, k2, k3, k4, k5, k6, k7, k8; 1030 real r, t, ans, z; 1031 1032 k1 = a; 1033 k2 = b - 1.0L; 1034 k3 = a; 1035 k4 = a + 1.0L; 1036 k5 = 1.0L; 1037 k6 = a + b; 1038 k7 = a + 1.0L; 1039 k8 = a + 2.0L; 1040 1041 pkm2 = 0.0L; 1042 qkm2 = 1.0L; 1043 pkm1 = 1.0L; 1044 qkm1 = 1.0L; 1045 z = x / (1.0L-x); 1046 ans = 1.0L; 1047 r = 1.0L; 1048 int n = 0; 1049 const real thresh = 3.0L * real.epsilon; 1050 do { 1051 1052 xk = -( z * k1 * k2 )/( k3 * k4 ); 1053 pk = pkm1 + pkm2 * xk; 1054 qk = qkm1 + qkm2 * xk; 1055 pkm2 = pkm1; 1056 pkm1 = pk; 1057 qkm2 = qkm1; 1058 qkm1 = qk; 1059 1060 xk = ( z * k5 * k6 )/( k7 * k8 ); 1061 pk = pkm1 + pkm2 * xk; 1062 qk = qkm1 + qkm2 * xk; 1063 pkm2 = pkm1; 1064 pkm1 = pk; 1065 qkm2 = qkm1; 1066 qkm1 = qk; 1067 1068 if( qk != 0.0L ) 1069 r = pk/qk; 1070 if( r != 0.0L ) { 1071 t = fabs( (ans - r)/r ); 1072 ans = r; 1073 } else 1074 t = 1.0L; 1075 1076 if( t < thresh ) 1077 return ans; 1078 k1 += 1.0L; 1079 k2 -= 1.0L; 1080 k3 += 2.0L; 1081 k4 += 2.0L; 1082 k5 += 1.0L; 1083 k6 += 1.0L; 1084 k7 += 2.0L; 1085 k8 += 2.0L; 1086 1087 if( (fabs(qk) + fabs(pk)) > BETA_BIG ) { 1088 pkm2 *= BETA_BIGINV; 1089 pkm1 *= BETA_BIGINV; 1090 qkm2 *= BETA_BIGINV; 1091 qkm1 *= BETA_BIGINV; 1092 } 1093 if( (fabs(qk) < BETA_BIGINV) || (fabs(pk) < BETA_BIGINV) ) { 1094 pkm2 *= BETA_BIG; 1095 pkm1 *= BETA_BIG; 1096 qkm2 *= BETA_BIG; 1097 qkm1 *= BETA_BIG; 1098 } 1099 } while( ++n < 400 ); 1100// loss of precision has occurred 1101//mtherr( "incbetl", PLOSS ); 1102 return ans; 1103} 1104 1105/* Power series for incomplete gamma integral. 1106 Use when b*x is small. */ 1107real betaDistPowerSeries(real a, real b, real x ) 1108{ 1109 real ai = 1.0L / a; 1110 real u = (1.0L - b) * x; 1111 real v = u / (a + 1.0L); 1112 real t1 = v; 1113 real t = u; 1114 real n = 2.0L; 1115 real s = 0.0L; 1116 real z = real.epsilon * ai; 1117 while( fabs(v) > z ) { 1118 u = (n - b) * x / n; 1119 t *= u; 1120 v = t / (a + n); 1121 s += v; 1122 n += 1.0L; 1123 } 1124 s += t1; 1125 s += ai; 1126 1127 u = a * log(x); 1128 if ( (a+b) < MAXGAMMA && fabs(u) < MAXLOG ) { 1129 t = gamma(a+b)/(gamma(a)*gamma(b)); 1130 s = s * t * pow(x,a); 1131 } else { 1132 t = logGamma(a+b) - logGamma(a) - logGamma(b) + u + log(s); 1133 1134 if( t < MINLOG ) { 1135 s = 0.0L; 1136 } else 1137 s = exp(t); 1138 } 1139 return s; 1140} 1141 1142} 1143 1144/*************************************** 1145 * Incomplete gamma integral and its complement 1146 * 1147 * These functions are defined by 1148 * 1149 * gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a) 1150 * 1151 * gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x) 1152 * = ($(INTEGRATE x, ∞) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a) 1153 * 1154 * In this implementation both arguments must be positive. 1155 * The integral is evaluated by either a power series or 1156 * continued fraction expansion, depending on the relative 1157 * values of a and x. 1158 */ 1159real gammaIncomplete(real a, real x ) 1160in { 1161 assert(x >= 0); 1162 assert(a > 0); 1163} 1164body { 1165 /* left tail of incomplete gamma function: 1166 * 1167 * inf. k 1168 * a -x - x 1169 * x e > ---------- 1170 * - - 1171 * k=0 | (a+k+1) 1172 * 1173 */ 1174 if (x==0) 1175 return 0.0L; 1176 1177 if ( (x > 1.0L) && (x > a ) ) 1178 return 1.0L - gammaIncompleteCompl(a,x); 1179 1180 real ax = a * log(x) - x - logGamma(a); 1181/+ 1182 if( ax < MINLOGL ) return 0; // underflow 1183 // { mtherr( "igaml", UNDERFLOW ); return( 0.0L ); } 1184+/ 1185 ax = exp(ax); 1186 1187 /* power series */ 1188 real r = a; 1189 real c = 1.0L; 1190 real ans = 1.0L; 1191 1192 do { 1193 r += 1.0L; 1194 c *= x/r; 1195 ans += c; 1196 } while( c/ans > real.epsilon ); 1197 1198 return ans * ax/a; 1199} 1200 1201/** ditto */ 1202real gammaIncompleteCompl(real a, real x ) 1203in { 1204 assert(x >= 0); 1205 assert(a > 0); 1206} 1207body { 1208 if (x==0) 1209 return 1.0L; 1210 if ( (x < 1.0L) || (x < a) ) 1211 return 1.0L - gammaIncomplete(a,x); 1212 1213 // DAC (Cephes bug fix): This is necessary to avoid 1214 // spurious nans, eg 1215 // log(x)-x = NaN when x = real.infinity 1216 const real MAXLOGL = 1.1356523406294143949492E4L; 1217 if (x > MAXLOGL) 1218 return igammaTemmeLarge(a, x); 1219 1220 real ax = a * log(x) - x - logGamma(a); 1221//const real MINLOGL = -1.1355137111933024058873E4L; 1222// if ( ax < MINLOGL ) return 0; // underflow; 1223 ax = exp(ax); 1224 1225 1226 /* continued fraction */ 1227 real y = 1.0L - a; 1228 real z = x + y + 1.0L; 1229 real c = 0.0L; 1230 1231 real pk, qk, t; 1232 1233 real pkm2 = 1.0L; 1234 real qkm2 = x; 1235 real pkm1 = x + 1.0L; 1236 real qkm1 = z * x; 1237 real ans = pkm1/qkm1; 1238 1239 do { 1240 c += 1.0L; 1241 y += 1.0L; 1242 z += 2.0L; 1243 real yc = y * c; 1244 pk = pkm1 * z - pkm2 * yc; 1245 qk = qkm1 * z - qkm2 * yc; 1246 if( qk != 0.0L ) { 1247 real r = pk/qk; 1248 t = fabs( (ans - r)/r ); 1249 ans = r; 1250 } else { 1251 t = 1.0L; 1252 } 1253 pkm2 = pkm1; 1254 pkm1 = pk; 1255 qkm2 = qkm1; 1256 qkm1 = qk; 1257 1258 const real BIG = 9.223372036854775808e18L; 1259 1260 if ( fabs(pk) > BIG ) { 1261 pkm2 /= BIG; 1262 pkm1 /= BIG; 1263 qkm2 /= BIG; 1264 qkm1 /= BIG; 1265 } 1266 } while ( t > real.epsilon ); 1267 1268 return ans * ax; 1269} 1270 1271/** Inverse of complemented incomplete gamma integral 1272 * 1273 * Given a and p, the function finds x such that 1274 * 1275 * gammaIncompleteCompl( a, x ) = p. 1276 * 1277 * Starting with the approximate value x = a $(POWER t, 3), where 1278 * t = 1 - d - normalDistributionInv(p) sqrt(d), 1279 * and d = 1/9a, 1280 * the routine performs up to 10 Newton iterations to find the 1281 * root of incompleteGammaCompl(a,x) - p = 0. 1282 */ 1283real gammaIncompleteComplInv(real a, real p) 1284in { 1285 assert(p>=0 && p<= 1); 1286 assert(a>0); 1287} 1288body { 1289 if (p==0) return real.infinity; 1290 1291 real y0 = p; 1292 const real MAXLOGL = 1.1356523406294143949492E4L; 1293 real x0, x1, x, yl, yh, y, d, lgm, dithresh; 1294 int i, dir; 1295 1296 /* bound the solution */ 1297 x0 = real.max; 1298 yl = 0.0L; 1299 x1 = 0.0L; 1300 yh = 1.0L; 1301 dithresh = 4.0 * real.epsilon; 1302 1303 /* approximation to inverse function */ 1304 d = 1.0L/(9.0L*a); 1305 y = 1.0L - d - normalDistributionInvImpl(y0) * sqrt(d); 1306 x = a * y * y * y; 1307 1308 lgm = logGamma(a); 1309 1310 for( i=0; i<10; i++ ) { 1311 if( x > x0 || x < x1 ) 1312 goto ihalve; 1313 y = gammaIncompleteCompl(a,x); 1314 if ( y < yl || y > yh ) 1315 goto ihalve; 1316 if ( y < y0 ) { 1317 x0 = x; 1318 yl = y; 1319 } else { 1320 x1 = x; 1321 yh = y; 1322 } 1323 /* compute the derivative of the function at this point */ 1324 d = (a - 1.0L) * log(x0) - x0 - lgm; 1325 if ( d < -MAXLOGL ) 1326 goto ihalve; 1327 d = -exp(d); 1328 /* compute the step to the next approximation of x */ 1329 d = (y - y0)/d; 1330 x = x - d; 1331 if ( i < 3 ) continue; 1332 if ( fabs(d/x) < dithresh ) return x; 1333 } 1334 1335 /* Resort to interval halving if Newton iteration did not converge. */ 1336ihalve: 1337 d = 0.0625L; 1338 if ( x0 == real.max ) { 1339 if( x <= 0.0L ) 1340 x = 1.0L; 1341 while( x0 == real.max ) { 1342 x = (1.0L + d) * x; 1343 y = gammaIncompleteCompl( a, x ); 1344 if ( y < y0 ) { 1345 x0 = x; 1346 yl = y; 1347 break; 1348 } 1349 d = d + d; 1350 } 1351 } 1352 d = 0.5L; 1353 dir = 0; 1354 1355 for( i=0; i<400; i++ ) { 1356 x = x1 + d * (x0 - x1); 1357 y = gammaIncompleteCompl( a, x ); 1358 lgm = (x0 - x1)/(x1 + x0); 1359 if ( fabs(lgm) < dithresh ) 1360 break; 1361 lgm = (y - y0)/y0; 1362 if ( fabs(lgm) < dithresh ) 1363 break; 1364 if ( x <= 0.0L ) 1365 break; 1366 if ( y > y0 ) { 1367 x1 = x; 1368 yh = y; 1369 if ( dir < 0 ) { 1370 dir = 0; 1371 d = 0.5L; 1372 } else if ( dir > 1 ) 1373 d = 0.5L * d + 0.5L; 1374 else 1375 d = (y0 - yl)/(yh - yl); 1376 dir += 1; 1377 } else { 1378 x0 = x; 1379 yl = y; 1380 if ( dir > 0 ) { 1381 dir = 0; 1382 d = 0.5L; 1383 } else if ( dir < -1 ) 1384 d = 0.5L * d; 1385 else 1386 d = (y0 - yl)/(yh - yl); 1387 dir -= 1; 1388 } 1389 } 1390 /+ 1391 if( x == 0.0L ) 1392 mtherr( "igamil", UNDERFLOW ); 1393 +/ 1394 return x; 1395} 1396 1397unittest { 1398//Values from Excel's GammaInv(1-p, x, 1) 1399assert(fabs(gammaIncompleteComplInv(1, 0.5) - 0.693147188044814) < 0.00000005); 1400assert(fabs(gammaIncompleteComplInv(12, 0.99) - 5.42818075054289) < 0.00000005); 1401assert(fabs(gammaIncompleteComplInv(100, 0.8) - 91.5013985848288L) < 0.000005); 1402assert(gammaIncomplete(1, 0)==0); 1403assert(gammaIncompleteCompl(1, 0)==1); 1404assert(gammaIncomplete(4545, real.infinity)==1); 1405 1406// Values from Excel's (1-GammaDist(x, alpha, 1, TRUE)) 1407 1408assert(fabs(1.0L-gammaIncompleteCompl(0.5, 2) - 0.954499729507309L) < 0.00000005); 1409assert(fabs(gammaIncomplete(0.5, 2) - 0.954499729507309L) < 0.00000005); 1410// Fixed Cephes bug: 1411assert(gammaIncompleteCompl(384, real.infinity)==0); 1412assert(gammaIncompleteComplInv(3, 0)==real.infinity); 1413// Fixed a bug that caused gammaIncompleteCompl to return a wrong value when 1414// x was larger than a, but not by much, and both were large: 1415// The value is from WolframAlpha (Gamma[100000, 100001, inf] / Gamma[100000]) 1416assert(fabs(gammaIncompleteCompl(100000, 100001) - 0.49831792109) < 0.000000000005); 1417} 1418 1419 1420/** Digamma function 1421* 1422* The digamma function is the logarithmic derivative of the gamma function. 1423* 1424* digamma(x) = d/dx logGamma(x) 1425* 1426*/ 1427real digamma(real x) 1428{ 1429 // Based on CEPHES, Stephen L. Moshier. 1430 1431 // DAC: These values are Bn / n for n=2,4,6,8,10,12,14. 1432 static immutable real [7] Bn_n = [ 1433 1.0L/(6*2), -1.0L/(30*4), 1.0L/(42*6), -1.0L/(30*8), 1434 5.0L/(66*10), -691.0L/(2730*12), 7.0L/(6*14) ]; 1435 1436 real p, q, nz, s, w, y, z; 1437 long i, n; 1438 int negative; 1439 1440 negative = 0; 1441 nz = 0.0; 1442 1443 if ( x <= 0.0 ) { 1444 negative = 1; 1445 q = x; 1446 p = floor(q); 1447 if( p == q ) { 1448 return real.nan; // singularity. 1449 } 1450 /* Remove the zeros of tan(PI x) 1451 * by subtracting the nearest integer from x 1452 */ 1453 nz = q - p; 1454 if ( nz != 0.5 ) { 1455 if ( nz > 0.5 ) { 1456 p += 1.0; 1457 nz = q - p; 1458 } 1459 nz = PI/tan(PI*nz); 1460 } else { 1461 nz = 0.0; 1462 } 1463 x = 1.0 - x; 1464 } 1465 1466 // check for small positive integer 1467 if ((x <= 13.0) && (x == floor(x)) ) { 1468 y = 0.0; 1469 n = lrint(x); 1470 // DAC: CEPHES bugfix. Cephes did this in reverse order, which 1471 // created a larger roundoff error. 1472 for (i=n-1; i>0; --i) { 1473 y+=1.0L/i; 1474 } 1475 y -= EULERGAMMA; 1476 goto done; 1477 } 1478 1479 s = x; 1480 w = 0.0; 1481 while ( s < 10.0 ) { 1482 w += 1.0/s; 1483 s += 1.0; 1484 } 1485 1486 if ( s < 1.0e17 ) { 1487 z = 1.0/(s * s); 1488 y = z * poly(z, Bn_n); 1489 } else 1490 y = 0.0; 1491 1492 y = log(s) - 0.5L/s - y - w; 1493 1494done: 1495 if ( negative ) { 1496 y -= nz; 1497 } 1498 return y; 1499} 1500 1501version(unittest) import core.stdc.stdio; 1502unittest { 1503 // Exact values 1504 assert(digamma(1)== -EULERGAMMA); 1505 assert(feqrel(digamma(0.25), -PI/2 - 3* LN2 - EULERGAMMA) >= real.mant_dig-7); 1506 assert(feqrel(digamma(1.0L/6), -PI/2 *sqrt(3.0L) - 2* LN2 -1.5*log(3.0L) - EULERGAMMA) >= real.mant_dig-7); 1507 assert(digamma(-5).isNaN); 1508 assert(feqrel(digamma(2.5), -EULERGAMMA - 2*LN2 + 2.0 + 2.0L/3) >= real.mant_dig-9); 1509 assert(isIdentical(digamma(NaN(0xABC)), NaN(0xABC))); 1510 1511 for (int k=1; k<40; ++k) { 1512 real y=0; 1513 for (int u=k; u>=1; --u) { 1514 y += 1.0L/u; 1515 } 1516 assert(feqrel(digamma(k+1), -EULERGAMMA + y) >= real.mant_dig-2); 1517 } 1518}