|
132
|
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 = Γ
|
|
|
17 * INTEGRATE = $(BIG ∫<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 π) 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, ∞) $(POWER t, z-1)$(POWER e, -t) dt
|
|
|
171 *
|
|
|
172 * $(TABLE_SV
|
|
|
173 * $(SVH x, $(GAMMA)(x) )
|
|
|
174 * $(SV $(NAN), $(NAN) )
|
|
|
175 * $(SV ±0.0, ±∞)
|
|
|
176 * $(SV integer > 0, (x-1)! )
|
|
|
177 * $(SV integer < 0, $(NAN) )
|
|
|
178 * $(SV +∞, +∞ )
|
|
|
179 * $(SV -∞, $(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, +∞ )
|
|
|
314 * $(SV ±∞, +∞ )
|
|
|
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) = (Γ(x) Γ(y))/Γ(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) = Γ(a+b)/(Γ(a) Γ(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, ∞) $(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
|
