libphobos/src/std/internal/math/gammafunction.d

   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 =  &#915;
  17  *  INTEGRATE = $(BIG &#8747;<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 &pi;) 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, &infin;) $(POWER t, z-1)$(POWER e, -t) dt
 271  *
 272  *  $(TABLE_SV
 273  *    $(SVH  x,          $(GAMMA)(x) )
 274  *    $(SV  $(NAN),      $(NAN)      )
 275  *    $(SV  &plusmn;0.0, &plusmn;&infin;)
 276  *    $(SV integer > 0,  (x-1)!      )
 277  *    $(SV integer < 0,  $(NAN)      )
 278  *    $(SV +&infin;,     +&infin;    )
 279  *    $(SV -&infin;,     $(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,    +&infin;    )
 434  *    $(SV &plusmn;&infin;, +&infin;    )
 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) = &Gamma;(a+b)/(&Gamma;(a) &Gamma;(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, &infin;) $(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 }