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