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