Mercurial > projects > ldc

comparison tango/tango/math/GammaFunction.d @ 132:1700239cab2e trunk

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