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