View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  //  Copyright John Maddock 2006-7, 2013-20.
19  //  Copyright Paul A. Bristow 2007, 2013-14.
20  //  Copyright Nikhar Agrawal 2013-14
21  //  Copyright Christopher Kormanyos 2013-14, 2020
22  //  Use, modification and distribution are subject to the
23  //  Boost Software License, Version 1.0. (See accompanying file
24  //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
25  
26  package org.apache.commons.numbers.gamma;
27  
28  import java.util.function.DoubleSupplier;
29  import java.util.function.Supplier;
30  import org.apache.commons.numbers.core.Sum;
31  import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction;
32  import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction.Coefficient;
33  
34  /**
35   * Implementation of the
36   * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">Regularized Gamma functions</a> and
37   * <a href="https://mathworld.wolfram.com/IncompleteGammaFunction.html">Incomplete Gamma functions</a>.
38   *
39   * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
40   * {@code c++} implementation {@code <boost/math/special_functions/gamma.hpp>}.
41   * All work is copyright to the original authors and subject to the Boost Software License.
42   *
43   * @see
44   * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma.html">
45   * Boost C++ Gamma functions</a>
46   */
47  final class BoostGamma {
48      //
49      // Code ported from Boost 1.77.0
50      //
51      // boost/math/special_functions/gamma.hpp
52      // boost/math/special_functions/detail/igamma_large.hpp
53      // boost/math/special_functions/lanczos.hpp
54      //
55      // Original code comments are preserved.
56      //
57      // Changes to the Boost implementation:
58      // - Update method names to replace underscores with camel case
59      // - Explicitly inline the polynomial function evaluation
60      //   using Horner's method (https://en.wikipedia.org/wiki/Horner%27s_method)
61      // - Remove checks for under/overflow. In this implementation no error is raised
62      //   for overflow (infinity is returned) or underflow (sub-normal or zero is returned).
63      //   This follows the conventions in java.lang.Math for the same conditions.
64      // - Removed the pointer p_derivative in the gammaIncompleteImp. This is used
65      //   in the Boost code for the gamma_(p|q)_inv functions for a derivative
66      //   based inverse function. This is currently not supported.
67      // - Added extended precision arithmetic for some series summations or other computations.
68      //   The Boost default policy is to evaluate in long double for a double result. Extended
69      //   precision is not possible for the entire computation but has been used where
70      //   possible for some terms to reduce errors. Error reduction verified on the test data.
71      // - Altered the tgamma(x) function to use the double-precision NSWC Library of Mathematics
72      //   Subroutines when the error is lower. This is for the non Lanczos code.
73      // - Altered the condition used for the asymptotic approximation method to avoid
74      //   loss of precision in the series summation when a ~ z.
75      // - Altered series generators to use integer counters added to the double term
76      //   replacing directly incrementing a double term. When the term is large it cannot
77      //   be incremented: 1e16 + 1 == 1e16.
78      // - Removed unreachable code branch in tgammaDeltaRatioImpLanczos when z + delta == z.
79      //
80      // Note:
81      // The major source of error is in the function regularisedGammaPrefix when computing
82      // (z^a)(e^-z)/tgamma(a) with extreme input to the power and exponential terms.
83      // An extended precision pow and exp function returning a quad length result would
84      // be required to reduce error for these arguments. Tests using the Dfp class
85      // from o.a.c.math4.legacy.core have been demonstrated to effectively eliminate the
86      // errors from the power terms and improve accuracy on the current test data.
87      // In the interest of performance the Dfp class is not used in this version.
88  
89      /** Default epsilon value for relative error.
90       * This is equal to the Boost constant {@code boost::math::tools::epsilon<double>()}. */
91      private static final double EPSILON = 0x1.0p-52;
92      /** Value for the sqrt of the epsilon for relative error.
93       * This is equal to the Boost constant {@code boost::math::tools::root_epsilon<double>()}. */
94      private static final double ROOT_EPSILON = 1.4901161193847656E-8;
95      /** Approximate value for ln(Double.MAX_VALUE).
96       * This is equal to the Boost constant {@code boost::math::tools::log_max_value<double>()}.
97       * No term {@code x} should be used in {@code exp(x)} if {@code x > LOG_MAX_VALUE} to avoid
98       * overflow. */
99      private static final int LOG_MAX_VALUE = 709;
100     /** Approximate value for ln(Double.MIN_VALUE).
101      * This is equal to the Boost constant {@code boost::math::tools::log_min_value<double>()}.
102      * No term {@code x} should be used in {@code exp(x)} if {@code x < LOG_MIN_VALUE} to avoid
103      * underflow to sub-normal or zero. */
104     private static final int LOG_MIN_VALUE = -708;
105     /** The largest factorial that can be represented as a double.
106      * This is equal to the Boost constant {@code boost::math::max_factorial<double>::value}. */
107     private static final int MAX_FACTORIAL = 170;
108     /** The largest integer value for gamma(z) that can be represented as a double. */
109     private static final int MAX_GAMMA_Z = MAX_FACTORIAL + 1;
110     /** ln(sqrt(2 pi)). Computed to 25-digits precision. */
111     private static final double LOG_ROOT_TWO_PI = 0.9189385332046727417803297;
112     /** ln(pi). Computed to 25-digits precision. */
113     private static final double LOG_PI = 1.144729885849400174143427;
114     /** Euler's constant. */
115     private static final double EULER = 0.5772156649015328606065120900824024310;
116     /** The threshold value for choosing the Lanczos approximation. */
117     private static final int LANCZOS_THRESHOLD = 20;
118     /** 2^53. */
119     private static final double TWO_POW_53 = 0x1.0p53;
120 
121     /** All factorials that can be represented as a double. Size = 171. */
122     private static final double[] FACTORIAL = {
123         1,
124         1,
125         2,
126         6,
127         24,
128         120,
129         720,
130         5040,
131         40320,
132         362880.0,
133         3628800.0,
134         39916800.0,
135         479001600.0,
136         6227020800.0,
137         87178291200.0,
138         1307674368000.0,
139         20922789888000.0,
140         355687428096000.0,
141         6402373705728000.0,
142         121645100408832000.0,
143         0.243290200817664e19,
144         0.5109094217170944e20,
145         0.112400072777760768e22,
146         0.2585201673888497664e23,
147         0.62044840173323943936e24,
148         0.15511210043330985984e26,
149         0.403291461126605635584e27,
150         0.10888869450418352160768e29,
151         0.304888344611713860501504e30,
152         0.8841761993739701954543616e31,
153         0.26525285981219105863630848e33,
154         0.822283865417792281772556288e34,
155         0.26313083693369353016721801216e36,
156         0.868331761881188649551819440128e37,
157         0.29523279903960414084761860964352e39,
158         0.103331479663861449296666513375232e41,
159         0.3719933267899012174679994481508352e42,
160         0.137637530912263450463159795815809024e44,
161         0.5230226174666011117600072241000742912e45,
162         0.203978820811974433586402817399028973568e47,
163         0.815915283247897734345611269596115894272e48,
164         0.3345252661316380710817006205344075166515e50,
165         0.1405006117752879898543142606244511569936e52,
166         0.6041526306337383563735513206851399750726e53,
167         0.265827157478844876804362581101461589032e55,
168         0.1196222208654801945619631614956577150644e57,
169         0.5502622159812088949850305428800254892962e58,
170         0.2586232415111681806429643551536119799692e60,
171         0.1241391559253607267086228904737337503852e62,
172         0.6082818640342675608722521633212953768876e63,
173         0.3041409320171337804361260816606476884438e65,
174         0.1551118753287382280224243016469303211063e67,
175         0.8065817517094387857166063685640376697529e68,
176         0.427488328406002556429801375338939964969e70,
177         0.2308436973392413804720927426830275810833e72,
178         0.1269640335365827592596510084756651695958e74,
179         0.7109985878048634518540456474637249497365e75,
180         0.4052691950487721675568060190543232213498e77,
181         0.2350561331282878571829474910515074683829e79,
182         0.1386831185456898357379390197203894063459e81,
183         0.8320987112741390144276341183223364380754e82,
184         0.507580213877224798800856812176625227226e84,
185         0.3146997326038793752565312235495076408801e86,
186         0.1982608315404440064116146708361898137545e88,
187         0.1268869321858841641034333893351614808029e90,
188         0.8247650592082470666723170306785496252186e91,
189         0.5443449390774430640037292402478427526443e93,
190         0.3647111091818868528824985909660546442717e95,
191         0.2480035542436830599600990418569171581047e97,
192         0.1711224524281413113724683388812728390923e99,
193         0.1197857166996989179607278372168909873646e101,
194         0.8504785885678623175211676442399260102886e102,
195         0.6123445837688608686152407038527467274078e104,
196         0.4470115461512684340891257138125051110077e106,
197         0.3307885441519386412259530282212537821457e108,
198         0.2480914081139539809194647711659403366093e110,
199         0.188549470166605025498793226086114655823e112,
200         0.1451830920282858696340707840863082849837e114,
201         0.1132428117820629783145752115873204622873e116,
202         0.8946182130782975286851441715398316520698e117,
203         0.7156945704626380229481153372318653216558e119,
204         0.5797126020747367985879734231578109105412e121,
205         0.4753643337012841748421382069894049466438e123,
206         0.3945523969720658651189747118012061057144e125,
207         0.3314240134565353266999387579130131288001e127,
208         0.2817104114380550276949479442260611594801e129,
209         0.2422709538367273238176552320344125971528e131,
210         0.210775729837952771721360051869938959523e133,
211         0.1854826422573984391147968456455462843802e135,
212         0.1650795516090846108121691926245361930984e137,
213         0.1485715964481761497309522733620825737886e139,
214         0.1352001527678402962551665687594951421476e141,
215         0.1243841405464130725547532432587355307758e143,
216         0.1156772507081641574759205162306240436215e145,
217         0.1087366156656743080273652852567866010042e147,
218         0.103299784882390592625997020993947270954e149,
219         0.9916779348709496892095714015418938011582e150,
220         0.9619275968248211985332842594956369871234e152,
221         0.942689044888324774562618574305724247381e154,
222         0.9332621544394415268169923885626670049072e156,
223         0.9332621544394415268169923885626670049072e158,
224         0.9425947759838359420851623124482936749562e160,
225         0.9614466715035126609268655586972595484554e162,
226         0.990290071648618040754671525458177334909e164,
227         0.1029901674514562762384858386476504428305e167,
228         0.1081396758240290900504101305800329649721e169,
229         0.1146280563734708354534347384148349428704e171,
230         0.1226520203196137939351751701038733888713e173,
231         0.132464181945182897449989183712183259981e175,
232         0.1443859583202493582204882102462797533793e177,
233         0.1588245541522742940425370312709077287172e179,
234         0.1762952551090244663872161047107075788761e181,
235         0.1974506857221074023536820372759924883413e183,
236         0.2231192748659813646596607021218715118256e185,
237         0.2543559733472187557120132004189335234812e187,
238         0.2925093693493015690688151804817735520034e189,
239         0.339310868445189820119825609358857320324e191,
240         0.396993716080872089540195962949863064779e193,
241         0.4684525849754290656574312362808384164393e195,
242         0.5574585761207605881323431711741977155627e197,
243         0.6689502913449127057588118054090372586753e199,
244         0.8094298525273443739681622845449350829971e201,
245         0.9875044200833601362411579871448208012564e203,
246         0.1214630436702532967576624324188129585545e206,
247         0.1506141741511140879795014161993280686076e208,
248         0.1882677176888926099743767702491600857595e210,
249         0.237217324288004688567714730513941708057e212,
250         0.3012660018457659544809977077527059692324e214,
251         0.3856204823625804217356770659234636406175e216,
252         0.4974504222477287440390234150412680963966e218,
253         0.6466855489220473672507304395536485253155e220,
254         0.8471580690878820510984568758152795681634e222,
255         0.1118248651196004307449963076076169029976e225,
256         0.1487270706090685728908450891181304809868e227,
257         0.1992942746161518876737324194182948445223e229,
258         0.269047270731805048359538766214698040105e231,
259         0.3659042881952548657689727220519893345429e233,
260         0.5012888748274991661034926292112253883237e235,
261         0.6917786472619488492228198283114910358867e237,
262         0.9615723196941089004197195613529725398826e239,
263         0.1346201247571752460587607385894161555836e242,
264         0.1898143759076170969428526414110767793728e244,
265         0.2695364137888162776588507508037290267094e246,
266         0.3854370717180072770521565736493325081944e248,
267         0.5550293832739304789551054660550388118e250,
268         0.80479260574719919448490292577980627711e252,
269         0.1174997204390910823947958271638517164581e255,
270         0.1727245890454638911203498659308620231933e257,
271         0.2556323917872865588581178015776757943262e259,
272         0.380892263763056972698595524350736933546e261,
273         0.571338395644585459047893286526105400319e263,
274         0.8627209774233240431623188626544191544816e265,
275         0.1311335885683452545606724671234717114812e268,
276         0.2006343905095682394778288746989117185662e270,
277         0.308976961384735088795856467036324046592e272,
278         0.4789142901463393876335775239063022722176e274,
279         0.7471062926282894447083809372938315446595e276,
280         0.1172956879426414428192158071551315525115e279,
281         0.1853271869493734796543609753051078529682e281,
282         0.2946702272495038326504339507351214862195e283,
283         0.4714723635992061322406943211761943779512e285,
284         0.7590705053947218729075178570936729485014e287,
285         0.1229694218739449434110178928491750176572e290,
286         0.2004401576545302577599591653441552787813e292,
287         0.3287218585534296227263330311644146572013e294,
288         0.5423910666131588774984495014212841843822e296,
289         0.9003691705778437366474261723593317460744e298,
290         0.1503616514864999040201201707840084015944e301,
291         0.2526075744973198387538018869171341146786e303,
292         0.4269068009004705274939251888899566538069e305,
293         0.7257415615307998967396728211129263114717e307,
294     };
295 
296     /**
297      * 53-bit precision implementation of the Lanczos approximation.
298      *
299      * <p>This implementation is in partial fraction form with the leading constant
300      * of \( \sqrt{2\pi} \) absorbed into the sum.
301      *
302      * <p>It is related to the Gamma function by the following equation
303      * \[
304      * \Gamma(z) = \frac{(z + g - 0.5)^{z - 0.5}}{e^{z + g - 0.5}} \mathrm{lanczos}(z)
305      * \]
306      * where \( g \) is the Lanczos constant.
307      *
308      * <h2>Warning</h2>
309      *
310      * <p>This is not a substitute for {@link LanczosApproximation}. The approximation is
311      * written in partial fraction form with the leading constants absorbed by the
312      * coefficients in the sum.
313      *
314      * @see <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/lanczos.html">
315      * Boost Lanczos Approximation</a>
316      */
317     static final class Lanczos {
318         // Optimal values for G for each N are taken from
319         // http://web.mala.bc.ca/pughg/phdThesis/phdThesis.pdf,
320         // as are the theoretical error bounds.
321         //
322         // Constants calculated using the method described by Godfrey
323         // http://my.fit.edu/~gabdo/gamma.txt and elaborated by Toth at
324         // http://www.rskey.org/gamma.htm using NTL::RR at 1000 bit precision.
325 
326         //
327         // Lanczos Coefficients for N=13 G=6.024680040776729583740234375
328         // Max experimental error (with arbitrary precision arithmetic) 1.196214e-17
329         // Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
330         //
331 
332         /**
333          * Lanczos constant G.
334          */
335         static final double G = 6.024680040776729583740234375;
336 
337         /**
338          * Lanczos constant G - half.
339          *
340          * <p>Note: The form {@code (g - 0.5)} is used when computing the gamma function.
341          */
342         static final double GMH = 5.524680040776729583740234375;
343 
344         /** Common denominator used for the rational evaluation. */
345         private static final int[] DENOM = {
346             0,
347             39916800,
348             120543840,
349             150917976,
350             105258076,
351             45995730,
352             13339535,
353             2637558,
354             357423,
355             32670,
356             1925,
357             66,
358             1
359         };
360 
361         /** Private constructor. */
362         private Lanczos() {
363             // intentionally empty.
364         }
365 
366         /**
367          * Computes the Lanczos approximation.
368          *
369          * @param z Argument.
370          * @return the Lanczos approximation.
371          */
372         static double lanczosSum(double z) {
373             final double[] num = {
374                 23531376880.41075968857200767445163675473,
375                 42919803642.64909876895789904700198885093,
376                 35711959237.35566804944018545154716670596,
377                 17921034426.03720969991975575445893111267,
378                 6039542586.35202800506429164430729792107,
379                 1439720407.311721673663223072794912393972,
380                 248874557.8620541565114603864132294232163,
381                 31426415.58540019438061423162831820536287,
382                 2876370.628935372441225409051620849613599,
383                 186056.2653952234950402949897160456992822,
384                 8071.672002365816210638002902272250613822,
385                 210.8242777515793458725097339207133627117,
386                 2.506628274631000270164908177133837338626
387             };
388             return evaluateRational(num, DENOM, z);
389         }
390 
391         /**
392          * Computes the Lanczos approximation scaled by {@code exp(g)}.
393          *
394          * @param z Argument.
395          * @return the scaled Lanczos approximation.
396          */
397         static double lanczosSumExpGScaled(double z) {
398             // As above with numerator divided by exp(g) = 413.509...
399             final double[] num = {
400                 56906521.91347156388090791033559122686859,
401                 103794043.1163445451906271053616070238554,
402                 86363131.28813859145546927288977868422342,
403                 43338889.32467613834773723740590533316085,
404                 14605578.08768506808414169982791359218571,
405                 3481712.15498064590882071018964774556468,
406                 601859.6171681098786670226533699352302507,
407                 75999.29304014542649875303443598909137092,
408                 6955.999602515376140356310115515198987526,
409                 449.9445569063168119446858607650988409623,
410                 19.51992788247617482847860966235652136208,
411                 0.5098416655656676188125178644804694509993,
412                 0.006061842346248906525783753964555936883222
413             };
414             return evaluateRational(num, DENOM, z);
415         }
416 
417         /**
418          * Evaluate the rational number as two polynomials.
419          *
420          * <p>Adapted from {@code boost/math/tools/detail/rational_horner3_13.hpp}.
421          * Note: There are 3 variations of the unrolled rational evaluation.
422          * These methods change the order based on the sign of x. This
423          * should be used for the Lanczos code as this comment in
424          * {@code boost/math/tools/rational.hpp} notes:
425          *
426          * <blockquote>
427          * However, there
428          * are some tricks we can use to prevent overflow that might otherwise
429          * occur in polynomial evaluation, if z is large.  This is important
430          * in our Lanczos code for example.
431          * </blockquote>
432          *
433          * @param a Coefficients of the numerator polynomial
434          * @param b Coefficients of the denominator polynomial
435          * @param x Value
436          * @return the rational number
437          */
438         private static double evaluateRational(double[] a, int[] b, double x) {
439             // The choice of algorithm in Boost is based on the compiler
440             // to suite the available optimisations.
441             //
442             // Tests against rational_horner1_13.hpp which uses a first order
443             // Horner method (no x*x term) show only minor variations in
444             // error. rational_horner2_13.hpp has the same second order Horner
445             // method with different code layout of the same sum.
446 
447             // rational_horner3_13.hpp
448             if (x <= 1) {
449                 final double x2 = x * x;
450                 double t0 = a[12] * x2 + a[10];
451                 double t1 = a[11] * x2 + a[9];
452                 double t2 = b[12] * x2 + b[10];
453                 double t3 = b[11] * x2 + b[9];
454                 t0 *= x2;
455                 t1 *= x2;
456                 t2 *= x2;
457                 t3 *= x2;
458                 t0 += a[8];
459                 t1 += a[7];
460                 t2 += b[8];
461                 t3 += b[7];
462                 t0 *= x2;
463                 t1 *= x2;
464                 t2 *= x2;
465                 t3 *= x2;
466                 t0 += a[6];
467                 t1 += a[5];
468                 t2 += b[6];
469                 t3 += b[5];
470                 t0 *= x2;
471                 t1 *= x2;
472                 t2 *= x2;
473                 t3 *= x2;
474                 t0 += a[4];
475                 t1 += a[3];
476                 t2 += b[4];
477                 t3 += b[3];
478                 t0 *= x2;
479                 t1 *= x2;
480                 t2 *= x2;
481                 t3 *= x2;
482                 t0 += a[2];
483                 t1 += a[1];
484                 t2 += b[2];
485                 t3 += b[1];
486                 t0 *= x2;
487                 t2 *= x2;
488                 t0 += a[0];
489                 t2 += b[0];
490                 t1 *= x;
491                 t3 *= x;
492                 return (t0 + t1) / (t2 + t3);
493             }
494             final double z = 1 / x;
495             final double z2 = 1 / (x * x);
496             double t0 = a[0] * z2 + a[2];
497             double t1 = a[1] * z2 + a[3];
498             double t2 = b[0] * z2 + b[2];
499             double t3 = b[1] * z2 + b[3];
500             t0 *= z2;
501             t1 *= z2;
502             t2 *= z2;
503             t3 *= z2;
504             t0 += a[4];
505             t1 += a[5];
506             t2 += b[4];
507             t3 += b[5];
508             t0 *= z2;
509             t1 *= z2;
510             t2 *= z2;
511             t3 *= z2;
512             t0 += a[6];
513             t1 += a[7];
514             t2 += b[6];
515             t3 += b[7];
516             t0 *= z2;
517             t1 *= z2;
518             t2 *= z2;
519             t3 *= z2;
520             t0 += a[8];
521             t1 += a[9];
522             t2 += b[8];
523             t3 += b[9];
524             t0 *= z2;
525             t1 *= z2;
526             t2 *= z2;
527             t3 *= z2;
528             t0 += a[10];
529             t1 += a[11];
530             t2 += b[10];
531             t3 += b[11];
532             t0 *= z2;
533             t2 *= z2;
534             t0 += a[12];
535             t2 += b[12];
536             t1 *= z;
537             t3 *= z;
538             return (t0 + t1) / (t2 + t3);
539         }
540 
541         // Not implemented:
542         // lanczos_sum_near_1
543         // lanczos_sum_near_2
544     }
545 
546     /** Private constructor. */
547     private BoostGamma() {
548         // intentionally empty.
549     }
550 
551     /**
552      * All factorials that are representable as a double.
553      * This data is exposed for testing.
554      *
555      * @return factorials
556      */
557     static double[] getFactorials() {
558         return FACTORIAL.clone();
559     }
560 
561     /**
562      * Returns the factorial of n.
563      * This is unchecked as an index out of bound exception will occur
564      * if the value n is not within the range [0, 170].
565      * This function is exposed for use in {@link BoostBeta}.
566      *
567      * @param n Argument n (must be in [0, 170])
568      * @return n!
569      */
570     static double uncheckedFactorial(int n) {
571         return FACTORIAL[n];
572     }
573 
574     /**
575      * Gamma function.
576      *
577      * <p>For small {@code z} this is based on the <em>NSWC Library of Mathematics
578      * Subroutines</em> double precision implementation, {@code DGAMMA}.
579      *
580      * <p>For large {@code z} this is an implementation of the Boost C++ tgamma
581      * function with Lanczos support.
582      *
583      * <p>Integers are handled using a look-up table of factorials.
584      *
585      * <p>Note: The Boost C++ implementation uses the Lanczos sum for all {@code z}.
586      * When promotion of double to long double is not available this has larger
587      * errors than the double precision specific NSWC implementation. For larger
588      * {@code z} the Boost C++ Lanczos implementation incorporates the sqrt(2 pi)
589      * factor and has lower error than the implementation using the
590      * {@link LanczosApproximation} class.
591      *
592      * @param z Argument z
593      * @return gamma value
594      */
595     static double tgamma(double z) {
596         // Handle integers
597         if (Math.rint(z) == z) {
598             if (z <= 0) {
599                 // Pole error
600                 return Double.NaN;
601             }
602             if (z <= MAX_GAMMA_Z) {
603                 // Gamma(n) = (n-1)!
604                 return FACTORIAL[(int) z - 1];
605             }
606             // Overflow
607             return Double.POSITIVE_INFINITY;
608         }
609 
610         if (Math.abs(z) <= LANCZOS_THRESHOLD) {
611             // Small z
612             // NSWC Library of Mathematics Subroutines
613             // Note:
614             // This does not benefit from using extended precision to track the sum (t).
615             // Extended precision on the product reduces the error but the majority
616             // of error remains in InvGamma1pm1.
617 
618             if (z >= 1) {
619                 /*
620                  * From the recurrence relation
621                  * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
622                  * then
623                  * Gamma(t) = 1 / [1 + InvGamma1pm1.value(t - 1)],
624                  * where t = x - n. This means that t must satisfy
625                  * -0.5 <= t - 1 <= 1.5.
626                  */
627                 double prod = 1;
628                 double t = z;
629                 while (t > 2.5) {
630                     t -= 1;
631                     prod *= t;
632                 }
633                 return prod / (1 + InvGamma1pm1.value(t - 1));
634             }
635             /*
636              * From the recurrence relation
637              * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
638              * then
639              * Gamma(x + n + 1) = 1 / [1 + InvGamma1pm1.value(x + n)],
640              * which requires -0.5 <= x + n <= 1.5.
641              */
642             double prod = z;
643             double t = z;
644             while (t < -0.5) {
645                 t += 1;
646                 prod *= t;
647             }
648             return 1 / (prod * (1 + InvGamma1pm1.value(t)));
649         }
650 
651         // Large non-integer z
652         // Boost C++ tgamma implementation
653 
654         if (z < 0) {
655             /*
656              * From the reflection formula
657              * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
658              * and the recurrence relation
659              * Gamma(1 - x) = -x * Gamma(-x),
660              * it is found
661              * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
662              */
663             return -Math.PI / (sinpx(z) * tgamma(-z));
664         } else if (z > MAX_GAMMA_Z + 1) {
665             // Addition to the Boost code: Simple overflow detection
666             return Double.POSITIVE_INFINITY;
667         }
668 
669         double result = Lanczos.lanczosSum(z);
670         final double zgh = z + Lanczos.GMH;
671         final double lzgh = Math.log(zgh);
672         if (z * lzgh > LOG_MAX_VALUE) {
673             // we're going to overflow unless this is done with care:
674 
675             // Updated
676             // Check for overflow removed:
677             // if (lzgh * z / 2 > LOG_MAX_VALUE) ... overflow
678             // This is replaced by checking z > MAX_FACTORIAL + 2
679 
680             final double hp = Math.pow(zgh, (z / 2) - 0.25);
681             result *= hp / Math.exp(zgh);
682             // Check for overflow has been removed:
683             // if (Double.MAX_VALUE / hp < result) ... overflow
684             result *= hp;
685         } else {
686             result *= Math.pow(zgh, z - 0.5) / Math.exp(zgh);
687         }
688 
689         return result;
690     }
691 
692     /**
693      * Ad hoc function calculates x * sin(pi * x), taking extra care near when x is
694      * near a whole number.
695      *
696      * @param x Value (assumed to be negative)
697      * @return x * sin(pi * x)
698      */
699     static double sinpx(double x) {
700         int sign = 1;
701         // This is always called with a negative
702         // if (x < 0)
703         x = -x;
704         double fl = Math.floor(x);
705         double dist;
706         if (isOdd(fl)) {
707             fl += 1;
708             dist = fl - x;
709             sign = -sign;
710         } else {
711             dist = x - fl;
712         }
713         if (dist > 0.5f) {
714             dist = 1 - dist;
715         }
716         final double result = Math.sin(dist * Math.PI);
717         return sign * x * result;
718     }
719 
720     /**
721      * Checks if the value is odd.
722      *
723      * @param v Value (assumed to be positive and an integer)
724      * @return true if odd
725      */
726     private static boolean isOdd(double v) {
727         // Note:
728         // Any value larger than 2^53 should be even.
729         // If the input is positive then truncation of extreme doubles (>2^63)
730         // to the primitive long creates an odd value: 2^63-1.
731         // This is corrected by inverting the sign of v and the extreme is even: -2^63.
732         // This function is never called when the argument is this large
733         // as this is a pole error in tgamma so the effect is never observed.
734         // However the isOdd function is correct for all positive finite v.
735         return (((long) -v) & 0x1) == 1;
736     }
737 
738     /**
739      * Log Gamma function.
740      * Defined as the natural logarithm of the absolute value of tgamma(z).
741      *
742      * @param z Argument z
743      * @return log gamma value
744      */
745     static double lgamma(double z) {
746         return lgamma(z, null);
747     }
748 
749     /**
750      * Log Gamma function.
751      * Defined as the natural logarithm of the absolute value of tgamma(z).
752      *
753      * @param z Argument z
754      * @param sign If a non-zero length array the first index is set on output to the sign of tgamma(z)
755      * @return log gamma value
756      */
757     static double lgamma(double z, int[] sign) {
758         double result;
759         int sresult = 1;
760         if (z <= -ROOT_EPSILON) {
761             // reflection formula:
762             if (Math.rint(z) == z) {
763                 // Pole error
764                 return Double.NaN;
765             }
766 
767             double t = sinpx(z);
768             z = -z;
769             if (t < 0) {
770                 t = -t;
771             } else {
772                 sresult = -sresult;
773             }
774 
775             // This summation can have large magnitudes with opposite signs.
776             // Use an extended precision sum to reduce cancellation.
777             result = Sum.of(-lgamma(z)).add(-Math.log(t)).add(LOG_PI).getAsDouble();
778 
779         } else if (z < ROOT_EPSILON) {
780             if (z == 0) {
781                 // Pole error
782                 return Double.NaN;
783             }
784             if (4 * Math.abs(z) < EPSILON) {
785                 result = -Math.log(Math.abs(z));
786             } else {
787                 result = Math.log(Math.abs(1 / z - EULER));
788             }
789             if (z < 0) {
790                 sresult = -1;
791             }
792         } else if (z < 15) {
793             result = lgammaSmall(z, z - 1, z - 2);
794         // The z > 3 condition is always true
795         //} else if (z > 3 && z < 100) {
796         } else if (z < 100) {
797             // taking the log of tgamma reduces the error, no danger of overflow here:
798             result = Math.log(tgamma(z));
799         } else {
800             // regular evaluation:
801             final double zgh = z + Lanczos.GMH;
802             result = Math.log(zgh) - 1;
803             result *= z - 0.5f;
804             //
805             // Only add on the lanczos sum part if we're going to need it:
806             //
807             if (result * EPSILON < 20) {
808                 result += Math.log(Lanczos.lanczosSumExpGScaled(z));
809             }
810         }
811 
812         if (nonZeroLength(sign)) {
813             sign[0] = sresult;
814         }
815         return result;
816     }
817 
818     /**
819      * Log Gamma function for small z.
820      *
821      * @param z Argument z
822      * @param zm1 {@code z - 1}
823      * @param zm2 {@code z - 2}
824      * @return log gamma value
825      */
826     private static double lgammaSmall(double z, double zm1, double zm2) {
827         // This version uses rational approximations for small
828         // values of z accurate enough for 64-bit mantissas
829         // (80-bit long doubles), works well for 53-bit doubles as well.
830 
831         // Updated to use an extended precision sum
832         final Sum result = Sum.create();
833 
834         // Note:
835         // Removed z < EPSILON branch.
836         // The function is called
837         // from lgamma:
838         //   ROOT_EPSILON <= z < 15
839         // from tgamma1pm1:
840         //   1.5 <= z < 2
841         //   1 <= z < 3
842 
843         if ((zm1 == 0) || (zm2 == 0)) {
844             // nothing to do, result is zero....
845             return 0;
846         } else if (z > 2) {
847             //
848             // Begin by performing argument reduction until
849             // z is in [2,3):
850             //
851             if (z >= 3) {
852                 do {
853                     z -= 1;
854                     result.add(Math.log(z));
855                 } while (z >= 3);
856                 // Update zm2, we need it below:
857                 zm2 = z - 2;
858             }
859 
860             //
861             // Use the following form:
862             //
863             // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
864             //
865             // where R(z-2) is a rational approximation optimised for
866             // low absolute error - as long as its absolute error
867             // is small compared to the constant Y - then any rounding
868             // error in its computation will get wiped out.
869             //
870             // R(z-2) has the following properties:
871             //
872             // At double: Max error found:                    4.231e-18
873             // At long double: Max error found:               1.987e-21
874             // Maximum Deviation Found (approximation error): 5.900e-24
875             //
876             double P;
877             P = -0.324588649825948492091e-4;
878             P = -0.541009869215204396339e-3 + P * zm2;
879             P = -0.259453563205438108893e-3 + P * zm2;
880             P =  0.172491608709613993966e-1 + P * zm2;
881             P =  0.494103151567532234274e-1 + P * zm2;
882             P =   0.25126649619989678683e-1 + P * zm2;
883             P = -0.180355685678449379109e-1 + P * zm2;
884             double Q;
885             Q = -0.223352763208617092964e-6;
886             Q =  0.224936291922115757597e-3 + Q * zm2;
887             Q =   0.82130967464889339326e-2 + Q * zm2;
888             Q =  0.988504251128010129477e-1 + Q * zm2;
889             Q =   0.541391432071720958364e0 + Q * zm2;
890             Q =   0.148019669424231326694e1 + Q * zm2;
891             Q =   0.196202987197795200688e1 + Q * zm2;
892             Q =                       0.1e1 + Q * zm2;
893 
894             final float Y = 0.158963680267333984375e0f;
895 
896             final double r = zm2 * (z + 1);
897             final double R = P / Q;
898 
899             result.addProduct(r, Y).addProduct(r, R);
900         } else {
901             //
902             // If z is less than 1 use recurrence to shift to
903             // z in the interval [1,2]:
904             //
905             if (z < 1) {
906                 result.add(-Math.log(z));
907                 zm2 = zm1;
908                 zm1 = z;
909                 z += 1;
910             }
911             //
912             // Two approximations, one for z in [1,1.5] and
913             // one for z in [1.5,2]:
914             //
915             if (z <= 1.5) {
916                 //
917                 // Use the following form:
918                 //
919                 // lgamma(z) = (z-1)(z-2)(Y + R(z-1
920                 //
921                 // where R(z-1) is a rational approximation optimised for
922                 // low absolute error - as long as its absolute error
923                 // is small compared to the constant Y - then any rounding
924                 // error in its computation will get wiped out.
925                 //
926                 // R(z-1) has the following properties:
927                 //
928                 // At double precision: Max error found:                1.230011e-17
929                 // At 80-bit long double precision:   Max error found:  5.631355e-21
930                 // Maximum Deviation Found:                             3.139e-021
931                 // Expected Error Term:                                 3.139e-021
932 
933                 //
934                 final float Y = 0.52815341949462890625f;
935 
936                 double P;
937                 P = -0.100346687696279557415e-2;
938                 P = -0.240149820648571559892e-1 + P * zm1;
939                 P =  -0.158413586390692192217e0 + P * zm1;
940                 P =  -0.406567124211938417342e0 + P * zm1;
941                 P =  -0.414983358359495381969e0 + P * zm1;
942                 P = -0.969117530159521214579e-1 + P * zm1;
943                 P =  0.490622454069039543534e-1 + P * zm1;
944                 double Q;
945                 Q = 0.195768102601107189171e-2;
946                 Q = 0.577039722690451849648e-1 + Q * zm1;
947                 Q =  0.507137738614363510846e0 + Q * zm1;
948                 Q =  0.191415588274426679201e1 + Q * zm1;
949                 Q =  0.348739585360723852576e1 + Q * zm1;
950                 Q =  0.302349829846463038743e1 + Q * zm1;
951                 Q =                      0.1e1 + Q * zm1;
952 
953                 final double r = P / Q;
954                 final double prefix = zm1 * zm2;
955 
956                 result.addProduct(prefix, Y).addProduct(prefix, r);
957             } else {
958                 //
959                 // Use the following form:
960                 //
961                 // lgamma(z) = (2-z)(1-z)(Y + R(2-z
962                 //
963                 // where R(2-z) is a rational approximation optimised for
964                 // low absolute error - as long as its absolute error
965                 // is small compared to the constant Y - then any rounding
966                 // error in its computation will get wiped out.
967                 //
968                 // R(2-z) has the following properties:
969                 //
970                 // At double precision, max error found:              1.797565e-17
971                 // At 80-bit long double precision, max error found:  9.306419e-21
972                 // Maximum Deviation Found:                           2.151e-021
973                 // Expected Error Term:                               2.150e-021
974                 //
975                 final float Y = 0.452017307281494140625f;
976 
977                 final double mzm2 = -zm2;
978                 double P;
979                 P =  0.431171342679297331241e-3;
980                 P = -0.850535976868336437746e-2 + P * mzm2;
981                 P =  0.542809694055053558157e-1 + P * mzm2;
982                 P =  -0.142440390738631274135e0 + P * mzm2;
983                 P =   0.144216267757192309184e0 + P * mzm2;
984                 P = -0.292329721830270012337e-1 + P * mzm2;
985                 double Q;
986                 Q = -0.827193521891290553639e-6;
987                 Q = -0.100666795539143372762e-2 + Q * mzm2;
988                 Q =   0.25582797155975869989e-1 + Q * mzm2;
989                 Q =  -0.220095151814995745555e0 + Q * mzm2;
990                 Q =   0.846973248876495016101e0 + Q * mzm2;
991                 Q =  -0.150169356054485044494e1 + Q * mzm2;
992                 Q =                       0.1e1 + Q * mzm2;
993                 final double r = zm2 * zm1;
994                 final double R = P / Q;
995 
996                 result.addProduct(r, Y).addProduct(r, R);
997             }
998         }
999         return result.getAsDouble();
1000     }
1001 
1002     /**
1003      * Calculates tgamma(1+dz)-1.
1004      *
1005      * @param dz Argument
1006      * @return tgamma(1+dz)-1
1007      */
1008     static double tgamma1pm1(double dz) {
1009         //
1010         // This helper calculates tgamma(1+dz)-1 without cancellation errors,
1011         // used by the upper incomplete gamma with z < 1:
1012         //
1013         final double result;
1014         if (dz < 0) {
1015             if (dz < -0.5) {
1016                 // Best method is simply to subtract 1 from tgamma:
1017                 result = tgamma(1 + dz) - 1;
1018             } else {
1019                 // Use expm1 on lgamma:
1020                 result = Math.expm1(-Math.log1p(dz) + lgammaSmall(dz + 2, dz + 1, dz));
1021             }
1022         } else {
1023             if (dz < 2) {
1024                 // Use expm1 on lgamma:
1025                 result = Math.expm1(lgammaSmall(dz + 1, dz, dz - 1));
1026             } else {
1027                 // Best method is simply to subtract 1 from tgamma:
1028                 result = tgamma(1 + dz) - 1;
1029             }
1030         }
1031 
1032         return result;
1033     }
1034 
1035     /**
1036      * Full upper incomplete gamma.
1037      *
1038      * @param a Argument a
1039      * @param x Argument x
1040      * @return upper gamma value
1041      */
1042     static double tgamma(double a, double x) {
1043         return gammaIncompleteImp(a, x, false, true, Policy.getDefault());
1044     }
1045 
1046     /**
1047      * Full upper incomplete gamma.
1048      *
1049      * @param a Argument a
1050      * @param x Argument x
1051      * @param policy Function evaluation policy
1052      * @return upper gamma value
1053      */
1054     static double tgamma(double a, double x, Policy policy) {
1055         return gammaIncompleteImp(a, x, false, true, policy);
1056     }
1057 
1058     /**
1059      * Full lower incomplete gamma.
1060      *
1061      * @param a Argument a
1062      * @param x Argument x
1063      * @return lower gamma value
1064      */
1065     static double tgammaLower(double a, double x) {
1066         return gammaIncompleteImp(a, x, false, false, Policy.getDefault());
1067     }
1068 
1069     /**
1070      * Full lower incomplete gamma.
1071      *
1072      * @param a Argument a
1073      * @param x Argument x
1074      * @param policy Function evaluation policy
1075      * @return lower gamma value
1076      */
1077     static double tgammaLower(double a, double x, Policy policy) {
1078         return gammaIncompleteImp(a, x, false, false, policy);
1079     }
1080 
1081     /**
1082      * Regularised upper incomplete gamma.
1083      *
1084      * @param a Argument a
1085      * @param x Argument x
1086      * @return q
1087      */
1088     static double gammaQ(double a, double x) {
1089         return gammaIncompleteImp(a, x, true, true, Policy.getDefault());
1090     }
1091 
1092     /**
1093      * Regularised upper incomplete gamma.
1094      *
1095      * @param a Argument a
1096      * @param x Argument x
1097      * @param policy Function evaluation policy
1098      * @return q
1099      */
1100     static double gammaQ(double a, double x, Policy policy) {
1101         return gammaIncompleteImp(a, x, true, true, policy);
1102     }
1103 
1104     /**
1105      * Regularised lower incomplete gamma.
1106      *
1107      * @param a Argument a
1108      * @param x Argument x
1109      * @return p
1110      */
1111     static double gammaP(double a, double x) {
1112         return gammaIncompleteImp(a, x, true, false, Policy.getDefault());
1113     }
1114 
1115     /**
1116      * Regularised lower incomplete gamma.
1117      *
1118      * @param a Argument a
1119      * @param x Argument x
1120      * @param policy Function evaluation policy
1121      * @return p
1122      */
1123     static double gammaP(double a, double x, Policy policy) {
1124         return gammaIncompleteImp(a, x, true, false, policy);
1125     }
1126 
1127     /**
1128      * Derivative of the regularised lower incomplete gamma.
1129      * <p>\( \frac{e^{-x} x^{a-1}}{\Gamma(a)} \)
1130      *
1131      * <p>Adapted from {@code boost::math::detail::gamma_p_derivative_imp}
1132      *
1133      * @param a Argument a
1134      * @param x Argument x
1135      * @return p derivative
1136      */
1137     static double gammaPDerivative(double a, double x) {
1138         //
1139         // Usual error checks first:
1140         //
1141         if (Double.isNaN(a) || Double.isNaN(x) || a <= 0 || x < 0) {
1142             return Double.NaN;
1143         }
1144         //
1145         // Now special cases:
1146         //
1147         if (x == 0) {
1148             if (a > 1) {
1149                 return 0;
1150             }
1151             return (a == 1) ? 1 : Double.POSITIVE_INFINITY;
1152         }
1153         //
1154         // Normal case:
1155         //
1156         double f1 = regularisedGammaPrefix(a, x);
1157         if (f1 == 0) {
1158             // Underflow in calculation, use logs instead:
1159             f1 = a * Math.log(x) - x - lgamma(a) - Math.log(x);
1160             f1 = Math.exp(f1);
1161         } else {
1162             // Will overflow when (x < 1) && (Double.MAX_VALUE * x < f1).
1163             // There is no exception for this case so just return the result.
1164             f1 /= x;
1165         }
1166 
1167         return f1;
1168     }
1169 
1170     /**
1171      * Main incomplete gamma entry point, handles all four incomplete gammas.
1172      * Adapted from {@code boost::math::detail::gamma_incomplete_imp}.
1173      *
1174      * <p>The Boost code has a pointer {@code p_derivative} that can be set to the
1175      * value of the derivative. This is used for the inverse incomplete
1176      * gamma functions {@code gamma_(p|q)_inv_imp}. It is not required for the forward
1177      * evaluation functions.
1178      *
1179      * @param a Argument a
1180      * @param x Argument x
1181      * @param normalised true to compute the regularised value
1182      * @param invert true to compute the upper value Q (default is lower value P)
1183      * @param pol Function evaluation policy
1184      * @return gamma value
1185      */
1186     private static double gammaIncompleteImp(double a, double x,
1187             boolean normalised, boolean invert, Policy pol) {
1188         if (Double.isNaN(a) || Double.isNaN(x) || a <= 0 || x < 0) {
1189             return Double.NaN;
1190         }
1191 
1192         double result = 0;
1193 
1194         if (a >= MAX_FACTORIAL && !normalised) {
1195             //
1196             // When we're computing the non-normalized incomplete gamma
1197             // and a is large the result is rather hard to compute unless
1198             // we use logs. There are really two options - if x is a long
1199             // way from a in value then we can reliably use methods 2 and 4
1200             // below in logarithmic form and go straight to the result.
1201             // Otherwise we let the regularized gamma take the strain
1202             // (the result is unlikely to underflow in the central region anyway)
1203             // and combine with lgamma in the hopes that we get a finite result.
1204             //
1205 
1206             if (invert && (a * 4 < x)) {
1207                 // This is method 4 below, done in logs:
1208                 result = a * Math.log(x) - x;
1209                 result += Math.log(upperGammaFraction(a, x, pol));
1210             } else if (!invert && (a > 4 * x)) {
1211                 // This is method 2 below, done in logs:
1212                 result = a * Math.log(x) - x;
1213                 result += Math.log(lowerGammaSeries(a, x, 0, pol) / a);
1214             } else {
1215                 result = gammaIncompleteImp(a, x, true, invert, pol);
1216                 if (result == 0) {
1217                     if (invert) {
1218                         // Try http://functions.wolfram.com/06.06.06.0039.01
1219                         result = 1 + 1 / (12 * a) + 1 / (288 * a * a);
1220                         result = Math.log(result) - a + (a - 0.5f) * Math.log(a) + LOG_ROOT_TWO_PI;
1221                     } else {
1222                         // This is method 2 below, done in logs, we're really outside the
1223                         // range of this method, but since the result is almost certainly
1224                         // infinite, we should probably be OK:
1225                         result = a * Math.log(x) - x;
1226                         result += Math.log(lowerGammaSeries(a, x, 0, pol) / a);
1227                     }
1228                 } else {
1229                     result = Math.log(result) + lgamma(a);
1230                 }
1231             }
1232             // If result is > log(MAX_VALUE) the result will overflow.
1233             // There is no exception for this case so just return the result.
1234             return Math.exp(result);
1235         }
1236 
1237         final boolean isInt;
1238         final boolean isHalfInt;
1239         // Update. x must be safe for exp(-x). Change to -x > LOG_MIN_VALUE.
1240         final boolean isSmallA = (a < 30) && (a <= x + 1) && (-x > LOG_MIN_VALUE);
1241         if (isSmallA) {
1242             final double fa = Math.floor(a);
1243             isInt = fa == a;
1244             isHalfInt = !isInt && (Math.abs(fa - a) == 0.5f);
1245         } else {
1246             isInt = isHalfInt = false;
1247         }
1248 
1249         final int evalMethod;
1250 
1251         if (isInt && (x > 0.6)) {
1252             // calculate Q via finite sum:
1253             invert = !invert;
1254             evalMethod = 0;
1255         } else if (isHalfInt && (x > 0.2)) {
1256             // calculate Q via finite sum for half integer a:
1257             invert = !invert;
1258             evalMethod = 1;
1259         } else if ((x < ROOT_EPSILON) && (a > 1)) {
1260             evalMethod = 6;
1261         } else if ((x > 1000) && (a < x * 0.75f)) {
1262             // Note:
1263             // The branch is used in Boost when:
1264             // ((x > 1000) && ((a < x) || (Math.abs(a - 50) / x < 1)))
1265             //
1266             // This case was added after Boost 1_68_0.
1267             // See: https://github.com/boostorg/math/issues/168
1268             //
1269             // When using only double precision for the evaluation
1270             // it is a source of error when a ~ z as the asymptotic approximation
1271             // sums terms t_n+1 = t_n * (a - n - 1) / z starting from t_0 = 1.
1272             // These terms are close to 1 when a ~ z and the sum has many terms
1273             // with reduced precision.
1274             // This has been updated to allow only cases with fast convergence.
1275             // It will be used when x -> infinity and a << x.
1276 
1277             // calculate Q via asymptotic approximation:
1278             invert = !invert;
1279             evalMethod = 7;
1280 
1281         } else if (x < 0.5) {
1282             //
1283             // Changeover criterion chosen to give a changeover at Q ~ 0.33
1284             //
1285             if (-0.4 / Math.log(x) < a) {
1286                 // Compute P
1287                 evalMethod = 2;
1288             } else {
1289                 evalMethod = 3;
1290             }
1291         } else if (x < 1.1) {
1292             //
1293             // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
1294             //
1295             if (x * 0.75f < a) {
1296                 // Compute P
1297                 evalMethod = 2;
1298             } else {
1299                 evalMethod = 3;
1300             }
1301         } else {
1302             //
1303             // Begin by testing whether we're in the "bad" zone
1304             // where the result will be near 0.5 and the usual
1305             // series and continued fractions are slow to converge:
1306             //
1307             boolean useTemme = false;
1308             if (normalised && (a > 20)) {
1309                 final double sigma = Math.abs((x - a) / a);
1310                 if (a > 200) {
1311                     //
1312                     // This limit is chosen so that we use Temme's expansion
1313                     // only if the result would be larger than about 10^-6.
1314                     // Below that the regular series and continued fractions
1315                     // converge OK, and if we use Temme's method we get increasing
1316                     // errors from the dominant erfc term as its (inexact) argument
1317                     // increases in magnitude.
1318                     //
1319                     if (20 / a > sigma * sigma) {
1320                         useTemme = true;
1321                     }
1322                 } else {
1323                     // Note in this zone we can't use Temme's expansion for
1324                     // types longer than an 80-bit real:
1325                     // it would require too many terms in the polynomials.
1326                     if (sigma < 0.4) {
1327                         useTemme = true;
1328                     }
1329                 }
1330             }
1331             if (useTemme) {
1332                 evalMethod = 5;
1333             } else {
1334                 //
1335                 // Regular case where the result will not be too close to 0.5.
1336                 //
1337                 // Changeover here occurs at P ~ Q ~ 0.5
1338                 // Note that series computation of P is about x2 faster than continued fraction
1339                 // calculation of Q, so try and use the CF only when really necessary,
1340                 // especially for small x.
1341                 //
1342                 if (x - (1 / (3 * x)) < a) {
1343                     evalMethod = 2;
1344                 } else {
1345                     evalMethod = 4;
1346                     invert = !invert;
1347                 }
1348             }
1349         }
1350 
1351         switch (evalMethod) {
1352         case 0:
1353             result = finiteGammaQ(a, x);
1354             if (!normalised) {
1355                 result *= tgamma(a);
1356             }
1357             break;
1358         case 1:
1359             result = finiteHalfGammaQ(a, x);
1360             if (!normalised) {
1361                 result *= tgamma(a);
1362             }
1363             break;
1364         case 2:
1365             // Compute P:
1366             result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
1367             if (result != 0) {
1368                 //
1369                 // If we're going to be inverting the result then we can
1370                 // reduce the number of series evaluations by quite
1371                 // a few iterations if we set an initial value for the
1372                 // series sum based on what we'll end up subtracting it from
1373                 // at the end.
1374                 // Have to be careful though that this optimization doesn't
1375                 // lead to spurious numeric overflow. Note that the
1376                 // scary/expensive overflow checks below are more often
1377                 // than not bypassed in practice for "sensible" input
1378                 // values:
1379                 //
1380 
1381                 double initValue = 0;
1382                 boolean optimisedInvert = false;
1383                 if (invert) {
1384                     initValue = normalised ? 1 : tgamma(a);
1385                     if (normalised || (result >= 1) || (Double.MAX_VALUE * result > initValue)) {
1386                         initValue /= result;
1387                         if (normalised || (a < 1) || (Double.MAX_VALUE / a > initValue)) {
1388                             initValue *= -a;
1389                             optimisedInvert = true;
1390                         } else {
1391                             initValue = 0;
1392                         }
1393                     } else {
1394                         initValue = 0;
1395                     }
1396                 }
1397                 result *= lowerGammaSeries(a, x, initValue, pol) / a;
1398                 if (optimisedInvert) {
1399                     invert = false;
1400                     result = -result;
1401                 }
1402             }
1403             break;
1404         case 3:
1405             // Compute Q:
1406             invert = !invert;
1407             final double[] g = {0};
1408             result = tgammaSmallUpperPart(a, x, pol, g, invert);
1409             invert = false;
1410             if (normalised) {
1411                 // Addition to the Boost code:
1412                 if (g[0] == Double.POSITIVE_INFINITY) {
1413                     // Very small a will overflow gamma(a). Resort to logs.
1414                     // This method requires improvement as the error is very large.
1415                     // It is better than returning zero for a non-zero result.
1416                     result = Math.exp(Math.log(result) - lgamma(a));
1417                 } else {
1418                     result /= g[0];
1419                 }
1420             }
1421             break;
1422         case 4:
1423             // Compute Q:
1424             result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
1425             if (result != 0) {
1426                 result *= upperGammaFraction(a, x, pol);
1427             }
1428             break;
1429         case 5:
1430             // Call 53-bit version
1431             result = igammaTemmeLarge(a, x);
1432             if (x >= a) {
1433                 invert = !invert;
1434             }
1435             break;
1436         case 6:
1437             // x is so small that P is necessarily very small too, use
1438             // http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
1439             if (normalised) {
1440                 // If tgamma overflows then result = 0
1441                 result = Math.pow(x, a) / tgamma(a + 1);
1442             } else {
1443                 result = Math.pow(x, a) / a;
1444             }
1445             result *= 1 - a * x / (a + 1);
1446             break;
1447         case 7:
1448         default:
1449             // x is large,
1450             // Compute Q:
1451             result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
1452             result /= x;
1453             if (result != 0) {
1454                 result *= incompleteTgammaLargeX(a, x, pol);
1455             }
1456             break;
1457         }
1458 
1459         if (normalised && (result > 1)) {
1460             result = 1;
1461         }
1462         if (invert) {
1463             final double gam = normalised ? 1 : tgamma(a);
1464             result = gam - result;
1465         }
1466 
1467         return result;
1468     }
1469 
1470     /**
1471      * Upper gamma fraction.
1472      * Multiply result by z^a * e^-z to get the full
1473      * upper incomplete integral.  Divide by tgamma(z)
1474      * to normalise.
1475      *
1476      * @param a Argument a
1477      * @param z Argument z
1478      * @param pol Function evaluation policy
1479      * @return upper gamma fraction
1480      */
1481     // This is package-private for testing
1482     static double upperGammaFraction(double a, double z, Policy pol) {
1483         final double eps = pol.getEps();
1484         final int maxIterations = pol.getMaxIterations();
1485 
1486         // This is computing:
1487         //              1
1488         // ------------------------------
1489         // b0 + a1 / (b1 +     a2       )
1490         //                 -------------
1491         //                 b2 +    a3
1492         //                      --------
1493         //                      b3 + ...
1494         //
1495         // b0 = z - a + 1
1496         // a1 = a - 1
1497         //
1498         // It can be done several ways with variations in accuracy.
1499         // The current implementation has the best accuracy and matches the Boost code.
1500 
1501         final double zma1 = z - a + 1;
1502 
1503         final Supplier<Coefficient> gen = new Supplier<Coefficient>() {
1504             /** Iteration. */
1505             private int k;
1506 
1507             @Override
1508             public Coefficient get() {
1509                 ++k;
1510                 return Coefficient.of(k * (a - k), zma1 + 2.0 * k);
1511             }
1512         };
1513 
1514         return 1 / GeneralizedContinuedFraction.value(zma1, gen, eps, maxIterations);
1515     }
1516 
1517     /**
1518      * Upper gamma fraction for integer a.
1519      * Called when {@code a < 30} and {@code -x > LOG_MIN_VALUE}.
1520      *
1521      * @param a Argument a (assumed to be small)
1522      * @param x Argument x
1523      * @return upper gamma fraction
1524      */
1525     private static double finiteGammaQ(double a, double x) {
1526         //
1527         // Calculates normalised Q when a is an integer:
1528         //
1529 
1530         // Update:
1531         // Assume -x > log min value and no underflow to zero.
1532 
1533         double sum = Math.exp(-x);
1534         double term = sum;
1535         for (int n = 1; n < a; ++n) {
1536             term /= n;
1537             term *= x;
1538             sum += term;
1539         }
1540         return sum;
1541     }
1542 
1543     /**
1544      * Upper gamma fraction for half integer a.
1545      * Called when {@code a < 30} and {@code -x > LOG_MIN_VALUE}.
1546      *
1547      * @param a Argument a (assumed to be small)
1548      * @param x Argument x
1549      * @return upper gamma fraction
1550      */
1551     private static double finiteHalfGammaQ(double a, double x) {
1552         //
1553         // Calculates normalised Q when a is a half-integer:
1554         //
1555 
1556         // Update:
1557         // Assume -x > log min value:
1558         // erfc(sqrt(708)) = erfc(26.6) => erfc has a non-zero value
1559 
1560         double e = BoostErf.erfc(Math.sqrt(x));
1561         if (a > 1) {
1562             double term = Math.exp(-x) / Math.sqrt(Math.PI * x);
1563             term *= x;
1564             term /= 0.5;
1565             double sum = term;
1566             for (int n = 2; n < a; ++n) {
1567                 term /= n - 0.5;
1568                 term *= x;
1569                 sum += term;
1570             }
1571             e += sum;
1572         }
1573         return e;
1574     }
1575 
1576     /**
1577      * Lower gamma series.
1578      * Multiply result by ((z^a) * (e^-z) / a) to get the full
1579      * lower incomplete integral. Then divide by tgamma(a)
1580      * to get the normalised value.
1581      *
1582      * @param a Argument a
1583      * @param z Argument z
1584      * @param initValue Initial value
1585      * @param pol Function evaluation policy
1586      * @return lower gamma series
1587      */
1588     // This is package-private for testing
1589     static double lowerGammaSeries(double a, double z, double initValue, Policy pol) {
1590         final double eps = pol.getEps();
1591         final int maxIterations = pol.getMaxIterations();
1592 
1593         // Lower gamma series representation.
1594         final DoubleSupplier gen = new DoubleSupplier() {
1595             /** Next result. */
1596             private double result = 1;
1597             /** Iteration. */
1598             private int n;
1599 
1600             @Override
1601             public double getAsDouble() {
1602                 final double r = result;
1603                 n++;
1604                 result *= z / (a + n);
1605                 return r;
1606             }
1607         };
1608 
1609         return BoostTools.kahanSumSeries(gen, eps, maxIterations, initValue);
1610     }
1611 
1612     /**
1613      * Upper gamma fraction for very small a.
1614      *
1615      * @param a Argument a (assumed to be small)
1616      * @param x Argument x
1617      * @param pol Function evaluation policy
1618      * @param pgam set to value of gamma(a) on output
1619      * @param invert true to invert the result
1620      * @return upper gamma fraction
1621      */
1622     private static double tgammaSmallUpperPart(double a, double x, Policy pol, double[] pgam, boolean invert) {
1623         //
1624         // Compute the full upper fraction (Q) when a is very small:
1625         //
1626         double result;
1627         result = tgamma1pm1(a);
1628 
1629         // Note: Replacing this with tgamma(a) does not reduce error on current test data.
1630 
1631         // gamma(1+z) = z gamma(z)
1632         // pgam[0] == gamma(a)
1633         pgam[0] = (result + 1) / a;
1634 
1635         double p = BoostMath.powm1(x, a);
1636         result -= p;
1637         result /= a;
1638         // Removed subtraction of 10 from this value
1639         final int maxIter = pol.getMaxIterations();
1640         p += 1;
1641         final double initValue = invert ? pgam[0] : 0;
1642 
1643         // Series representation for upper fraction when z is small.
1644         final DoubleSupplier gen = new DoubleSupplier() {
1645             /** Result term. */
1646             private double result = -x;
1647             /** Argument x (this is negated on purpose). */
1648             private final double z = -x;
1649             /** Iteration. */
1650             private int n = 1;
1651 
1652             @Override
1653             public double getAsDouble() {
1654                 final double r = result / (a + n);
1655                 n++;
1656                 result = result * z / n;
1657                 return r;
1658             }
1659         };
1660 
1661         result = -p * BoostTools.kahanSumSeries(gen, pol.getEps(), maxIter, (initValue - result) / p);
1662         if (invert) {
1663             result = -result;
1664         }
1665         return result;
1666     }
1667 
1668     /**
1669      * Calculate power term prefix (z^a)(e^-z) used in the non-normalised
1670      * incomplete gammas.
1671      *
1672      * @param a Argument a
1673      * @param z Argument z
1674      * @return incomplete gamma prefix
1675      */
1676     static double fullIgammaPrefix(double a, double z) {
1677         if (z > Double.MAX_VALUE) {
1678             return 0;
1679         }
1680         final double alz = a * Math.log(z);
1681         final double prefix;
1682 
1683         if (z >= 1) {
1684             if ((alz < LOG_MAX_VALUE) && (-z > LOG_MIN_VALUE)) {
1685                 prefix = Math.pow(z, a) * Math.exp(-z);
1686             } else if (a >= 1) {
1687                 prefix = Math.pow(z / Math.exp(z / a), a);
1688             } else {
1689                 prefix = Math.exp(alz - z);
1690             }
1691         } else {
1692             if (alz > LOG_MIN_VALUE) {
1693                 prefix = Math.pow(z, a) * Math.exp(-z);
1694             } else {
1695                 // Updated to remove unreachable final branch using Math.exp(alz - z).
1696                 // This branch requires (z / a < LOG_MAX_VALUE) to avoid overflow in exp.
1697                 // At this point:
1698                 // 1. log(z) is negative;
1699                 // 2. a * log(z) <= -708 requires a > -708 / log(z).
1700                 // For any z < 1: -708 / log(z) is > z. Thus a is > z and
1701                 // z / a < LOG_MAX_VALUE is always true.
1702                 prefix = Math.pow(z / Math.exp(z / a), a);
1703             }
1704         }
1705         // Removed overflow check. Return infinity if it occurs.
1706         return prefix;
1707     }
1708 
1709     /**
1710      * Compute (z^a)(e^-z)/tgamma(a).
1711      * <p>Most of the error occurs in this function.
1712      *
1713      * @param a Argument a
1714      * @param z Argument z
1715      * @return regularized gamma prefix
1716      */
1717     // This is package-private for testing
1718     static double regularisedGammaPrefix(double a, double z) {
1719         if (z >= Double.MAX_VALUE) {
1720             return 0;
1721         }
1722 
1723         // Update this condition from: a < 1
1724         if (a <= 1) {
1725             //
1726             // We have to treat a < 1 as a special case because our Lanczos
1727             // approximations are optimised against the factorials with a > 1,
1728             // and for high precision types especially (128-bit reals for example)
1729             // very small values of a can give rather erroneous results for gamma
1730             // unless we do this:
1731             //
1732             // Boost todo: is this still required? Lanczos approx should be better now?
1733             //
1734 
1735             // Update this condition from: z <= LOG_MIN_VALUE
1736             // Require exp(-z) to not underflow:
1737             // -z > log(min_value)
1738             if (-z <= LOG_MIN_VALUE) {
1739                 // Oh dear, have to use logs, should be free of cancellation errors though:
1740                 return Math.exp(a * Math.log(z) - z - lgamma(a));
1741             }
1742             // direct calculation, no danger of overflow as gamma(a) < 1/a
1743             // for small a.
1744             return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
1745         }
1746 
1747         // Update to the Boost code.
1748         // Use some of the logic from fullIgammaPrefix(a, z) to use the direct
1749         // computation if it is valid. Assuming pow and exp are accurate to 1 ULP it
1750         // puts most of the error in evaluation of tgamma(a). This is accurate
1751         // enough that this reduces max error on the current test data.
1752         //
1753         // Overflow cases fall-through to the Lanczos approximation that incorporates
1754         // the pow and exp terms used in the tgamma(a) computation with the terms
1755         // z^a and e^-z into a single evaluation of pow and exp. See equation 15:
1756         // https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html
1757         if (a <= MAX_GAMMA_Z) {
1758             final double alz1 = a * Math.log(z);
1759             if (z >= 1) {
1760                 if ((alz1 < LOG_MAX_VALUE) && (-z > LOG_MIN_VALUE)) {
1761                     return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
1762                 }
1763             } else if (alz1 > LOG_MIN_VALUE) {
1764                 return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
1765             }
1766         }
1767 
1768         //
1769         // For smallish a and x combining the power terms with the Lanczos approximation
1770         // gives the greatest accuracy
1771         //
1772 
1773         final double agh = a + Lanczos.GMH;
1774         double prefix;
1775 
1776         final double factor = Math.sqrt(agh / Math.E) / Lanczos.lanczosSumExpGScaled(a);
1777 
1778         // Update to the Boost code.
1779         // Lower threshold for large a from 150 to 128 and compute d on demand.
1780         // See NUMBERS-179.
1781         if (a > 128) {
1782             final double d = ((z - a) - Lanczos.GMH) / agh;
1783             if (Math.abs(d * d * a) <= 100) {
1784                 // special case for large a and a ~ z.
1785                 // When a and x are large, we end up with a very large exponent with a base near one:
1786                 // this will not be computed accurately via the pow function, and taking logs simply
1787                 // leads to cancellation errors.
1788                 prefix = a * SpecialMath.log1pmx(d) + z * -Lanczos.GMH / agh;
1789                 prefix = Math.exp(prefix);
1790                 return prefix * factor;
1791             }
1792         }
1793 
1794         //
1795         // general case.
1796         // direct computation is most accurate, but use various fallbacks
1797         // for different parts of the problem domain:
1798         //
1799 
1800         final double alz = a * Math.log(z / agh);
1801         final double amz = a - z;
1802         if ((Math.min(alz, amz) <= LOG_MIN_VALUE) || (Math.max(alz, amz) >= LOG_MAX_VALUE)) {
1803             final double amza = amz / a;
1804             if ((Math.min(alz, amz) / 2 > LOG_MIN_VALUE) && (Math.max(alz, amz) / 2 < LOG_MAX_VALUE)) {
1805                 // compute square root of the result and then square it:
1806                 final double sq = Math.pow(z / agh, a / 2) * Math.exp(amz / 2);
1807                 prefix = sq * sq;
1808             } else if ((Math.min(alz, amz) / 4 > LOG_MIN_VALUE) &&
1809                     (Math.max(alz, amz) / 4 < LOG_MAX_VALUE) && (z > a)) {
1810                 // compute the 4th root of the result then square it twice:
1811                 final double sq = Math.pow(z / agh, a / 4) * Math.exp(amz / 4);
1812                 prefix = sq * sq;
1813                 prefix *= prefix;
1814             } else if ((amza > LOG_MIN_VALUE) && (amza < LOG_MAX_VALUE)) {
1815                 prefix = Math.pow((z * Math.exp(amza)) / agh, a);
1816             } else {
1817                 prefix = Math.exp(alz + amz);
1818             }
1819         } else {
1820             prefix = Math.pow(z / agh, a) * Math.exp(amz);
1821         }
1822         prefix *= factor;
1823         return prefix;
1824     }
1825 
1826     /**
1827      * Implements the asymptotic expansions of the incomplete
1828      * gamma functions P(a, x) and Q(a, x), used when a is large and
1829      * x ~ a.
1830      *
1831      * <p>The primary reference is:
1832      * <pre>
1833      * "The Asymptotic Expansion of the Incomplete Gamma Functions"
1834      * N. M. Temme.
1835      * Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
1836      * </pre>
1837      *
1838      * <p>A different way of evaluating these expansions,
1839      * plus a lot of very useful background information is in:
1840      * <pre>
1841      * "A Set of Algorithms For the Incomplete Gamma Functions."
1842      * N. M. Temme.
1843      * Probability in the Engineering and Informational Sciences,
1844      * 8, 1994, 291.
1845      * </pre>
1846      *
1847      * <p>An alternative implementation is in:
1848      * <pre>
1849      * "Computation of the Incomplete Gamma Function Ratios and their Inverse."
1850      * A. R. Didonato and A. H. Morris.
1851      * ACM TOMS, Vol 12, No 4, Dec 1986, p377.
1852      * </pre>
1853      *
1854      * <p>This is a port of the function accurate for 53-bit mantissas
1855      * (IEEE double precision or 10^-17). To understand the code, refer to Didonato
1856      * and Morris, from Eq 17 and 18 onwards.
1857      *
1858      * <p>The coefficients used here are not taken from Didonato and Morris:
1859      * the domain over which these expansions are used is slightly different
1860      * to theirs, and their constants are not quite accurate enough for
1861      * 128-bit long doubles.  Instead the coefficients were calculated
1862      * using the methods described by Temme p762 from Eq 3.8 onwards.
1863      * The values obtained agree with those obtained by Didonato and Morris
1864      * (at least to the first 30 digits that they provide).
1865      * At double precision the degrees of polynomial required for full
1866      * machine precision are close to those recommended to Didonato and Morris,
1867      * but of course many more terms are needed for larger types.
1868      *
1869      * <p>Adapted from {@code boost/math/special_functions/detail/igamma_large.hpp}.
1870      *
1871      * @param a the a
1872      * @param x the x
1873      * @return the double
1874      */
1875     // This is package-private for testing
1876     static double igammaTemmeLarge(double a, double x) {
1877         final double sigma = (x - a) / a;
1878         final double phi = -SpecialMath.log1pmx(sigma);
1879         final double y = a * phi;
1880         double z = Math.sqrt(2 * phi);
1881         if (x < a) {
1882             z = -z;
1883         }
1884 
1885         // The following polynomials are evaluated with a loop
1886         // with Horner's method. Variations exist using
1887         // a second order Horner's method with an unrolled loop.
1888         // These are chosen in Boost based on the C++ compiler.
1889         // For example:
1890         // boost/math/tools/detail/polynomial_horner1_15.hpp
1891         // boost/math/tools/detail/polynomial_horner2_15.hpp
1892         // boost/math/tools/detail/polynomial_horner3_15.hpp
1893 
1894         final double[] workspace = new double[10];
1895 
1896         final double[] C0 = {
1897             -0.33333333333333333,
1898             0.083333333333333333,
1899             -0.014814814814814815,
1900             0.0011574074074074074,
1901             0.0003527336860670194,
1902             -0.00017875514403292181,
1903             0.39192631785224378e-4,
1904             -0.21854485106799922e-5,
1905             -0.185406221071516e-5,
1906             0.8296711340953086e-6,
1907             -0.17665952736826079e-6,
1908             0.67078535434014986e-8,
1909             0.10261809784240308e-7,
1910             -0.43820360184533532e-8,
1911             0.91476995822367902e-9,
1912         };
1913         workspace[0] = BoostTools.evaluatePolynomial(C0, z);
1914 
1915         final double[] C1 = {
1916             -0.0018518518518518519,
1917             -0.0034722222222222222,
1918             0.0026455026455026455,
1919             -0.00099022633744855967,
1920             0.00020576131687242798,
1921             -0.40187757201646091e-6,
1922             -0.18098550334489978e-4,
1923             0.76491609160811101e-5,
1924             -0.16120900894563446e-5,
1925             0.46471278028074343e-8,
1926             0.1378633446915721e-6,
1927             -0.5752545603517705e-7,
1928             0.11951628599778147e-7,
1929         };
1930         workspace[1] = BoostTools.evaluatePolynomial(C1, z);
1931 
1932         final double[] C2 = {
1933             0.0041335978835978836,
1934             -0.0026813271604938272,
1935             0.00077160493827160494,
1936             0.20093878600823045e-5,
1937             -0.00010736653226365161,
1938             0.52923448829120125e-4,
1939             -0.12760635188618728e-4,
1940             0.34235787340961381e-7,
1941             0.13721957309062933e-5,
1942             -0.6298992138380055e-6,
1943             0.14280614206064242e-6,
1944         };
1945         workspace[2] = BoostTools.evaluatePolynomial(C2, z);
1946 
1947         final double[] C3 = {
1948             0.00064943415637860082,
1949             0.00022947209362139918,
1950             -0.00046918949439525571,
1951             0.00026772063206283885,
1952             -0.75618016718839764e-4,
1953             -0.23965051138672967e-6,
1954             0.11082654115347302e-4,
1955             -0.56749528269915966e-5,
1956             0.14230900732435884e-5,
1957         };
1958         workspace[3] = BoostTools.evaluatePolynomial(C3, z);
1959 
1960         final double[] C4 = {
1961             -0.0008618882909167117,
1962             0.00078403922172006663,
1963             -0.00029907248030319018,
1964             -0.14638452578843418e-5,
1965             0.66414982154651222e-4,
1966             -0.39683650471794347e-4,
1967             0.11375726970678419e-4,
1968         };
1969         workspace[4] = BoostTools.evaluatePolynomial(C4, z);
1970 
1971         final double[] C5 = {
1972             -0.00033679855336635815,
1973             -0.69728137583658578e-4,
1974             0.00027727532449593921,
1975             -0.00019932570516188848,
1976             0.67977804779372078e-4,
1977             0.1419062920643967e-6,
1978             -0.13594048189768693e-4,
1979             0.80184702563342015e-5,
1980             -0.22914811765080952e-5,
1981         };
1982         workspace[5] = BoostTools.evaluatePolynomial(C5, z);
1983 
1984         final double[] C6 = {
1985             0.00053130793646399222,
1986             -0.00059216643735369388,
1987             0.00027087820967180448,
1988             0.79023532326603279e-6,
1989             -0.81539693675619688e-4,
1990             0.56116827531062497e-4,
1991             -0.18329116582843376e-4,
1992         };
1993         workspace[6] = BoostTools.evaluatePolynomial(C6, z);
1994 
1995         final double[] C7 = {
1996             0.00034436760689237767,
1997             0.51717909082605922e-4,
1998             -0.00033493161081142236,
1999             0.0002812695154763237,
2000             -0.00010976582244684731,
2001         };
2002         workspace[7] = BoostTools.evaluatePolynomial(C7, z);
2003 
2004         final double[] C8 = {
2005             -0.00065262391859530942,
2006             0.00083949872067208728,
2007             -0.00043829709854172101,
2008         };
2009         workspace[8] = BoostTools.evaluatePolynomial(C8, z);
2010         workspace[9] = -0.00059676129019274625;
2011 
2012         double result = BoostTools.evaluatePolynomial(workspace, 1 / a);
2013         result *= Math.exp(-y) / Math.sqrt(2 * Math.PI * a);
2014         if (x < a) {
2015             result = -result;
2016         }
2017 
2018         result += BoostErf.erfc(Math.sqrt(y)) / 2;
2019 
2020         return result;
2021     }
2022 
2023     /**
2024      * Incomplete tgamma for large X.
2025      *
2026      * <p>This summation is a source of error as the series starts at 1 and descends to zero.
2027      * It can have thousands of iterations when a and z are large and close in value.
2028      *
2029      * @param a Argument a
2030      * @param x Argument x
2031      * @param pol Function evaluation policy
2032      * @return incomplete tgamma
2033      */
2034     // This is package-private for testing
2035     static double incompleteTgammaLargeX(double a, double x, Policy pol) {
2036         final double eps = pol.getEps();
2037         final int maxIterations = pol.getMaxIterations();
2038 
2039         // Asymptotic approximation for large argument, see: https://dlmf.nist.gov/8.11#E2.
2040         final DoubleSupplier gen = new DoubleSupplier() {
2041             /** Result term. */
2042             private double term = 1;
2043             /** Iteration. */
2044             private int n;
2045 
2046             @Override
2047             public double getAsDouble() {
2048                 final double result = term;
2049                 n++;
2050                 term *= (a - n) / x;
2051                 return result;
2052             }
2053         };
2054 
2055         return BoostTools.kahanSumSeries(gen, eps, maxIterations);
2056     }
2057 
2058     /**
2059      * Return true if the array is not null and has non-zero length.
2060      *
2061      * @param array Array
2062      * @return true if a non-zero length array
2063      */
2064     private static boolean nonZeroLength(int[] array) {
2065         return array != null && array.length != 0;
2066     }
2067 
2068     /**
2069      * Ratio of gamma functions.
2070      *
2071      * <p>\[ tgamma_ratio(z, delta) = \frac{\Gamma(z)}{\Gamma(z + delta)} \]
2072      *
2073      * <p>Adapted from {@code tgamma_delta_ratio_imp}. The use of
2074      * {@code max_factorial<double>::value == 170} has been replaced with
2075      * {@code MAX_GAMMA_Z == 171}. This threshold is used when it is possible
2076      * to call the gamma function without overflow.
2077      *
2078      * @param z Argument z
2079      * @param delta The difference
2080      * @return gamma ratio
2081      */
2082     static double tgammaDeltaRatio(double z, double delta) {
2083         final double zDelta = z + delta;
2084         if (Double.isNaN(zDelta)) {
2085             // One or both arguments are NaN
2086             return Double.NaN;
2087         }
2088         if (z <= 0 || zDelta <= 0) {
2089             // This isn't very sophisticated, or accurate, but it does work:
2090             return tgamma(z) / tgamma(zDelta);
2091         }
2092 
2093         // Note: Directly calling tgamma(z) / tgamma(z + delta) if possible
2094         // without overflow is not more accurate
2095 
2096         if (Math.rint(delta) == delta) {
2097             if (delta == 0) {
2098                 return 1;
2099             }
2100             //
2101             // If both z and delta are integers, see if we can just use table lookup
2102             // of the factorials to get the result:
2103             //
2104             if (Math.rint(z) == z &&
2105                 z <= MAX_GAMMA_Z && zDelta <= MAX_GAMMA_Z) {
2106                 return FACTORIAL[(int) z - 1] / FACTORIAL[(int) zDelta - 1];
2107             }
2108             if (Math.abs(delta) < 20) {
2109                 //
2110                 // delta is a small integer, we can use a finite product:
2111                 //
2112                 if (delta < 0) {
2113                     z -= 1;
2114                     double result = z;
2115                     for (int d = (int) (delta + 1); d != 0; d++) {
2116                         z -= 1;
2117                         result *= z;
2118                     }
2119                     return result;
2120                 }
2121                 double result = 1 / z;
2122                 for (int d = (int) (delta - 1); d != 0; d--) {
2123                     z += 1;
2124                     result /= z;
2125                 }
2126                 return result;
2127             }
2128         }
2129         return tgammaDeltaRatioImpLanczos(z, delta);
2130     }
2131 
2132     /**
2133      * Ratio of gamma functions using Lanczos support.
2134      *
2135      * <p>\[ tgamma_delta_ratio(z, delta) = \frac{\Gamma(z)}{\Gamma(z + delta)}
2136      *
2137      * <p>Adapted from {@code tgamma_delta_ratio_imp_lanczos}. The use of
2138      * {@code max_factorial<double>::value == 170} has been replaced with
2139      * {@code MAX_GAMMA_Z == 171}. This threshold is used when it is possible
2140      * to use the precomputed factorial table.
2141      *
2142      * @param z Argument z
2143      * @param delta The difference
2144      * @return gamma ratio
2145      */
2146     private static double tgammaDeltaRatioImpLanczos(double z, double delta) {
2147         if (z < EPSILON) {
2148             //
2149             // We get spurious numeric overflow unless we're very careful, this
2150             // can occur either inside Lanczos::lanczos_sum(z) or in the
2151             // final combination of terms, to avoid this, split the product up
2152             // into 2 (or 3) parts:
2153             //
2154             // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta
2155             // z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial
2156             //
2157             if (MAX_GAMMA_Z < delta) {
2158                 double ratio = tgammaDeltaRatioImpLanczos(delta, MAX_GAMMA_Z - delta);
2159                 ratio *= z;
2160                 ratio *= FACTORIAL[MAX_FACTORIAL];
2161                 return 1 / ratio;
2162             }
2163             return 1 / (z * tgamma(z + delta));
2164         }
2165         final double zgh = z + Lanczos.GMH;
2166         double result;
2167         if (z + delta == z) {
2168             // Here:
2169             // lanczosSum(z) / lanczosSum(z + delta) == 1
2170 
2171             // Update to the Boost code to remove unreachable code:
2172             // Given z + delta == z then |delta / z| < EPSILON; and z < zgh
2173             // assume (abs(delta / zgh) < EPSILON) and remove unreachable
2174             // else branch which sets result = 1
2175 
2176             // We have:
2177             // result = exp((0.5 - z) * log1p(delta / zgh));
2178             // 0.5 - z == -z
2179             // log1p(delta / zgh) = delta / zgh = delta / z
2180             // multiplying we get -delta.
2181 
2182             // Note:
2183             // This can be different from
2184             // exp((0.5 - z) * log1p(delta / zgh)) when z is small.
2185             // In this case the result is exp(small) and the result
2186             // is within 1 ULP of 1.0. This is left as the original
2187             // Boost method using exp(-delta).
2188 
2189             result = Math.exp(-delta);
2190         } else {
2191             if (Math.abs(delta) < 10) {
2192                 result = Math.exp((0.5 - z) * Math.log1p(delta / zgh));
2193             } else {
2194                 result = Math.pow(zgh / (zgh + delta), z - 0.5);
2195             }
2196             // Split the calculation up to avoid spurious overflow:
2197             result *= Lanczos.lanczosSum(z) / Lanczos.lanczosSum(z + delta);
2198         }
2199         result *= Math.pow(Math.E / (zgh + delta), delta);
2200         return result;
2201     }
2202 
2203     /**
2204      * Ratio of gamma functions.
2205      *
2206      * <p>\[ tgamma_ratio(x, y) = \frac{\Gamma(x)}{\Gamma(y)}
2207      *
2208      * <p>Adapted from {@code tgamma_ratio_imp}. The use of
2209      * {@code max_factorial<double>::value == 170} has been replaced with
2210      * {@code MAX_GAMMA_Z == 171}. This threshold is used when it is possible
2211      * to call the gamma function without overflow.
2212      *
2213      * @param x Argument x (must be positive finite)
2214      * @param y Argument y (must be positive finite)
2215      * @return gamma ratio (or nan)
2216      */
2217     static double tgammaRatio(double x, double y) {
2218         if (x <= 0 || !Double.isFinite(x) || y <= 0 || !Double.isFinite(y)) {
2219             return Double.NaN;
2220         }
2221         if (x <= Double.MIN_NORMAL) {
2222             // Special case for denorms...Ugh.
2223             return TWO_POW_53 * tgammaRatio(x * TWO_POW_53, y);
2224         }
2225 
2226         if (x <= MAX_GAMMA_Z && y <= MAX_GAMMA_Z) {
2227             // Rather than subtracting values, lets just call the gamma functions directly:
2228             return tgamma(x) / tgamma(y);
2229         }
2230         double prefix = 1;
2231         if (x < 1) {
2232             if (y < 2 * MAX_GAMMA_Z) {
2233                 // We need to sidestep on x as well, otherwise we'll underflow
2234                 // before we get to factor in the prefix term:
2235                 prefix /= x;
2236                 x += 1;
2237                 while (y >= MAX_GAMMA_Z) {
2238                     y -= 1;
2239                     prefix /= y;
2240                 }
2241                 return prefix * tgamma(x) / tgamma(y);
2242             }
2243             //
2244             // result is almost certainly going to underflow to zero, try logs just in case:
2245             //
2246             return Math.exp(lgamma(x) - lgamma(y));
2247         }
2248         if (y < 1) {
2249             if (x < 2 * MAX_GAMMA_Z) {
2250                 // We need to sidestep on y as well, otherwise we'll overflow
2251                 // before we get to factor in the prefix term:
2252                 prefix *= y;
2253                 y += 1;
2254                 while (x >= MAX_GAMMA_Z) {
2255                     x -= 1;
2256                     prefix *= x;
2257                 }
2258                 return prefix * tgamma(x) / tgamma(y);
2259             }
2260             //
2261             // Result will almost certainly overflow, try logs just in case:
2262             //
2263             return Math.exp(lgamma(x) - lgamma(y));
2264         }
2265         //
2266         // Regular case, x and y both large and similar in magnitude:
2267         //
2268         return tgammaDeltaRatio(x, y - x);
2269     }
2270 }