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