/**
* Implementation of the gamma and beta functions, and their integrals.
*
* License: $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0).
* Copyright: Based on the CEPHES math library, which is
* Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com).
* Authors: Stephen L. Moshier (original C code). Conversion to D by Don Clugston
*
*
Macros:
* TABLE_SV = <table border="1" cellpadding="4" cellspacing="0">
* <caption>Special Values</caption>
* $0</table>
* SVH = $(TR $(TH $1) $(TH $2))
* SV = $(TR $(TD $1) $(TD $2))
* GAMMA = Γ
* INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>)
* POWER = $1<sup>$2</sup>
* NAN = $(RED NAN)
*/
module std.internal.math.gammafunction;
import std.internal.math.errorfunction;
import std.math;
import core.math : fabs, sin, sqrt;
pure:
nothrow:
@safe:
@nogc:
private {
enum real SQRT2PI = 2.50662827463100050242E0L; // sqrt(2pi)
immutable real EULERGAMMA = 0.57721_56649_01532_86060_65120_90082_40243_10421_59335_93992L; /** Euler-Mascheroni constant 0.57721566.. */
// Polynomial approximations for gamma and loggamma.
immutable real[8] GammaNumeratorCoeffs = [ 1.0L,
0x1.acf42d903366539ep-1L, 0x1.73a991c8475f1aeap-2L, 0x1.c7e918751d6b2a92p-4L,
0x1.86d162cca32cfe86p-6L, 0x1.0c378e2e6eaf7cd8p-8L, 0x1.dc5c66b7d05feb54p-12L,
0x1.616457b47e448694p-15L
];
immutable real[9] GammaDenominatorCoeffs = [ 1.0L,
0x1.a8f9faae5d8fc8bp-2L, -0x1.cb7895a6756eebdep-3L, -0x1.7b9bab006d30652ap-5L,
0x1.c671af78f312082ep-6L, -0x1.a11ebbfaf96252dcp-11L, -0x1.447b4d2230a77ddap-10L,
0x1.ec1d45bb85e06696p-13L,-0x1.d4ce24d05bd0a8e6p-17L
];
immutable real[9] GammaSmallCoeffs = [ 1.0L,
0x1.2788cfc6fb618f52p-1L, -0x1.4fcf4026afa2f7ecp-1L, -0x1.5815e8fa24d7e306p-5L,
0x1.5512320aea2ad71ap-3L, -0x1.59af0fb9d82e216p-5L, -0x1.3b4b61d3bfdf244ap-7L,
0x1.d9358e9d9d69fd34p-8L, -0x1.38fc4bcbada775d6p-10L
];
immutable real[9] GammaSmallNegCoeffs = [ -1.0L,
0x1.2788cfc6fb618f54p-1L, 0x1.4fcf4026afa2bc4cp-1L, -0x1.5815e8fa2468fec8p-5L,
-0x1.5512320baedaf4b6p-3L, -0x1.59af0fa283baf07ep-5L, 0x1.3b4a70de31e05942p-7L,
0x1.d9398be3bad13136p-8L, 0x1.291b73ee05bcbba2p-10L
];
immutable real[7] logGammaStirlingCoeffs = [
0x1.5555555555553f98p-4L, -0x1.6c16c16c07509b1p-9L, 0x1.a01a012461cbf1e4p-11L,
-0x1.3813089d3f9d164p-11L, 0x1.b911a92555a277b8p-11L, -0x1.ed0a7b4206087b22p-10L,
0x1.402523859811b308p-8L
];
immutable real[7] logGammaNumerator = [
-0x1.0edd25913aaa40a2p+23L, -0x1.31c6ce2e58842d1ep+24L, -0x1.f015814039477c3p+23L,
-0x1.74ffe40c4b184b34p+22L, -0x1.0d9c6d08f9eab55p+20L, -0x1.54c6b71935f1fc88p+16L,
-0x1.0e761b42932b2aaep+11L
];
immutable real[8] logGammaDenominator = [
-0x1.4055572d75d08c56p+24L, -0x1.deeb6013998e4d76p+24L, -0x1.106f7cded5dcc79ep+24L,
-0x1.25e17184848c66d2p+22L, -0x1.301303b99a614a0ap+19L, -0x1.09e76ab41ae965p+15L,
-0x1.00f95ced9e5f54eep+9L, 1.0L
];
/*
* Helper function: Gamma function computed by Stirling's formula.
*
* Stirling's formula for the gamma function is:
*
* $(GAMMA)(x) = sqrt(2 π) x<sup>x-0.5</sup> exp(-x) (1 + 1/x P(1/x))
*
*/
real gammaStirling(real x)
{
// CEPHES code Copyright 1994 by Stephen L. Moshier
static immutable real[9] SmallStirlingCoeffs = [
0x1.55555555555543aap-4L, 0x1.c71c71c720dd8792p-9L, -0x1.5f7268f0b5907438p-9L,
-0x1.e13cd410e0477de6p-13L, 0x1.9b0f31643442616ep-11L, 0x1.2527623a3472ae08p-14L,
-0x1.37f6bc8ef8b374dep-11L,-0x1.8c968886052b872ap-16L, 0x1.76baa9c6d3eeddbcp-11L
];
static immutable real[7] LargeStirlingCoeffs = [ 1.0L,
8.33333333333333333333E-2L, 3.47222222222222222222E-3L,
-2.68132716049382716049E-3L, -2.29472093621399176955E-4L,
7.84039221720066627474E-4L, 6.97281375836585777429E-5L
];
real w = 1.0L/x;
real y = exp(x);
if ( x > 1024.0L )
{
// For large x, use rational coefficients from the analytical expansion.
w = poly(w, LargeStirlingCoeffs);
// Avoid overflow in pow()
real v = pow( x, 0.5L * x - 0.25L );
y = v * (v / y);
}
else
{
w = 1.0L + w * poly( w, SmallStirlingCoeffs);
static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
{
// Avoid overflow in pow() for 64-bit reals
if (x > 143.0)
{
real v = pow( x, 0.5 * x - 0.25 );
y = v * (v / y);
}
else
{
y = pow( x, x - 0.5 ) / y;
}
}
else
{
y = pow( x, x - 0.5L ) / y;
}
}
y = SQRT2PI * y * w;
return y;
}
/*
* Helper function: Incomplete gamma function computed by Temme's expansion.
*
* This is a port of igamma_temme_large from Boost.
*
*/
real igammaTemmeLarge(real a, real x)
{
static immutable real[][13] coef = [
[ -0.333333333333333333333L, 0.0833333333333333333333L,
-0.0148148148148148148148L, 0.00115740740740740740741L,
0.000352733686067019400353L, -0.0001787551440329218107L,
0.39192631785224377817e-4L, -0.218544851067999216147e-5L,
-0.18540622107151599607e-5L, 0.829671134095308600502e-6L,
-0.176659527368260793044e-6L, 0.670785354340149858037e-8L,
0.102618097842403080426e-7L, -0.438203601845335318655e-8L,
0.914769958223679023418e-9L, -0.255141939949462497669e-10L,
-0.583077213255042506746e-10L, 0.243619480206674162437e-10L,
-0.502766928011417558909e-11L ],
[ -0.00185185185185185185185L, -0.00347222222222222222222L,
0.00264550264550264550265L, -0.000990226337448559670782L,
0.000205761316872427983539L, -0.40187757201646090535e-6L,
-0.18098550334489977837e-4L, 0.764916091608111008464e-5L,
-0.161209008945634460038e-5L, 0.464712780280743434226e-8L,
0.137863344691572095931e-6L, -0.575254560351770496402e-7L,
0.119516285997781473243e-7L, -0.175432417197476476238e-10L,
-0.100915437106004126275e-8L, 0.416279299184258263623e-9L,
-0.856390702649298063807e-10L ],
[ 0.00413359788359788359788L, -0.00268132716049382716049L,
0.000771604938271604938272L, 0.200938786008230452675e-5L,
-0.000107366532263651605215L, 0.529234488291201254164e-4L,
-0.127606351886187277134e-4L, 0.342357873409613807419e-7L,
0.137219573090629332056e-5L, -0.629899213838005502291e-6L,
0.142806142060642417916e-6L, -0.204770984219908660149e-9L,
-0.140925299108675210533e-7L, 0.622897408492202203356e-8L,
-0.136704883966171134993e-8L ],
[ 0.000649434156378600823045L, 0.000229472093621399176955L,
-0.000469189494395255712128L, 0.000267720632062838852962L,
-0.756180167188397641073e-4L, -0.239650511386729665193e-6L,
0.110826541153473023615e-4L, -0.56749528269915965675e-5L,
0.142309007324358839146e-5L, -0.278610802915281422406e-10L,
-0.169584040919302772899e-6L, 0.809946490538808236335e-7L,
-0.191111684859736540607e-7L ],
[ -0.000861888290916711698605L, 0.000784039221720066627474L,
-0.000299072480303190179733L, -0.146384525788434181781e-5L,
0.664149821546512218666e-4L, -0.396836504717943466443e-4L,
0.113757269706784190981e-4L, 0.250749722623753280165e-9L,
-0.169541495365583060147e-5L, 0.890750753220530968883e-6L,
-0.229293483400080487057e-6L ],
[ -0.000336798553366358150309L, -0.697281375836585777429e-4L,
0.000277275324495939207873L, -0.000199325705161888477003L,
0.679778047793720783882e-4L, 0.141906292064396701483e-6L,
-0.135940481897686932785e-4L, 0.801847025633420153972e-5L,
-0.229148117650809517038e-5L ],
[ 0.000531307936463992223166L, -0.000592166437353693882865L,
0.000270878209671804482771L, 0.790235323266032787212e-6L,
-0.815396936756196875093e-4L, 0.561168275310624965004e-4L,
-0.183291165828433755673e-4L, -0.307961345060330478256e-8L,
0.346515536880360908674e-5L, -0.20291327396058603727e-5L,
0.57887928631490037089e-6L ],
[ 0.000344367606892377671254L, 0.517179090826059219337e-4L,
-0.000334931610811422363117L, 0.000281269515476323702274L,
-0.000109765822446847310235L, -0.127410090954844853795e-6L,
0.277444515115636441571e-4L, -0.182634888057113326614e-4L,
0.578769494973505239894e-5L ],
[ -0.000652623918595309418922L, 0.000839498720672087279993L,
-0.000438297098541721005061L, -0.696909145842055197137e-6L,
0.000166448466420675478374L, -0.000127835176797692185853L,
0.462995326369130429061e-4L ],
[ -0.000596761290192746250124L, -0.720489541602001055909e-4L,
0.000678230883766732836162L, -0.0006401475260262758451L,
0.000277501076343287044992L ],
[ 0.00133244544948006563713L, -0.0019144384985654775265L,
0.00110893691345966373396L ],
[ 0.00157972766073083495909L, 0.000162516262783915816899L,
-0.00206334210355432762645L, 0.00213896861856890981541L,
-0.00101085593912630031708L ],
[ -0.00407251211951401664727L, 0.00640336283380806979482L,
-0.00404101610816766177474L ]
];
// avoid nans when one of the arguments is inf:
if (x == real.infinity && a != real.infinity)
return 0;
if (x != real.infinity && a == real.infinity)
return 1;
real sigma = (x - a) / a;
real phi = sigma - log(sigma + 1);
real y = a * phi;
real z = sqrt(2 * phi);
if (x < a)
z = -z;
real[13] workspace;
foreach (i; 0 .. coef.length)
workspace[i] = poly(z, coef[i]);
real result = poly(1 / a, workspace);
result *= exp(-y) / sqrt(2 * PI * a);
if (x < a)
result = -result;
result += erfc(sqrt(y)) / 2;
return result;
}
} // private
public:
/// The maximum value of x for which gamma(x) < real.infinity.
static if (floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
enum real MAXGAMMA = 1755.5483429L;
else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
enum real MAXGAMMA = 1755.5483429L;
else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended53)
enum real MAXGAMMA = 1755.5483429L;
else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
enum real MAXGAMMA = 171.6243769L;
else
static assert(0, "missing MAXGAMMA for other real types");
/*****************************************************
* The Gamma function, $(GAMMA)(x)
*
* $(GAMMA)(x) is a generalisation of the factorial function
* to real and complex numbers.
* Like x!, $(GAMMA)(x+1) = x*$(GAMMA)(x).
*
* Mathematically, if z.re > 0 then
* $(GAMMA)(z) = $(INTEGRATE 0, ∞) $(POWER t, z-1)$(POWER e, -t) dt
*
* $(TABLE_SV
* $(SVH x, $(GAMMA)(x) )
* $(SV $(NAN), $(NAN) )
* $(SV ±0.0, ±∞)
* $(SV integer > 0, (x-1)! )
* $(SV integer < 0, $(NAN) )
* $(SV +∞, +∞ )
* $(SV -∞, $(NAN) )
* )
*/
real gamma(real x)
{
/* Based on code from the CEPHES library.
* CEPHES code Copyright 1994 by Stephen L. Moshier
*
* Arguments |x| <= 13 are reduced by recurrence and the function
* approximated by a rational function of degree 7/8 in the
* interval (2,3). Large arguments are handled by Stirling's
* formula. Large negative arguments are made positive using
* a reflection formula.
*/
real q, z;
if (isNaN(x)) return x;
if (x == -x.infinity) return real.nan;
if ( fabs(x) > MAXGAMMA ) return real.infinity;
if (x == 0) return 1.0 / x; // +- infinity depending on sign of x, create an exception.
q = fabs(x);
if ( q > 13.0L )
{
// Large arguments are handled by Stirling's
// formula. Large negative arguments are made positive using
// the reflection formula.
if ( x < 0.0L )
{
if (x < -1/real.epsilon)
{
// Large negatives lose all precision
return real.nan;
}
int sgngam = 1; // sign of gamma.
long intpart = cast(long)(q);
if (q == intpart)
return real.nan; // poles for all integers <0.
real p = intpart;
if ( (intpart & 1) == 0 )
sgngam = -1;
z = q - p;
if ( z > 0.5L )
{
p += 1.0L;
z = q - p;
}
z = q * sin( PI * z );
z = fabs(z) * gammaStirling(q);
if ( z <= PI/real.max ) return sgngam * real.infinity;
return sgngam * PI/z;
}
else
{
return gammaStirling(x);
}
}
// Arguments |x| <= 13 are reduced by recurrence and the function
// approximated by a rational function of degree 7/8 in the
// interval (2,3).
z = 1.0L;
while ( x >= 3.0L )
{
x -= 1.0L;
z *= x;
}
while ( x < -0.03125L )
{
z /= x;
x += 1.0L;
}
if ( x <= 0.03125L )
{
if ( x == 0.0L )
return real.nan;
else
{
if ( x < 0.0L )
{
x = -x;
return z / (x * poly( x, GammaSmallNegCoeffs ));
}
else
{
return z / (x * poly( x, GammaSmallCoeffs ));
}
}
}
while ( x < 2.0L )
{
z /= x;
x += 1.0L;
}
if ( x == 2.0L ) return z;
x -= 2.0L;
return z * poly( x, GammaNumeratorCoeffs ) / poly( x, GammaDenominatorCoeffs );
}
@safe unittest
{
// gamma(n) = factorial(n-1) if n is an integer.
real fact = 1.0L;
for (int i=1; fact<real.max; ++i)
{
// Require exact equality for small factorials
if (i<14) assert(gamma(i*1.0L) == fact);
assert(feqrel(gamma(i*1.0L), fact) >= real.mant_dig-15);
fact *= (i*1.0L);
}
assert(gamma(0.0) == real.infinity);
assert(gamma(-0.0) == -real.infinity);
assert(isNaN(gamma(-1.0)));
assert(isNaN(gamma(-15.0)));
assert(isIdentical(gamma(NaN(0xABC)), NaN(0xABC)));
assert(gamma(real.infinity) == real.infinity);
assert(gamma(real.max) == real.infinity);
assert(isNaN(gamma(-real.infinity)));
assert(gamma(real.min_normal*real.epsilon) == real.infinity);
assert(gamma(MAXGAMMA)< real.infinity);
assert(gamma(MAXGAMMA*2) == real.infinity);
// Test some high-precision values (50 decimal digits)
real SQRT_PI = 1.77245385090551602729816748334114518279754945612238L;
assert(feqrel(gamma(0.5L), SQRT_PI) >= real.mant_dig-1);
assert(feqrel(gamma(17.25L), 4.224986665692703551570937158682064589938e13L) >= real.mant_dig-4);
assert(feqrel(gamma(1.0 / 3.0L), 2.67893853470774763365569294097467764412868937795730L) >= real.mant_dig-2);
assert(feqrel(gamma(0.25L),
3.62560990822190831193068515586767200299516768288006L) >= real.mant_dig-1);
assert(feqrel(gamma(1.0 / 5.0L),
4.59084371199880305320475827592915200343410999829340L) >= real.mant_dig-1);
}
/*****************************************************
* Natural logarithm of gamma function.
*
* Returns the base e (2.718...) logarithm of the absolute
* value of the gamma function of the argument.
*
* For reals, logGamma is equivalent to log(fabs(gamma(x))).
*
* $(TABLE_SV
* $(SVH x, logGamma(x) )
* $(SV $(NAN), $(NAN) )
* $(SV integer <= 0, +∞ )
* $(SV ±∞, +∞ )
* )
*/
real logGamma(real x)
{
/* Based on code from the CEPHES library.
* CEPHES code Copyright 1994 by Stephen L. Moshier
*
* For arguments greater than 33, the logarithm of the gamma
* function is approximated by the logarithmic version of
* Stirling's formula using a polynomial approximation of
* degree 4. Arguments between -33 and +33 are reduced by
* recurrence to the interval [2,3] of a rational approximation.
* The cosecant reflection formula is employed for arguments
* less than -33.
*/
real q, w, z, f, nx;
if (isNaN(x)) return x;
if (fabs(x) == x.infinity) return x.infinity;
if ( x < -34.0L )
{
q = -x;
w = logGamma(q);
real p = floor(q);
if ( p == q )
return real.infinity;
int intpart = cast(int)(p);
real sgngam = 1;
if ( (intpart & 1) == 0 )
sgngam = -1;
z = q - p;
if ( z > 0.5L )
{
p += 1.0L;
z = p - q;
}
z = q * sin( PI * z );
if ( z == 0.0L )
return sgngam * real.infinity;
/* z = LOGPI - logl( z ) - w; */
z = log( PI/z ) - w;
return z;
}
if ( x < 13.0L )
{
z = 1.0L;
nx = floor( x + 0.5L );
f = x - nx;
while ( x >= 3.0L )
{
nx -= 1.0L;
x = nx + f;
z *= x;
}
while ( x < 2.0L )
{
if ( fabs(x) <= 0.03125L )
{
if ( x == 0.0L )
return real.infinity;
if ( x < 0.0L )
{
x = -x;
q = z / (x * poly( x, GammaSmallNegCoeffs));
} else
q = z / (x * poly( x, GammaSmallCoeffs));
return log( fabs(q) );
}
z /= nx + f;
nx += 1.0L;
x = nx + f;
}
z = fabs(z);
if ( x == 2.0L )
return log(z);
x = (nx - 2.0L) + f;
real p = x * rationalPoly( x, logGammaNumerator, logGammaDenominator);
return log(z) + p;
}
// const real MAXLGM = 1.04848146839019521116e+4928L;
// if ( x > MAXLGM ) return sgngaml * real.infinity;
const real LOGSQRT2PI = 0.91893853320467274178L; // log( sqrt( 2*pi ) )
q = ( x - 0.5L ) * log(x) - x + LOGSQRT2PI;
if (x > 1.0e10L) return q;
real p = 1.0L / (x*x);
q += poly( p, logGammaStirlingCoeffs ) / x;
return q ;
}
@safe unittest
{
assert(isIdentical(logGamma(NaN(0xDEF)), NaN(0xDEF)));
assert(logGamma(real.infinity) == real.infinity);
assert(logGamma(-1.0) == real.infinity);
assert(logGamma(0.0) == real.infinity);
assert(logGamma(-50.0) == real.infinity);
assert(isIdentical(0.0L, logGamma(1.0L)));
assert(isIdentical(0.0L, logGamma(2.0L)));
assert(logGamma(real.min_normal*real.epsilon) == real.infinity);
assert(logGamma(-real.min_normal*real.epsilon) == real.infinity);
// x, correct loggamma(x), correct d/dx loggamma(x).
immutable static real[] testpoints = [
8.0L, 8.525146484375L + 1.48766904143001655310E-5, 2.01564147795560999654E0L,
8.99993896484375e-1L, 6.6375732421875e-2L + 5.11505711292524166220E-6L, -7.54938684259372234258E-1,
7.31597900390625e-1L, 2.2369384765625e-1 + 5.21506341809849792422E-6L, -1.13355566660398608343E0L,
2.31639862060546875e-1L, 1.3686676025390625L + 1.12609441752996145670E-5L, -4.56670961813812679012E0,
1.73162841796875L, -8.88214111328125e-2L + 3.36207740803753034508E-6L, 2.33339034686200586920E-1L,
1.23162841796875L, -9.3902587890625e-2L + 1.28765089229009648104E-5L, -2.49677345775751390414E-1L,
7.3786976294838206464e19L, 3.301798506038663053312e21L - 1.656137564136932662487046269677E5L,
4.57477139169563904215E1L,
1.08420217248550443401E-19L, 4.36682586669921875e1L + 1.37082843669932230418E-5L,
-9.22337203685477580858E18L,
1.0L, 0.0L, -5.77215664901532860607E-1L,
2.0L, 0.0L, 4.22784335098467139393E-1L,
-0.5L, 1.2655029296875L + 9.19379714539648894580E-6L, 3.64899739785765205590E-2L,
-1.5L, 8.6004638671875e-1L + 6.28657731014510932682E-7L, 7.03156640645243187226E-1L,
-2.5L, -5.6243896484375E-2L + 1.79986700949327405470E-7, 1.10315664064524318723E0L,
-3.5L, -1.30902099609375L + 1.43111007079536392848E-5L, 1.38887092635952890151E0L
];
// TODO: test derivatives as well.
for (int i=0; i<testpoints.length; i+=3)
{
assert( feqrel(logGamma(testpoints[i]), testpoints[i+1]) > real.mant_dig-5);
if (testpoints[i]<MAXGAMMA)
{
assert( feqrel(log(fabs(gamma(testpoints[i]))), testpoints[i+1]) > real.mant_dig-5);
}
}
assert(feqrel(logGamma(-50.2L),log(fabs(gamma(-50.2L)))) > real.mant_dig-2);
assert(feqrel(logGamma(-0.008L),log(fabs(gamma(-0.008L)))) > real.mant_dig-2);
assert(feqrel(logGamma(-38.8L),log(fabs(gamma(-38.8L)))) > real.mant_dig-4);
static if (real.mant_dig >= 64) // incl. 80-bit reals
assert(feqrel(logGamma(1500.0L),log(gamma(1500.0L))) > real.mant_dig-2);
else static if (real.mant_dig >= 53) // incl. 64-bit reals
assert(feqrel(logGamma(150.0L),log(gamma(150.0L))) > real.mant_dig-2);
}
private {
/*
* These value can be calculated like this:
* 1) Get exact real.max/min_normal/epsilon from compiler:
* writefln!"%a"(real.max/min_normal_epsilon)
* 2) Convert for Wolfram Alpha
* 0xf.fffffffffffffffp+16380 ==> (f.fffffffffffffff base 16) * 2^16380
* 3) Calculate result on wofram alpha:
* http://www.wolframalpha.com/input/?i=ln((1.ffffffffffffffffffffffffffff+base+16)+*+2%5E16383)+in+base+2
* 4) Convert to proper format:
* string mantissa = "1.011...";
* write(mantissa[0 .. 2]); mantissa = mantissa[2 .. $];
* for (size_t i = 0; i < mantissa.length/4; i++)
* {
* writef!"%x"(to!ubyte(mantissa[0 .. 4], 2)); mantissa = mantissa[4 .. $];
* }
*/
static if (floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
{
enum real MAXLOG = 0x1.62e42fefa39ef35793c7673007e6p+13L; // log(real.max)
enum real MINLOG = -0x1.6546282207802c89d24d65e96274p+13L; // log(real.min_normal*real.epsilon) = log(smallest denormal)
}
else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
{
enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max)
enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min_normal*real.epsilon) = log(smallest denormal)
}
else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended53)
{
enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max)
enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min_normal*real.epsilon) = log(smallest denormal)
}
else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
{
enum real MAXLOG = 0x1.62e42fefa39efp+9L; // log(real.max)
enum real MINLOG = -0x1.74385446d71c3p+9L; // log(real.min_normal*real.epsilon) = log(smallest denormal)
}
else
static assert(0, "missing MAXLOG and MINLOG for other real types");
enum real BETA_BIG = 9.223372036854775808e18L;
enum real BETA_BIGINV = 1.084202172485504434007e-19L;
}
/** Incomplete beta integral
*
* Returns incomplete beta integral of the arguments, evaluated
* from zero to x. The regularized incomplete beta function is defined as
*
* betaIncomplete(a, b, x) = Γ(a+b)/(Γ(a) Γ(b)) *
* $(INTEGRATE 0, x) $(POWER t, a-1)$(POWER (1-t),b-1) dt
*
* and is the same as the cumulative distribution function.
*
* The domain of definition is 0 <= x <= 1. In this
* implementation a and b are restricted to positive values.
* The integral from x to 1 may be obtained by the symmetry
* relation
*
* betaIncompleteCompl(a, b, x ) = betaIncomplete( b, a, 1-x )
*
* The integral is evaluated by a continued fraction expansion
* or, when b*x is small, by a power series.
*/
real betaIncomplete(real aa, real bb, real xx )
{
if ( !(aa>0 && bb>0) )
{
if ( isNaN(aa) ) return aa;
if ( isNaN(bb) ) return bb;
return real.nan; // domain error
}
if (!(xx>0 && xx<1.0))
{
if (isNaN(xx)) return xx;
if ( xx == 0.0L ) return 0.0;
if ( xx == 1.0L ) return 1.0;
return real.nan; // domain error
}
if ( (bb * xx) <= 1.0L && xx <= 0.95L)
{
return betaDistPowerSeries(aa, bb, xx);
}
real x;
real xc; // = 1 - x
real a, b;
int flag = 0;
/* Reverse a and b if x is greater than the mean. */
if ( xx > (aa/(aa+bb)) )
{
// here x > aa/(aa+bb) and (bb*x>1 or x>0.95)
flag = 1;
a = bb;
b = aa;
xc = xx;
x = 1.0L - xx;
}
else
{
a = aa;
b = bb;
xc = 1.0L - xx;
x = xx;
}
if ( flag == 1 && (b * x) <= 1.0L && x <= 0.95L)
{
// here xx > aa/(aa+bb) and ((bb*xx>1) or xx>0.95) and (aa*(1-xx)<=1) and xx > 0.05
return 1.0 - betaDistPowerSeries(a, b, x); // note loss of precision
}
real w;
// Choose expansion for optimal convergence
// One is for x * (a+b+2) < (a+1),
// the other is for x * (a+b+2) > (a+1).
real y = x * (a+b-2.0L) - (a-1.0L);
if ( y < 0.0L )
{
w = betaDistExpansion1( a, b, x );
}
else
{
w = betaDistExpansion2( a, b, x ) / xc;
}
/* Multiply w by the factor
a b
x (1-x) Gamma(a+b) / ( a Gamma(a) Gamma(b) ) . */
y = a * log(x);
real t = b * log(xc);
if ( (a+b) < MAXGAMMA && fabs(y) < MAXLOG && fabs(t) < MAXLOG )
{
t = pow(xc,b);
t *= pow(x,a);
t /= a;
t *= w;
t *= gamma(a+b) / (gamma(a) * gamma(b));
}
else
{
/* Resort to logarithms. */
y += t + logGamma(a+b) - logGamma(a) - logGamma(b);
y += log(w/a);
t = exp(y);
/+
// There seems to be a bug in Cephes at this point.
// Problems occur for y > MAXLOG, not y < MINLOG.
if ( y < MINLOG )
{
t = 0.0L;
}
else
{
t = exp(y);
}
+/
}
if ( flag == 1 )
{
/+ // CEPHES includes this code, but I think it is erroneous.
if ( t <= real.epsilon )
{
t = 1.0L - real.epsilon;
} else
+/
t = 1.0L - t;
}
return t;
}
/** Inverse of incomplete beta integral
*
* Given y, the function finds x such that
*
* betaIncomplete(a, b, x) == y
*
* Newton iterations or interval halving is used.
*/
real betaIncompleteInv(real aa, real bb, real yy0 )
{
real a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt;
int i, rflg, dir, nflg;
if (isNaN(yy0)) return yy0;
if (isNaN(aa)) return aa;
if (isNaN(bb)) return bb;
if ( yy0 <= 0.0L )
return 0.0L;
if ( yy0 >= 1.0L )
return 1.0L;
x0 = 0.0L;
yl = 0.0L;
x1 = 1.0L;
yh = 1.0L;
if ( aa <= 1.0L || bb <= 1.0L )
{
dithresh = 1.0e-7L;
rflg = 0;
a = aa;
b = bb;
y0 = yy0;
x = a/(a+b);
y = betaIncomplete( a, b, x );
nflg = 0;
goto ihalve;
}
else
{
nflg = 0;
dithresh = 1.0e-4L;
}
// approximation to inverse function
yp = -normalDistributionInvImpl( yy0 );
if ( yy0 > 0.5L )
{
rflg = 1;
a = bb;
b = aa;
y0 = 1.0L - yy0;
yp = -yp;
}
else
{
rflg = 0;
a = aa;
b = bb;
y0 = yy0;
}
lgm = (yp * yp - 3.0L)/6.0L;
x = 2.0L/( 1.0L/(2.0L * a-1.0L) + 1.0L/(2.0L * b - 1.0L) );
d = yp * sqrt( x + lgm ) / x
- ( 1.0L/(2.0L * b - 1.0L) - 1.0L/(2.0L * a - 1.0L) )
* (lgm + (5.0L/6.0L) - 2.0L/(3.0L * x));
d = 2.0L * d;
if ( d < MINLOG )
{
x = 1.0L;
goto under;
}
x = a/( a + b * exp(d) );
y = betaIncomplete( a, b, x );
yp = (y - y0)/y0;
if ( fabs(yp) < 0.2 )
goto newt;
/* Resort to interval halving if not close enough. */
ihalve:
dir = 0;
di = 0.5L;
for ( i=0; i<400; i++ )
{
if ( i != 0 )
{
x = x0 + di * (x1 - x0);
if ( x == 1.0L )
{
x = 1.0L - real.epsilon;
}
if ( x == 0.0L )
{
di = 0.5;
x = x0 + di * (x1 - x0);
if ( x == 0.0 )
goto under;
}
y = betaIncomplete( a, b, x );
yp = (x1 - x0)/(x1 + x0);
if ( fabs(yp) < dithresh )
goto newt;
yp = (y-y0)/y0;
if ( fabs(yp) < dithresh )
goto newt;
}
if ( y < y0 )
{
x0 = x;
yl = y;
if ( dir < 0 )
{
dir = 0;
di = 0.5L;
} else if ( dir > 3 )
di = 1.0L - (1.0L - di) * (1.0L - di);
else if ( dir > 1 )
di = 0.5L * di + 0.5L;
else
di = (y0 - y)/(yh - yl);
dir += 1;
if ( x0 > 0.95L )
{
if ( rflg == 1 )
{
rflg = 0;
a = aa;
b = bb;
y0 = yy0;
}
else
{
rflg = 1;
a = bb;
b = aa;
y0 = 1.0 - yy0;
}
x = 1.0L - x;
y = betaIncomplete( a, b, x );
x0 = 0.0;
yl = 0.0;
x1 = 1.0;
yh = 1.0;
goto ihalve;
}
}
else
{
x1 = x;
if ( rflg == 1 && x1 < real.epsilon )
{
x = 0.0L;
goto done;
}
yh = y;
if ( dir > 0 )
{
dir = 0;
di = 0.5L;
}
else if ( dir < -3 )
di = di * di;
else if ( dir < -1 )
di = 0.5L * di;
else
di = (y - y0)/(yh - yl);
dir -= 1;
}
}
if ( x0 >= 1.0L )
{
// partial loss of precision
x = 1.0L - real.epsilon;
goto done;
}
if ( x <= 0.0L )
{
under:
// underflow has occurred
x = real.min_normal * real.min_normal;
goto done;
}
newt:
if ( nflg )
{
goto done;
}
nflg = 1;
lgm = logGamma(a+b) - logGamma(a) - logGamma(b);
for ( i=0; i<15; i++ )
{
/* Compute the function at this point. */
if ( i != 0 )
y = betaIncomplete(a,b,x);
if ( y < yl )
{
x = x0;
y = yl;
}
else if ( y > yh )
{
x = x1;
y = yh;
}
else if ( y < y0 )
{
x0 = x;
yl = y;
}
else
{
x1 = x;
yh = y;
}
if ( x == 1.0L || x == 0.0L )
break;
/* Compute the derivative of the function at this point. */
d = (a - 1.0L) * log(x) + (b - 1.0L) * log(1.0L - x) + lgm;
if ( d < MINLOG )
{
goto done;
}
if ( d > MAXLOG )
{
break;
}
d = exp(d);
/* Compute the step to the next approximation of x. */
d = (y - y0)/d;
xt = x - d;
if ( xt <= x0 )
{
y = (x - x0) / (x1 - x0);
xt = x0 + 0.5L * y * (x - x0);
if ( xt <= 0.0L )
break;
}
if ( xt >= x1 )
{
y = (x1 - x) / (x1 - x0);
xt = x1 - 0.5L * y * (x1 - x);
if ( xt >= 1.0L )
break;
}
x = xt;
if ( fabs(d/x) < (128.0L * real.epsilon) )
goto done;
}
/* Did not converge. */
dithresh = 256.0L * real.epsilon;
goto ihalve;
done:
if ( rflg )
{
if ( x <= real.epsilon )
x = 1.0L - real.epsilon;
else
x = 1.0L - x;
}
return x;
}
@safe unittest { // also tested by the normal distribution
// check NaN propagation
assert(isIdentical(betaIncomplete(NaN(0xABC),2,3), NaN(0xABC)));
assert(isIdentical(betaIncomplete(7,NaN(0xABC),3), NaN(0xABC)));
assert(isIdentical(betaIncomplete(7,15,NaN(0xABC)), NaN(0xABC)));
assert(isIdentical(betaIncompleteInv(NaN(0xABC),1,17), NaN(0xABC)));
assert(isIdentical(betaIncompleteInv(2,NaN(0xABC),8), NaN(0xABC)));
assert(isIdentical(betaIncompleteInv(2,3, NaN(0xABC)), NaN(0xABC)));
assert(isNaN(betaIncomplete(-1, 2, 3)));
assert(betaIncomplete(1, 2, 0)==0);
assert(betaIncomplete(1, 2, 1)==1);
assert(isNaN(betaIncomplete(1, 2, 3)));
assert(betaIncompleteInv(1, 1, 0)==0);
assert(betaIncompleteInv(1, 1, 1)==1);
// Test against Mathematica betaRegularized[z,a,b]
// These arbitrary points are chosen to give good code coverage.
assert(feqrel(betaIncomplete(8, 10, 0.2L), 0.010_934_315_234_099_2L) >= real.mant_dig - 5);
assert(feqrel(betaIncomplete(2, 2.5L, 0.9L), 0.989_722_597_604_452_767_171_003_59L) >= real.mant_dig - 1);
static if (real.mant_dig >= 64) // incl. 80-bit reals
assert(feqrel(betaIncomplete(1000, 800, 0.5L), 1.179140859734704555102808541457164E-06L) >= real.mant_dig - 13);
else
assert(feqrel(betaIncomplete(1000, 800, 0.5L), 1.179140859734704555102808541457164E-06L) >= real.mant_dig - 14);
assert(feqrel(betaIncomplete(0.0001, 10000, 0.0001L), 0.999978059362107134278786L) >= real.mant_dig - 18);
assert(betaIncomplete(0.01L, 327726.7L, 0.545113L) == 1.0);
assert(feqrel(betaIncompleteInv(8, 10, 0.010_934_315_234_099_2L), 0.2L) >= real.mant_dig - 2);
assert(feqrel(betaIncomplete(0.01L, 498.437L, 0.0121433L), 0.99999664562033077636065L) >= real.mant_dig - 1);
assert(feqrel(betaIncompleteInv(5, 10, 0.2000002972865658842L), 0.229121208190918L) >= real.mant_dig - 3);
assert(feqrel(betaIncompleteInv(4, 7, 0.8000002209179505L), 0.483657360076904L) >= real.mant_dig - 3);
// Coverage tests. I don't have correct values for these tests, but
// these values cover most of the code, so they are useful for
// regression testing.
// Extensive testing failed to increase the coverage. It seems likely that about
// half the code in this function is unnecessary; there is potential for
// significant improvement over the original CEPHES code.
static if (real.mant_dig == 64) // 80-bit reals
{
assert(betaIncompleteInv(0.01L, 8e-48L, 5.45464e-20L) == 1-real.epsilon);
assert(betaIncompleteInv(0.01L, 8e-48L, 9e-26L) == 1-real.epsilon);
// Beware: a one-bit change in pow() changes almost all digits in the result!
assert(feqrel(
betaIncompleteInv(0x1.b3d151fbba0eb18p+1L, 1.2265e-19L, 2.44859e-18L),
0x1.c0110c8531d0952cp-1L
) > 10);
// This next case uncovered a one-bit difference in the FYL2X instruction
// between Intel and AMD processors. This difference gets magnified by 2^^38.
// WolframAlpha crashes attempting to calculate this.
assert(feqrel(betaIncompleteInv(0x1.ff1275ae5b939bcap-41L, 4.6713e18L, 0.0813601L),
0x1.f97749d90c7adba8p-63L) >= real.mant_dig - 39);
real a1 = 3.40483L;
assert(betaIncompleteInv(a1, 4.0640301659679627772e19L, 0.545113L) == 0x1.ba8c08108aaf5d14p-109L);
real b1 = 2.82847e-25L;
assert(feqrel(betaIncompleteInv(0.01L, b1, 9e-26L), 0x1.549696104490aa9p-830L) >= real.mant_dig-10);
// --- Problematic cases ---
// This is a situation where the series expansion fails to converge
assert( isNaN(betaIncompleteInv(0.12167L, 4.0640301659679627772e19L, 0.0813601L)));
// This next result is almost certainly erroneous.
// Mathematica states: "(cannot be determined by current methods)"
assert(betaIncomplete(1.16251e20L, 2.18e39L, 5.45e-20L) == -real.infinity);
// WolframAlpha gives no result for this, though indicates that it approximately 1.0 - 1.3e-9
assert(1 - betaIncomplete(0.01L, 328222, 4.0375e-5L) == 0x1.5f62926b4p-30L);
}
}
private {
// Implementation functions
// Continued fraction expansion #1 for incomplete beta integral
// Use when x < (a+1)/(a+b+2)
real betaDistExpansion1(real a, real b, real x )
{
real xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
real k1, k2, k3, k4, k5, k6, k7, k8;
real r, t, ans;
int n;
k1 = a;
k2 = a + b;
k3 = a;
k4 = a + 1.0L;
k5 = 1.0L;
k6 = b - 1.0L;
k7 = k4;
k8 = a + 2.0L;
pkm2 = 0.0L;
qkm2 = 1.0L;
pkm1 = 1.0L;
qkm1 = 1.0L;
ans = 1.0L;
r = 1.0L;
n = 0;
const real thresh = 3.0L * real.epsilon;
do
{
xk = -( x * k1 * k2 )/( k3 * k4 );
pk = pkm1 + pkm2 * xk;
qk = qkm1 + qkm2 * xk;
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
xk = ( x * k5 * k6 )/( k7 * k8 );
pk = pkm1 + pkm2 * xk;
qk = qkm1 + qkm2 * xk;
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
if ( qk != 0.0L )
r = pk/qk;
if ( r != 0.0L )
{
t = fabs( (ans - r)/r );
ans = r;
}
else
{
t = 1.0L;
}
if ( t < thresh )
return ans;
k1 += 1.0L;
k2 += 1.0L;
k3 += 2.0L;
k4 += 2.0L;
k5 += 1.0L;
k6 -= 1.0L;
k7 += 2.0L;
k8 += 2.0L;
if ( (fabs(qk) + fabs(pk)) > BETA_BIG )
{
pkm2 *= BETA_BIGINV;
pkm1 *= BETA_BIGINV;
qkm2 *= BETA_BIGINV;
qkm1 *= BETA_BIGINV;
}
if ( (fabs(qk) < BETA_BIGINV) || (fabs(pk) < BETA_BIGINV) )
{
pkm2 *= BETA_BIG;
pkm1 *= BETA_BIG;
qkm2 *= BETA_BIG;
qkm1 *= BETA_BIG;
}
}
while ( ++n < 400 );
// loss of precision has occurred
// mtherr( "incbetl", PLOSS );
return ans;
}
// Continued fraction expansion #2 for incomplete beta integral
// Use when x > (a+1)/(a+b+2)
real betaDistExpansion2(real a, real b, real x )
{
real xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
real k1, k2, k3, k4, k5, k6, k7, k8;
real r, t, ans, z;
k1 = a;
k2 = b - 1.0L;
k3 = a;
k4 = a + 1.0L;
k5 = 1.0L;
k6 = a + b;
k7 = a + 1.0L;
k8 = a + 2.0L;
pkm2 = 0.0L;
qkm2 = 1.0L;
pkm1 = 1.0L;
qkm1 = 1.0L;
z = x / (1.0L-x);
ans = 1.0L;
r = 1.0L;
int n = 0;
const real thresh = 3.0L * real.epsilon;
do
{
xk = -( z * k1 * k2 )/( k3 * k4 );
pk = pkm1 + pkm2 * xk;
qk = qkm1 + qkm2 * xk;
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
xk = ( z * k5 * k6 )/( k7 * k8 );
pk = pkm1 + pkm2 * xk;
qk = qkm1 + qkm2 * xk;
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
if ( qk != 0.0L )
r = pk/qk;
if ( r != 0.0L )
{
t = fabs( (ans - r)/r );
ans = r;
} else
t = 1.0L;
if ( t < thresh )
return ans;
k1 += 1.0L;
k2 -= 1.0L;
k3 += 2.0L;
k4 += 2.0L;
k5 += 1.0L;
k6 += 1.0L;
k7 += 2.0L;
k8 += 2.0L;
if ( (fabs(qk) + fabs(pk)) > BETA_BIG )
{
pkm2 *= BETA_BIGINV;
pkm1 *= BETA_BIGINV;
qkm2 *= BETA_BIGINV;
qkm1 *= BETA_BIGINV;
}
if ( (fabs(qk) < BETA_BIGINV) || (fabs(pk) < BETA_BIGINV) )
{
pkm2 *= BETA_BIG;
pkm1 *= BETA_BIG;
qkm2 *= BETA_BIG;
qkm1 *= BETA_BIG;
}
} while ( ++n < 400 );
// loss of precision has occurred
//mtherr( "incbetl", PLOSS );
return ans;
}
/* Power series for incomplete gamma integral.
Use when b*x is small. */
real betaDistPowerSeries(real a, real b, real x )
{
real ai = 1.0L / a;
real u = (1.0L - b) * x;
real v = u / (a + 1.0L);
real t1 = v;
real t = u;
real n = 2.0L;
real s = 0.0L;
real z = real.epsilon * ai;
while ( fabs(v) > z )
{
u = (n - b) * x / n;
t *= u;
v = t / (a + n);
s += v;
n += 1.0L;
}
s += t1;
s += ai;
u = a * log(x);
if ( (a+b) < MAXGAMMA && fabs(u) < MAXLOG )
{
t = gamma(a+b)/(gamma(a)*gamma(b));
s = s * t * pow(x,a);
}
else
{
t = logGamma(a+b) - logGamma(a) - logGamma(b) + u + log(s);
if ( t < MINLOG )
{
s = 0.0L;
} else
s = exp(t);
}
return s;
}
}
/***************************************
* Incomplete gamma integral and its complement
*
* These functions are defined by
*
* gammaIncomplete = ( $(INTEGRATE 0, x) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
*
* gammaIncompleteCompl(a,x) = 1 - gammaIncomplete(a,x)
* = ($(INTEGRATE x, ∞) $(POWER e, -t) $(POWER t, a-1) dt )/ $(GAMMA)(a)
*
* In this implementation both arguments must be positive.
* The integral is evaluated by either a power series or
* continued fraction expansion, depending on the relative
* values of a and x.
*/
real gammaIncomplete(real a, real x )
in
{
assert(x >= 0);
assert(a > 0);
}
do
{
/* left tail of incomplete gamma function:
*
* inf. k
* a -x - x
* x e > ----------
* - -
* k=0 | (a+k+1)
*
*/
if (x == 0)
return 0.0L;
if ( (x > 1.0L) && (x > a ) )
return 1.0L - gammaIncompleteCompl(a,x);
real ax = a * log(x) - x - logGamma(a);
/+
if ( ax < MINLOGL ) return 0; // underflow
// { mtherr( "igaml", UNDERFLOW ); return( 0.0L ); }
+/
ax = exp(ax);
/* power series */
real r = a;
real c = 1.0L;
real ans = 1.0L;
do
{
r += 1.0L;
c *= x/r;
ans += c;
} while ( c/ans > real.epsilon );
return ans * ax/a;
}
/** ditto */
real gammaIncompleteCompl(real a, real x )
in
{
assert(x >= 0);
assert(a > 0);
}
do
{
if (x == 0)
return 1.0L;
if ( (x < 1.0L) || (x < a) )
return 1.0L - gammaIncomplete(a,x);
// DAC (Cephes bug fix): This is necessary to avoid
// spurious nans, eg
// log(x)-x = NaN when x = real.infinity
const real MAXLOGL = 1.1356523406294143949492E4L;
if (x > MAXLOGL)
return igammaTemmeLarge(a, x);
real ax = a * log(x) - x - logGamma(a);
//const real MINLOGL = -1.1355137111933024058873E4L;
// if ( ax < MINLOGL ) return 0; // underflow;
ax = exp(ax);
/* continued fraction */
real y = 1.0L - a;
real z = x + y + 1.0L;
real c = 0.0L;
real pk, qk, t;
real pkm2 = 1.0L;
real qkm2 = x;
real pkm1 = x + 1.0L;
real qkm1 = z * x;
real ans = pkm1/qkm1;
do
{
c += 1.0L;
y += 1.0L;
z += 2.0L;
real yc = y * c;
pk = pkm1 * z - pkm2 * yc;
qk = qkm1 * z - qkm2 * yc;
if ( qk != 0.0L )
{
real r = pk/qk;
t = fabs( (ans - r)/r );
ans = r;
}
else
{
t = 1.0L;
}
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
const real BIG = 9.223372036854775808e18L;
if ( fabs(pk) > BIG )
{
pkm2 /= BIG;
pkm1 /= BIG;
qkm2 /= BIG;
qkm1 /= BIG;
}
} while ( t > real.epsilon );
return ans * ax;
}
/** Inverse of complemented incomplete gamma integral
*
* Given a and p, the function finds x such that
*
* gammaIncompleteCompl( a, x ) = p.
*
* Starting with the approximate value x = a $(POWER t, 3), where
* t = 1 - d - normalDistributionInv(p) sqrt(d),
* and d = 1/9a,
* the routine performs up to 10 Newton iterations to find the
* root of incompleteGammaCompl(a,x) - p = 0.
*/
real gammaIncompleteComplInv(real a, real p)
in
{
assert(p >= 0 && p <= 1);
assert(a>0);
}
do
{
if (p == 0) return real.infinity;
real y0 = p;
const real MAXLOGL = 1.1356523406294143949492E4L;
real x0, x1, x, yl, yh, y, d, lgm, dithresh;
int i, dir;
/* bound the solution */
x0 = real.max;
yl = 0.0L;
x1 = 0.0L;
yh = 1.0L;
dithresh = 4.0 * real.epsilon;
/* approximation to inverse function */
d = 1.0L/(9.0L*a);
y = 1.0L - d - normalDistributionInvImpl(y0) * sqrt(d);
x = a * y * y * y;
lgm = logGamma(a);
for ( i=0; i<10; i++ )
{
if ( x > x0 || x < x1 )
goto ihalve;
y = gammaIncompleteCompl(a,x);
if ( y < yl || y > yh )
goto ihalve;
if ( y < y0 )
{
x0 = x;
yl = y;
}
else
{
x1 = x;
yh = y;
}
/* compute the derivative of the function at this point */
d = (a - 1.0L) * log(x0) - x0 - lgm;
if ( d < -MAXLOGL )
goto ihalve;
d = -exp(d);
/* compute the step to the next approximation of x */
d = (y - y0)/d;
x = x - d;
if ( i < 3 ) continue;
if ( fabs(d/x) < dithresh ) return x;
}
/* Resort to interval halving if Newton iteration did not converge. */
ihalve:
d = 0.0625L;
if ( x0 == real.max )
{
if ( x <= 0.0L )
x = 1.0L;
while ( x0 == real.max )
{
x = (1.0L + d) * x;
y = gammaIncompleteCompl( a, x );
if ( y < y0 )
{
x0 = x;
yl = y;
break;
}
d = d + d;
}
}
d = 0.5L;
dir = 0;
for ( i=0; i<400; i++ )
{
x = x1 + d * (x0 - x1);
y = gammaIncompleteCompl( a, x );
lgm = (x0 - x1)/(x1 + x0);
if ( fabs(lgm) < dithresh )
break;
lgm = (y - y0)/y0;
if ( fabs(lgm) < dithresh )
break;
if ( x <= 0.0L )
break;
if ( y > y0 )
{
x1 = x;
yh = y;
if ( dir < 0 )
{
dir = 0;
d = 0.5L;
} else if ( dir > 1 )
d = 0.5L * d + 0.5L;
else
d = (y0 - yl)/(yh - yl);
dir += 1;
}
else
{
x0 = x;
yl = y;
if ( dir > 0 )
{
dir = 0;
d = 0.5L;
} else if ( dir < -1 )
d = 0.5L * d;
else
d = (y0 - yl)/(yh - yl);
dir -= 1;
}
}
/+
if ( x == 0.0L )
mtherr( "igamil", UNDERFLOW );
+/
return x;
}
@safe unittest
{
//Values from Excel's GammaInv(1-p, x, 1)
assert(fabs(gammaIncompleteComplInv(1, 0.5L) - 0.693147188044814L) < 0.00000005L);
assert(fabs(gammaIncompleteComplInv(12, 0.99L) - 5.42818075054289L) < 0.00000005L);
assert(fabs(gammaIncompleteComplInv(100, 0.8L) - 91.5013985848288L) < 0.000005L);
assert(gammaIncomplete(1, 0)==0);
assert(gammaIncompleteCompl(1, 0)==1);
assert(gammaIncomplete(4545, real.infinity)==1);
// Values from Excel's (1-GammaDist(x, alpha, 1, TRUE))
assert(fabs(1.0L-gammaIncompleteCompl(0.5L, 2) - 0.954499729507309L) < 0.00000005L);
assert(fabs(gammaIncomplete(0.5L, 2) - 0.954499729507309L) < 0.00000005L);
// Fixed Cephes bug:
assert(gammaIncompleteCompl(384, real.infinity)==0);
assert(gammaIncompleteComplInv(3, 0)==real.infinity);
// Fixed a bug that caused gammaIncompleteCompl to return a wrong value when
// x was larger than a, but not by much, and both were large:
// The value is from WolframAlpha (Gamma[100000, 100001, inf] / Gamma[100000])
static if (real.mant_dig >= 64) // incl. 80-bit reals
assert(fabs(gammaIncompleteCompl(100000, 100001) - 0.49831792109L) < 0.000000000005L);
else
assert(fabs(gammaIncompleteCompl(100000, 100001) - 0.49831792109L) < 0.00000005L);
}
// DAC: These values are Bn / n for n=2,4,6,8,10,12,14.
immutable real [7] Bn_n = [
1.0L/(6*2), -1.0L/(30*4), 1.0L/(42*6), -1.0L/(30*8),
5.0L/(66*10), -691.0L/(2730*12), 7.0L/(6*14) ];
/** Digamma function
*
* The digamma function is the logarithmic derivative of the gamma function.
*
* digamma(x) = d/dx logGamma(x)
*
* References:
* 1. Abramowitz, M., and Stegun, I. A. (1970).
* Handbook of mathematical functions. Dover, New York,
* pages 258-259, equations 6.3.6 and 6.3.18.
*/
real digamma(real x)
{
// Based on CEPHES, Stephen L. Moshier.
real p, q, nz, s, w, y, z;
long i, n;
int negative;
negative = 0;
nz = 0.0;
if ( x <= 0.0 )
{
negative = 1;
q = x;
p = floor(q);
if ( p == q )
{
return real.nan; // singularity.
}
/* Remove the zeros of tan(PI x)
* by subtracting the nearest integer from x
*/
nz = q - p;
if ( nz != 0.5 )
{
if ( nz > 0.5 )
{
p += 1.0;
nz = q - p;
}
nz = PI/tan(PI*nz);
}
else
{
nz = 0.0;
}
x = 1.0 - x;
}
// check for small positive integer
if ((x <= 13.0) && (x == floor(x)) )
{
y = 0.0;
n = lrint(x);
// DAC: CEPHES bugfix. Cephes did this in reverse order, which
// created a larger roundoff error.
for (i=n-1; i>0; --i)
{
y+=1.0L/i;
}
y -= EULERGAMMA;
goto done;
}
s = x;
w = 0.0;
while ( s < 10.0 )
{
w += 1.0/s;
s += 1.0;
}
if ( s < 1.0e17L )
{
z = 1.0/(s * s);
y = z * poly(z, Bn_n);
} else
y = 0.0;
y = log(s) - 0.5L/s - y - w;
done:
if ( negative )
{
y -= nz;
}
return y;
}
@safe unittest
{
// Exact values
assert(digamma(1.0)== -EULERGAMMA);
assert(feqrel(digamma(0.25), -PI/2 - 3* LN2 - EULERGAMMA) >= real.mant_dig-7);
assert(feqrel(digamma(1.0L/6), -PI/2 *sqrt(3.0L) - 2* LN2 -1.5*log(3.0L) - EULERGAMMA) >= real.mant_dig-7);
assert(digamma(-5.0).isNaN());
assert(feqrel(digamma(2.5), -EULERGAMMA - 2*LN2 + 2.0 + 2.0L/3) >= real.mant_dig-9);
assert(isIdentical(digamma(NaN(0xABC)), NaN(0xABC)));
for (int k=1; k<40; ++k)
{
real y=0;
for (int u=k; u >= 1; --u)
{
y += 1.0L/u;
}
assert(feqrel(digamma(k+1.0), -EULERGAMMA + y) >= real.mant_dig-2);
}
}
/** Log Minus Digamma function
*
* logmdigamma(x) = log(x) - digamma(x)
*
* References:
* 1. Abramowitz, M., and Stegun, I. A. (1970).
* Handbook of mathematical functions. Dover, New York,
* pages 258-259, equations 6.3.6 and 6.3.18.
*/
real logmdigamma(real x)
{
if (x <= 0.0)
{
if (x == 0.0)
{
return real.infinity;
}
return real.nan;
}
real s = x;
real w = 0.0;
while ( s < 10.0 )
{
w += 1.0/s;
s += 1.0;
}
real y;
if ( s < 1.0e17L )
{
immutable real z = 1.0/(s * s);
y = z * poly(z, Bn_n);
} else
y = 0.0;
return x == s ? y + 0.5L/s : (log(x/s) + 0.5L/s + y + w);
}
@safe unittest
{
assert(logmdigamma(-5.0).isNaN());
assert(isIdentical(logmdigamma(NaN(0xABC)), NaN(0xABC)));
assert(logmdigamma(0.0) == real.infinity);
for (auto x = 0.01; x < 1.0; x += 0.1)
assert(isClose(digamma(x), log(x) - logmdigamma(x)));
for (auto x = 1.0; x < 15.0; x += 1.0)
assert(isClose(digamma(x), log(x) - logmdigamma(x)));
}
/** Inverse of the Log Minus Digamma function
*
* Returns x such $(D log(x) - digamma(x) == y).
*
* References:
* 1. Abramowitz, M., and Stegun, I. A. (1970).
* Handbook of mathematical functions. Dover, New York,
* pages 258-259, equation 6.3.18.
*
* Authors: Ilya Yaroshenko
*/
real logmdigammaInverse(real y)
{
import std.numeric : findRoot;
// FIXME: should be returned back to enum.
// Fix requires CTFEable `log` on non-x86 targets (check both LDC and GDC).
immutable maxY = logmdigamma(real.min_normal);
assert(maxY > 0 && maxY <= real.max);
if (y >= maxY)
{
//lim x->0 (log(x)-digamma(x))*x == 1
return 1 / y;
}
if (y < 0)
{
return real.nan;
}
if (y < real.min_normal)
{
//6.3.18
return 0.5 / y;
}
if (y > 0)
{
// x/2 <= logmdigamma(1 / x) <= x, x > 0
// calls logmdigamma ~6 times
return 1 / findRoot((real x) => logmdigamma(1 / x) - y, y, 2*y);
}
return y; //NaN
}
@safe unittest
{
import std.typecons;
//WolframAlpha, 22.02.2015
immutable Tuple!(real, real)[5] testData = [
tuple(1.0L, 0.615556766479594378978099158335549201923L),
tuple(1.0L/8, 4.15937801516894947161054974029150730555L),
tuple(1.0L/1024, 512.166612384991507850643277924243523243L),
tuple(0.000500083333325000003968249801594877323784632117L, 1000.0L),
tuple(1017.644138623741168814449776695062817947092468536L, 1.0L/1024),
];
foreach (test; testData)
assert(isClose(logmdigammaInverse(test[0]), test[1], 2e-15L));
assert(isClose(logmdigamma(logmdigammaInverse(1)), 1, 1e-15L));
assert(isClose(logmdigamma(logmdigammaInverse(real.min_normal)), real.min_normal, 1e-15L));
assert(isClose(logmdigamma(logmdigammaInverse(real.max/2)), real.max/2, 1e-15L));
assert(isClose(logmdigammaInverse(logmdigamma(1)), 1, 1e-15L));
assert(isClose(logmdigammaInverse(logmdigamma(real.min_normal)), real.min_normal, 1e-15L));
assert(isClose(logmdigammaInverse(logmdigamma(real.max/2)), real.max/2, 1e-15L));
}