PageRenderTime 179ms CodeModel.GetById 13ms app.highlight 152ms RepoModel.GetById 2ms app.codeStats 0ms

/std/internal/math/gammafunction.d

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