/std/internal/math/gammafunction.d

http://github.com/jcd/phobos · D · 1518 lines · 1096 code · 144 blank · 278 comment · 252 complexity · a20985dd8c78a99da3c02d40dcce0143 MD5 · raw file

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