BoostGamma.java

  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. //  Copyright John Maddock 2006-7, 2013-20.
  18. //  Copyright Paul A. Bristow 2007, 2013-14.
  19. //  Copyright Nikhar Agrawal 2013-14
  20. //  Copyright Christopher Kormanyos 2013-14, 2020
  21. //  Use, modification and distribution are subject to the
  22. //  Boost Software License, Version 1.0. (See accompanying file
  23. //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

  24. package org.apache.commons.numbers.gamma;

  25. import java.util.function.DoubleSupplier;
  26. import java.util.function.Supplier;
  27. import org.apache.commons.numbers.core.Sum;
  28. import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction;
  29. import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction.Coefficient;

  30. /**
  31.  * Implementation of the
  32.  * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">Regularized Gamma functions</a> and
  33.  * <a href="https://mathworld.wolfram.com/IncompleteGammaFunction.html">Incomplete Gamma functions</a>.
  34.  *
  35.  * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
  36.  * {@code c++} implementation {@code <boost/math/special_functions/gamma.hpp>}.
  37.  * All work is copyright to the original authors and subject to the Boost Software License.
  38.  *
  39.  * @see
  40.  * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma.html">
  41.  * Boost C++ Gamma functions</a>
  42.  */
  43. final class BoostGamma {
  44.     //
  45.     // Code ported from Boost 1.77.0
  46.     //
  47.     // boost/math/special_functions/gamma.hpp
  48.     // boost/math/special_functions/detail/igamma_large.hpp
  49.     // boost/math/special_functions/lanczos.hpp
  50.     //
  51.     // Original code comments are preserved.
  52.     //
  53.     // Changes to the Boost implementation:
  54.     // - Update method names to replace underscores with camel case
  55.     // - Explicitly inline the polynomial function evaluation
  56.     //   using Horner's method (https://en.wikipedia.org/wiki/Horner%27s_method)
  57.     // - Remove checks for under/overflow. In this implementation no error is raised
  58.     //   for overflow (infinity is returned) or underflow (sub-normal or zero is returned).
  59.     //   This follows the conventions in java.lang.Math for the same conditions.
  60.     // - Removed the pointer p_derivative in the gammaIncompleteImp. This is used
  61.     //   in the Boost code for the gamma_(p|q)_inv functions for a derivative
  62.     //   based inverse function. This is currently not supported.
  63.     // - Added extended precision arithmetic for some series summations or other computations.
  64.     //   The Boost default policy is to evaluate in long double for a double result. Extended
  65.     //   precision is not possible for the entire computation but has been used where
  66.     //   possible for some terms to reduce errors. Error reduction verified on the test data.
  67.     // - Altered the tgamma(x) function to use the double-precision NSWC Library of Mathematics
  68.     //   Subroutines when the error is lower. This is for the non Lanczos code.
  69.     // - Altered the condition used for the asymptotic approximation method to avoid
  70.     //   loss of precision in the series summation when a ~ z.
  71.     // - Altered series generators to use integer counters added to the double term
  72.     //   replacing directly incrementing a double term. When the term is large it cannot
  73.     //   be incremented: 1e16 + 1 == 1e16.
  74.     // - Removed unreachable code branch in tgammaDeltaRatioImpLanczos when z + delta == z.
  75.     //
  76.     // Note:
  77.     // The major source of error is in the function regularisedGammaPrefix when computing
  78.     // (z^a)(e^-z)/tgamma(a) with extreme input to the power and exponential terms.
  79.     // An extended precision pow and exp function returning a quad length result would
  80.     // be required to reduce error for these arguments. Tests using the Dfp class
  81.     // from o.a.c.math4.legacy.core have been demonstrated to effectively eliminate the
  82.     // errors from the power terms and improve accuracy on the current test data.
  83.     // In the interest of performance the Dfp class is not used in this version.

  84.     /** Default epsilon value for relative error.
  85.      * This is equal to the Boost constant {@code boost::math::tools::epsilon<double>()}. */
  86.     private static final double EPSILON = 0x1.0p-52;
  87.     /** Value for the sqrt of the epsilon for relative error.
  88.      * This is equal to the Boost constant {@code boost::math::tools::root_epsilon<double>()}. */
  89.     private static final double ROOT_EPSILON = 1.4901161193847656E-8;
  90.     /** Approximate value for ln(Double.MAX_VALUE).
  91.      * This is equal to the Boost constant {@code boost::math::tools::log_max_value<double>()}.
  92.      * No term {@code x} should be used in {@code exp(x)} if {@code x > LOG_MAX_VALUE} to avoid
  93.      * overflow. */
  94.     private static final int LOG_MAX_VALUE = 709;
  95.     /** Approximate value for ln(Double.MIN_VALUE).
  96.      * This is equal to the Boost constant {@code boost::math::tools::log_min_value<double>()}.
  97.      * No term {@code x} should be used in {@code exp(x)} if {@code x < LOG_MIN_VALUE} to avoid
  98.      * underflow to sub-normal or zero. */
  99.     private static final int LOG_MIN_VALUE = -708;
  100.     /** The largest factorial that can be represented as a double.
  101.      * This is equal to the Boost constant {@code boost::math::max_factorial<double>::value}. */
  102.     private static final int MAX_FACTORIAL = 170;
  103.     /** The largest integer value for gamma(z) that can be represented as a double. */
  104.     private static final int MAX_GAMMA_Z = MAX_FACTORIAL + 1;
  105.     /** ln(sqrt(2 pi)). Computed to 25-digits precision. */
  106.     private static final double LOG_ROOT_TWO_PI = 0.9189385332046727417803297;
  107.     /** ln(pi). Computed to 25-digits precision. */
  108.     private static final double LOG_PI = 1.144729885849400174143427;
  109.     /** Euler's constant. */
  110.     private static final double EULER = 0.5772156649015328606065120900824024310;
  111.     /** The threshold value for choosing the Lanczos approximation. */
  112.     private static final int LANCZOS_THRESHOLD = 20;
  113.     /** 2^53. */
  114.     private static final double TWO_POW_53 = 0x1.0p53;

  115.     /** All factorials that can be represented as a double. Size = 171. */
  116.     private static final double[] FACTORIAL = {
  117.         1,
  118.         1,
  119.         2,
  120.         6,
  121.         24,
  122.         120,
  123.         720,
  124.         5040,
  125.         40320,
  126.         362880.0,
  127.         3628800.0,
  128.         39916800.0,
  129.         479001600.0,
  130.         6227020800.0,
  131.         87178291200.0,
  132.         1307674368000.0,
  133.         20922789888000.0,
  134.         355687428096000.0,
  135.         6402373705728000.0,
  136.         121645100408832000.0,
  137.         0.243290200817664e19,
  138.         0.5109094217170944e20,
  139.         0.112400072777760768e22,
  140.         0.2585201673888497664e23,
  141.         0.62044840173323943936e24,
  142.         0.15511210043330985984e26,
  143.         0.403291461126605635584e27,
  144.         0.10888869450418352160768e29,
  145.         0.304888344611713860501504e30,
  146.         0.8841761993739701954543616e31,
  147.         0.26525285981219105863630848e33,
  148.         0.822283865417792281772556288e34,
  149.         0.26313083693369353016721801216e36,
  150.         0.868331761881188649551819440128e37,
  151.         0.29523279903960414084761860964352e39,
  152.         0.103331479663861449296666513375232e41,
  153.         0.3719933267899012174679994481508352e42,
  154.         0.137637530912263450463159795815809024e44,
  155.         0.5230226174666011117600072241000742912e45,
  156.         0.203978820811974433586402817399028973568e47,
  157.         0.815915283247897734345611269596115894272e48,
  158.         0.3345252661316380710817006205344075166515e50,
  159.         0.1405006117752879898543142606244511569936e52,
  160.         0.6041526306337383563735513206851399750726e53,
  161.         0.265827157478844876804362581101461589032e55,
  162.         0.1196222208654801945619631614956577150644e57,
  163.         0.5502622159812088949850305428800254892962e58,
  164.         0.2586232415111681806429643551536119799692e60,
  165.         0.1241391559253607267086228904737337503852e62,
  166.         0.6082818640342675608722521633212953768876e63,
  167.         0.3041409320171337804361260816606476884438e65,
  168.         0.1551118753287382280224243016469303211063e67,
  169.         0.8065817517094387857166063685640376697529e68,
  170.         0.427488328406002556429801375338939964969e70,
  171.         0.2308436973392413804720927426830275810833e72,
  172.         0.1269640335365827592596510084756651695958e74,
  173.         0.7109985878048634518540456474637249497365e75,
  174.         0.4052691950487721675568060190543232213498e77,
  175.         0.2350561331282878571829474910515074683829e79,
  176.         0.1386831185456898357379390197203894063459e81,
  177.         0.8320987112741390144276341183223364380754e82,
  178.         0.507580213877224798800856812176625227226e84,
  179.         0.3146997326038793752565312235495076408801e86,
  180.         0.1982608315404440064116146708361898137545e88,
  181.         0.1268869321858841641034333893351614808029e90,
  182.         0.8247650592082470666723170306785496252186e91,
  183.         0.5443449390774430640037292402478427526443e93,
  184.         0.3647111091818868528824985909660546442717e95,
  185.         0.2480035542436830599600990418569171581047e97,
  186.         0.1711224524281413113724683388812728390923e99,
  187.         0.1197857166996989179607278372168909873646e101,
  188.         0.8504785885678623175211676442399260102886e102,
  189.         0.6123445837688608686152407038527467274078e104,
  190.         0.4470115461512684340891257138125051110077e106,
  191.         0.3307885441519386412259530282212537821457e108,
  192.         0.2480914081139539809194647711659403366093e110,
  193.         0.188549470166605025498793226086114655823e112,
  194.         0.1451830920282858696340707840863082849837e114,
  195.         0.1132428117820629783145752115873204622873e116,
  196.         0.8946182130782975286851441715398316520698e117,
  197.         0.7156945704626380229481153372318653216558e119,
  198.         0.5797126020747367985879734231578109105412e121,
  199.         0.4753643337012841748421382069894049466438e123,
  200.         0.3945523969720658651189747118012061057144e125,
  201.         0.3314240134565353266999387579130131288001e127,
  202.         0.2817104114380550276949479442260611594801e129,
  203.         0.2422709538367273238176552320344125971528e131,
  204.         0.210775729837952771721360051869938959523e133,
  205.         0.1854826422573984391147968456455462843802e135,
  206.         0.1650795516090846108121691926245361930984e137,
  207.         0.1485715964481761497309522733620825737886e139,
  208.         0.1352001527678402962551665687594951421476e141,
  209.         0.1243841405464130725547532432587355307758e143,
  210.         0.1156772507081641574759205162306240436215e145,
  211.         0.1087366156656743080273652852567866010042e147,
  212.         0.103299784882390592625997020993947270954e149,
  213.         0.9916779348709496892095714015418938011582e150,
  214.         0.9619275968248211985332842594956369871234e152,
  215.         0.942689044888324774562618574305724247381e154,
  216.         0.9332621544394415268169923885626670049072e156,
  217.         0.9332621544394415268169923885626670049072e158,
  218.         0.9425947759838359420851623124482936749562e160,
  219.         0.9614466715035126609268655586972595484554e162,
  220.         0.990290071648618040754671525458177334909e164,
  221.         0.1029901674514562762384858386476504428305e167,
  222.         0.1081396758240290900504101305800329649721e169,
  223.         0.1146280563734708354534347384148349428704e171,
  224.         0.1226520203196137939351751701038733888713e173,
  225.         0.132464181945182897449989183712183259981e175,
  226.         0.1443859583202493582204882102462797533793e177,
  227.         0.1588245541522742940425370312709077287172e179,
  228.         0.1762952551090244663872161047107075788761e181,
  229.         0.1974506857221074023536820372759924883413e183,
  230.         0.2231192748659813646596607021218715118256e185,
  231.         0.2543559733472187557120132004189335234812e187,
  232.         0.2925093693493015690688151804817735520034e189,
  233.         0.339310868445189820119825609358857320324e191,
  234.         0.396993716080872089540195962949863064779e193,
  235.         0.4684525849754290656574312362808384164393e195,
  236.         0.5574585761207605881323431711741977155627e197,
  237.         0.6689502913449127057588118054090372586753e199,
  238.         0.8094298525273443739681622845449350829971e201,
  239.         0.9875044200833601362411579871448208012564e203,
  240.         0.1214630436702532967576624324188129585545e206,
  241.         0.1506141741511140879795014161993280686076e208,
  242.         0.1882677176888926099743767702491600857595e210,
  243.         0.237217324288004688567714730513941708057e212,
  244.         0.3012660018457659544809977077527059692324e214,
  245.         0.3856204823625804217356770659234636406175e216,
  246.         0.4974504222477287440390234150412680963966e218,
  247.         0.6466855489220473672507304395536485253155e220,
  248.         0.8471580690878820510984568758152795681634e222,
  249.         0.1118248651196004307449963076076169029976e225,
  250.         0.1487270706090685728908450891181304809868e227,
  251.         0.1992942746161518876737324194182948445223e229,
  252.         0.269047270731805048359538766214698040105e231,
  253.         0.3659042881952548657689727220519893345429e233,
  254.         0.5012888748274991661034926292112253883237e235,
  255.         0.6917786472619488492228198283114910358867e237,
  256.         0.9615723196941089004197195613529725398826e239,
  257.         0.1346201247571752460587607385894161555836e242,
  258.         0.1898143759076170969428526414110767793728e244,
  259.         0.2695364137888162776588507508037290267094e246,
  260.         0.3854370717180072770521565736493325081944e248,
  261.         0.5550293832739304789551054660550388118e250,
  262.         0.80479260574719919448490292577980627711e252,
  263.         0.1174997204390910823947958271638517164581e255,
  264.         0.1727245890454638911203498659308620231933e257,
  265.         0.2556323917872865588581178015776757943262e259,
  266.         0.380892263763056972698595524350736933546e261,
  267.         0.571338395644585459047893286526105400319e263,
  268.         0.8627209774233240431623188626544191544816e265,
  269.         0.1311335885683452545606724671234717114812e268,
  270.         0.2006343905095682394778288746989117185662e270,
  271.         0.308976961384735088795856467036324046592e272,
  272.         0.4789142901463393876335775239063022722176e274,
  273.         0.7471062926282894447083809372938315446595e276,
  274.         0.1172956879426414428192158071551315525115e279,
  275.         0.1853271869493734796543609753051078529682e281,
  276.         0.2946702272495038326504339507351214862195e283,
  277.         0.4714723635992061322406943211761943779512e285,
  278.         0.7590705053947218729075178570936729485014e287,
  279.         0.1229694218739449434110178928491750176572e290,
  280.         0.2004401576545302577599591653441552787813e292,
  281.         0.3287218585534296227263330311644146572013e294,
  282.         0.5423910666131588774984495014212841843822e296,
  283.         0.9003691705778437366474261723593317460744e298,
  284.         0.1503616514864999040201201707840084015944e301,
  285.         0.2526075744973198387538018869171341146786e303,
  286.         0.4269068009004705274939251888899566538069e305,
  287.         0.7257415615307998967396728211129263114717e307,
  288.     };

  289.     /**
  290.      * 53-bit precision implementation of the Lanczos approximation.
  291.      *
  292.      * <p>This implementation is in partial fraction form with the leading constant
  293.      * of \( \sqrt{2\pi} \) absorbed into the sum.
  294.      *
  295.      * <p>It is related to the Gamma function by the following equation
  296.      * \[
  297.      * \Gamma(z) = \frac{(z + g - 0.5)^{z - 0.5}}{e^{z + g - 0.5}} \mathrm{lanczos}(z)
  298.      * \]
  299.      * where \( g \) is the Lanczos constant.
  300.      *
  301.      * <h2>Warning</h2>
  302.      *
  303.      * <p>This is not a substitute for {@link LanczosApproximation}. The approximation is
  304.      * written in partial fraction form with the leading constants absorbed by the
  305.      * coefficients in the sum.
  306.      *
  307.      * @see <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/lanczos.html">
  308.      * Boost Lanczos Approximation</a>
  309.      */
  310.     static final class Lanczos {
  311.         // Optimal values for G for each N are taken from
  312.         // http://web.mala.bc.ca/pughg/phdThesis/phdThesis.pdf,
  313.         // as are the theoretical error bounds.
  314.         //
  315.         // Constants calculated using the method described by Godfrey
  316.         // http://my.fit.edu/~gabdo/gamma.txt and elaborated by Toth at
  317.         // http://www.rskey.org/gamma.htm using NTL::RR at 1000 bit precision.

  318.         //
  319.         // Lanczos Coefficients for N=13 G=6.024680040776729583740234375
  320.         // Max experimental error (with arbitrary precision arithmetic) 1.196214e-17
  321.         // Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
  322.         //

  323.         /**
  324.          * Lanczos constant G.
  325.          */
  326.         static final double G = 6.024680040776729583740234375;

  327.         /**
  328.          * Lanczos constant G - half.
  329.          *
  330.          * <p>Note: The form {@code (g - 0.5)} is used when computing the gamma function.
  331.          */
  332.         static final double GMH = 5.524680040776729583740234375;

  333.         /** Common denominator used for the rational evaluation. */
  334.         private static final int[] DENOM = {
  335.             0,
  336.             39916800,
  337.             120543840,
  338.             150917976,
  339.             105258076,
  340.             45995730,
  341.             13339535,
  342.             2637558,
  343.             357423,
  344.             32670,
  345.             1925,
  346.             66,
  347.             1
  348.         };

  349.         /** Private constructor. */
  350.         private Lanczos() {
  351.             // intentionally empty.
  352.         }

  353.         /**
  354.          * Computes the Lanczos approximation.
  355.          *
  356.          * @param z Argument.
  357.          * @return the Lanczos approximation.
  358.          */
  359.         static double lanczosSum(double z) {
  360.             final double[] num = {
  361.                 23531376880.41075968857200767445163675473,
  362.                 42919803642.64909876895789904700198885093,
  363.                 35711959237.35566804944018545154716670596,
  364.                 17921034426.03720969991975575445893111267,
  365.                 6039542586.35202800506429164430729792107,
  366.                 1439720407.311721673663223072794912393972,
  367.                 248874557.8620541565114603864132294232163,
  368.                 31426415.58540019438061423162831820536287,
  369.                 2876370.628935372441225409051620849613599,
  370.                 186056.2653952234950402949897160456992822,
  371.                 8071.672002365816210638002902272250613822,
  372.                 210.8242777515793458725097339207133627117,
  373.                 2.506628274631000270164908177133837338626
  374.             };
  375.             return evaluateRational(num, DENOM, z);
  376.         }

  377.         /**
  378.          * Computes the Lanczos approximation scaled by {@code exp(g)}.
  379.          *
  380.          * @param z Argument.
  381.          * @return the scaled Lanczos approximation.
  382.          */
  383.         static double lanczosSumExpGScaled(double z) {
  384.             // As above with numerator divided by exp(g) = 413.509...
  385.             final double[] num = {
  386.                 56906521.91347156388090791033559122686859,
  387.                 103794043.1163445451906271053616070238554,
  388.                 86363131.28813859145546927288977868422342,
  389.                 43338889.32467613834773723740590533316085,
  390.                 14605578.08768506808414169982791359218571,
  391.                 3481712.15498064590882071018964774556468,
  392.                 601859.6171681098786670226533699352302507,
  393.                 75999.29304014542649875303443598909137092,
  394.                 6955.999602515376140356310115515198987526,
  395.                 449.9445569063168119446858607650988409623,
  396.                 19.51992788247617482847860966235652136208,
  397.                 0.5098416655656676188125178644804694509993,
  398.                 0.006061842346248906525783753964555936883222
  399.             };
  400.             return evaluateRational(num, DENOM, z);
  401.         }

  402.         /**
  403.          * Evaluate the rational number as two polynomials.
  404.          *
  405.          * <p>Adapted from {@code boost/math/tools/detail/rational_horner3_13.hpp}.
  406.          * Note: There are 3 variations of the unrolled rational evaluation.
  407.          * These methods change the order based on the sign of x. This
  408.          * should be used for the Lanczos code as this comment in
  409.          * {@code boost/math/tools/rational.hpp} notes:
  410.          *
  411.          * <blockquote>
  412.          * However, there
  413.          * are some tricks we can use to prevent overflow that might otherwise
  414.          * occur in polynomial evaluation, if z is large.  This is important
  415.          * in our Lanczos code for example.
  416.          * </blockquote>
  417.          *
  418.          * @param a Coefficients of the numerator polynomial
  419.          * @param b Coefficients of the denominator polynomial
  420.          * @param x Value
  421.          * @return the rational number
  422.          */
  423.         private static double evaluateRational(double[] a, int[] b, double x) {
  424.             // The choice of algorithm in Boost is based on the compiler
  425.             // to suite the available optimisations.
  426.             //
  427.             // Tests against rational_horner1_13.hpp which uses a first order
  428.             // Horner method (no x*x term) show only minor variations in
  429.             // error. rational_horner2_13.hpp has the same second order Horner
  430.             // method with different code layout of the same sum.

  431.             // rational_horner3_13.hpp
  432.             if (x <= 1) {
  433.                 final double x2 = x * x;
  434.                 double t0 = a[12] * x2 + a[10];
  435.                 double t1 = a[11] * x2 + a[9];
  436.                 double t2 = b[12] * x2 + b[10];
  437.                 double t3 = b[11] * x2 + b[9];
  438.                 t0 *= x2;
  439.                 t1 *= x2;
  440.                 t2 *= x2;
  441.                 t3 *= x2;
  442.                 t0 += a[8];
  443.                 t1 += a[7];
  444.                 t2 += b[8];
  445.                 t3 += b[7];
  446.                 t0 *= x2;
  447.                 t1 *= x2;
  448.                 t2 *= x2;
  449.                 t3 *= x2;
  450.                 t0 += a[6];
  451.                 t1 += a[5];
  452.                 t2 += b[6];
  453.                 t3 += b[5];
  454.                 t0 *= x2;
  455.                 t1 *= x2;
  456.                 t2 *= x2;
  457.                 t3 *= x2;
  458.                 t0 += a[4];
  459.                 t1 += a[3];
  460.                 t2 += b[4];
  461.                 t3 += b[3];
  462.                 t0 *= x2;
  463.                 t1 *= x2;
  464.                 t2 *= x2;
  465.                 t3 *= x2;
  466.                 t0 += a[2];
  467.                 t1 += a[1];
  468.                 t2 += b[2];
  469.                 t3 += b[1];
  470.                 t0 *= x2;
  471.                 t2 *= x2;
  472.                 t0 += a[0];
  473.                 t2 += b[0];
  474.                 t1 *= x;
  475.                 t3 *= x;
  476.                 return (t0 + t1) / (t2 + t3);
  477.             }
  478.             final double z = 1 / x;
  479.             final double z2 = 1 / (x * x);
  480.             double t0 = a[0] * z2 + a[2];
  481.             double t1 = a[1] * z2 + a[3];
  482.             double t2 = b[0] * z2 + b[2];
  483.             double t3 = b[1] * z2 + b[3];
  484.             t0 *= z2;
  485.             t1 *= z2;
  486.             t2 *= z2;
  487.             t3 *= z2;
  488.             t0 += a[4];
  489.             t1 += a[5];
  490.             t2 += b[4];
  491.             t3 += b[5];
  492.             t0 *= z2;
  493.             t1 *= z2;
  494.             t2 *= z2;
  495.             t3 *= z2;
  496.             t0 += a[6];
  497.             t1 += a[7];
  498.             t2 += b[6];
  499.             t3 += b[7];
  500.             t0 *= z2;
  501.             t1 *= z2;
  502.             t2 *= z2;
  503.             t3 *= z2;
  504.             t0 += a[8];
  505.             t1 += a[9];
  506.             t2 += b[8];
  507.             t3 += b[9];
  508.             t0 *= z2;
  509.             t1 *= z2;
  510.             t2 *= z2;
  511.             t3 *= z2;
  512.             t0 += a[10];
  513.             t1 += a[11];
  514.             t2 += b[10];
  515.             t3 += b[11];
  516.             t0 *= z2;
  517.             t2 *= z2;
  518.             t0 += a[12];
  519.             t2 += b[12];
  520.             t1 *= z;
  521.             t3 *= z;
  522.             return (t0 + t1) / (t2 + t3);
  523.         }

  524.         // Not implemented:
  525.         // lanczos_sum_near_1
  526.         // lanczos_sum_near_2
  527.     }

  528.     /** Private constructor. */
  529.     private BoostGamma() {
  530.         // intentionally empty.
  531.     }

  532.     /**
  533.      * All factorials that are representable as a double.
  534.      * This data is exposed for testing.
  535.      *
  536.      * @return factorials
  537.      */
  538.     static double[] getFactorials() {
  539.         return FACTORIAL.clone();
  540.     }

  541.     /**
  542.      * Returns the factorial of n.
  543.      * This is unchecked as an index out of bound exception will occur
  544.      * if the value n is not within the range [0, 170].
  545.      * This function is exposed for use in {@link BoostBeta}.
  546.      *
  547.      * @param n Argument n (must be in [0, 170])
  548.      * @return n!
  549.      */
  550.     static double uncheckedFactorial(int n) {
  551.         return FACTORIAL[n];
  552.     }

  553.     /**
  554.      * Gamma function.
  555.      *
  556.      * <p>For small {@code z} this is based on the <em>NSWC Library of Mathematics
  557.      * Subroutines</em> double precision implementation, {@code DGAMMA}.
  558.      *
  559.      * <p>For large {@code z} this is an implementation of the Boost C++ tgamma
  560.      * function with Lanczos support.
  561.      *
  562.      * <p>Integers are handled using a look-up table of factorials.
  563.      *
  564.      * <p>Note: The Boost C++ implementation uses the Lanczos sum for all {@code z}.
  565.      * When promotion of double to long double is not available this has larger
  566.      * errors than the double precision specific NSWC implementation. For larger
  567.      * {@code z} the Boost C++ Lanczos implementation incorporates the sqrt(2 pi)
  568.      * factor and has lower error than the implementation using the
  569.      * {@link LanczosApproximation} class.
  570.      *
  571.      * @param z Argument z
  572.      * @return gamma value
  573.      */
  574.     static double tgamma(double z) {
  575.         // Handle integers
  576.         if (Math.rint(z) == z) {
  577.             if (z <= 0) {
  578.                 // Pole error
  579.                 return Double.NaN;
  580.             }
  581.             if (z <= MAX_GAMMA_Z) {
  582.                 // Gamma(n) = (n-1)!
  583.                 return FACTORIAL[(int) z - 1];
  584.             }
  585.             // Overflow
  586.             return Double.POSITIVE_INFINITY;
  587.         }

  588.         if (Math.abs(z) <= LANCZOS_THRESHOLD) {
  589.             // Small z
  590.             // NSWC Library of Mathematics Subroutines
  591.             // Note:
  592.             // This does not benefit from using extended precision to track the sum (t).
  593.             // Extended precision on the product reduces the error but the majority
  594.             // of error remains in InvGamma1pm1.

  595.             if (z >= 1) {
  596.                 /*
  597.                  * From the recurrence relation
  598.                  * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
  599.                  * then
  600.                  * Gamma(t) = 1 / [1 + InvGamma1pm1.value(t - 1)],
  601.                  * where t = x - n. This means that t must satisfy
  602.                  * -0.5 <= t - 1 <= 1.5.
  603.                  */
  604.                 double prod = 1;
  605.                 double t = z;
  606.                 while (t > 2.5) {
  607.                     t -= 1;
  608.                     prod *= t;
  609.                 }
  610.                 return prod / (1 + InvGamma1pm1.value(t - 1));
  611.             }
  612.             /*
  613.              * From the recurrence relation
  614.              * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
  615.              * then
  616.              * Gamma(x + n + 1) = 1 / [1 + InvGamma1pm1.value(x + n)],
  617.              * which requires -0.5 <= x + n <= 1.5.
  618.              */
  619.             double prod = z;
  620.             double t = z;
  621.             while (t < -0.5) {
  622.                 t += 1;
  623.                 prod *= t;
  624.             }
  625.             return 1 / (prod * (1 + InvGamma1pm1.value(t)));
  626.         }

  627.         // Large non-integer z
  628.         // Boost C++ tgamma implementation

  629.         if (z < 0) {
  630.             /*
  631.              * From the reflection formula
  632.              * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
  633.              * and the recurrence relation
  634.              * Gamma(1 - x) = -x * Gamma(-x),
  635.              * it is found
  636.              * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
  637.              */
  638.             return -Math.PI / (sinpx(z) * tgamma(-z));
  639.         } else if (z > MAX_GAMMA_Z + 1) {
  640.             // Addition to the Boost code: Simple overflow detection
  641.             return Double.POSITIVE_INFINITY;
  642.         }

  643.         double result = Lanczos.lanczosSum(z);
  644.         final double zgh = z + Lanczos.GMH;
  645.         final double lzgh = Math.log(zgh);
  646.         if (z * lzgh > LOG_MAX_VALUE) {
  647.             // we're going to overflow unless this is done with care:

  648.             // Updated
  649.             // Check for overflow removed:
  650.             // if (lzgh * z / 2 > LOG_MAX_VALUE) ... overflow
  651.             // This is replaced by checking z > MAX_FACTORIAL + 2

  652.             final double hp = Math.pow(zgh, (z / 2) - 0.25);
  653.             result *= hp / Math.exp(zgh);
  654.             // Check for overflow has been removed:
  655.             // if (Double.MAX_VALUE / hp < result) ... overflow
  656.             result *= hp;
  657.         } else {
  658.             result *= Math.pow(zgh, z - 0.5) / Math.exp(zgh);
  659.         }

  660.         return result;
  661.     }

  662.     /**
  663.      * Ad hoc function calculates x * sin(pi * x), taking extra care near when x is
  664.      * near a whole number.
  665.      *
  666.      * @param x Value (assumed to be negative)
  667.      * @return x * sin(pi * x)
  668.      */
  669.     static double sinpx(double x) {
  670.         int sign = 1;
  671.         // This is always called with a negative
  672.         // if (x < 0)
  673.         x = -x;
  674.         double fl = Math.floor(x);
  675.         double dist;
  676.         if (isOdd(fl)) {
  677.             fl += 1;
  678.             dist = fl - x;
  679.             sign = -sign;
  680.         } else {
  681.             dist = x - fl;
  682.         }
  683.         if (dist > 0.5f) {
  684.             dist = 1 - dist;
  685.         }
  686.         final double result = Math.sin(dist * Math.PI);
  687.         return sign * x * result;
  688.     }

  689.     /**
  690.      * Checks if the value is odd.
  691.      *
  692.      * @param v Value (assumed to be positive and an integer)
  693.      * @return true if odd
  694.      */
  695.     private static boolean isOdd(double v) {
  696.         // Note:
  697.         // Any value larger than 2^53 should be even.
  698.         // If the input is positive then truncation of extreme doubles (>2^63)
  699.         // to the primitive long creates an odd value: 2^63-1.
  700.         // This is corrected by inverting the sign of v and the extreme is even: -2^63.
  701.         // This function is never called when the argument is this large
  702.         // as this is a pole error in tgamma so the effect is never observed.
  703.         // However the isOdd function is correct for all positive finite v.
  704.         return (((long) -v) & 0x1) == 1;
  705.     }

  706.     /**
  707.      * Log Gamma function.
  708.      * Defined as the natural logarithm of the absolute value of tgamma(z).
  709.      *
  710.      * @param z Argument z
  711.      * @return log gamma value
  712.      */
  713.     static double lgamma(double z) {
  714.         return lgamma(z, null);
  715.     }

  716.     /**
  717.      * Log Gamma function.
  718.      * Defined as the natural logarithm of the absolute value of tgamma(z).
  719.      *
  720.      * @param z Argument z
  721.      * @param sign If a non-zero length array the first index is set on output to the sign of tgamma(z)
  722.      * @return log gamma value
  723.      */
  724.     static double lgamma(double z, int[] sign) {
  725.         double result;
  726.         int sresult = 1;
  727.         if (z <= -ROOT_EPSILON) {
  728.             // reflection formula:
  729.             if (Math.rint(z) == z) {
  730.                 // Pole error
  731.                 return Double.NaN;
  732.             }

  733.             double t = sinpx(z);
  734.             z = -z;
  735.             if (t < 0) {
  736.                 t = -t;
  737.             } else {
  738.                 sresult = -sresult;
  739.             }

  740.             // This summation can have large magnitudes with opposite signs.
  741.             // Use an extended precision sum to reduce cancellation.
  742.             result = Sum.of(-lgamma(z)).add(-Math.log(t)).add(LOG_PI).getAsDouble();

  743.         } else if (z < ROOT_EPSILON) {
  744.             if (z == 0) {
  745.                 // Pole error
  746.                 return Double.NaN;
  747.             }
  748.             if (4 * Math.abs(z) < EPSILON) {
  749.                 result = -Math.log(Math.abs(z));
  750.             } else {
  751.                 result = Math.log(Math.abs(1 / z - EULER));
  752.             }
  753.             if (z < 0) {
  754.                 sresult = -1;
  755.             }
  756.         } else if (z < 15) {
  757.             result = lgammaSmall(z, z - 1, z - 2);
  758.         // The z > 3 condition is always true
  759.         //} else if (z > 3 && z < 100) {
  760.         } else if (z < 100) {
  761.             // taking the log of tgamma reduces the error, no danger of overflow here:
  762.             result = Math.log(tgamma(z));
  763.         } else {
  764.             // regular evaluation:
  765.             final double zgh = z + Lanczos.GMH;
  766.             result = Math.log(zgh) - 1;
  767.             result *= z - 0.5f;
  768.             //
  769.             // Only add on the lanczos sum part if we're going to need it:
  770.             //
  771.             if (result * EPSILON < 20) {
  772.                 result += Math.log(Lanczos.lanczosSumExpGScaled(z));
  773.             }
  774.         }

  775.         if (nonZeroLength(sign)) {
  776.             sign[0] = sresult;
  777.         }
  778.         return result;
  779.     }

  780.     /**
  781.      * Log Gamma function for small z.
  782.      *
  783.      * @param z Argument z
  784.      * @param zm1 {@code z - 1}
  785.      * @param zm2 {@code z - 2}
  786.      * @return log gamma value
  787.      */
  788.     private static double lgammaSmall(double z, double zm1, double zm2) {
  789.         // This version uses rational approximations for small
  790.         // values of z accurate enough for 64-bit mantissas
  791.         // (80-bit long doubles), works well for 53-bit doubles as well.

  792.         // Updated to use an extended precision sum
  793.         final Sum result = Sum.create();

  794.         // Note:
  795.         // Removed z < EPSILON branch.
  796.         // The function is called
  797.         // from lgamma:
  798.         //   ROOT_EPSILON <= z < 15
  799.         // from tgamma1pm1:
  800.         //   1.5 <= z < 2
  801.         //   1 <= z < 3

  802.         if ((zm1 == 0) || (zm2 == 0)) {
  803.             // nothing to do, result is zero....
  804.             return 0;
  805.         } else if (z > 2) {
  806.             //
  807.             // Begin by performing argument reduction until
  808.             // z is in [2,3):
  809.             //
  810.             if (z >= 3) {
  811.                 do {
  812.                     z -= 1;
  813.                     result.add(Math.log(z));
  814.                 } while (z >= 3);
  815.                 // Update zm2, we need it below:
  816.                 zm2 = z - 2;
  817.             }

  818.             //
  819.             // Use the following form:
  820.             //
  821.             // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
  822.             //
  823.             // where R(z-2) is a rational approximation optimised for
  824.             // low absolute error - as long as its absolute error
  825.             // is small compared to the constant Y - then any rounding
  826.             // error in its computation will get wiped out.
  827.             //
  828.             // R(z-2) has the following properties:
  829.             //
  830.             // At double: Max error found:                    4.231e-18
  831.             // At long double: Max error found:               1.987e-21
  832.             // Maximum Deviation Found (approximation error): 5.900e-24
  833.             //
  834.             double P;
  835.             P = -0.324588649825948492091e-4;
  836.             P = -0.541009869215204396339e-3 + P * zm2;
  837.             P = -0.259453563205438108893e-3 + P * zm2;
  838.             P =  0.172491608709613993966e-1 + P * zm2;
  839.             P =  0.494103151567532234274e-1 + P * zm2;
  840.             P =   0.25126649619989678683e-1 + P * zm2;
  841.             P = -0.180355685678449379109e-1 + P * zm2;
  842.             double Q;
  843.             Q = -0.223352763208617092964e-6;
  844.             Q =  0.224936291922115757597e-3 + Q * zm2;
  845.             Q =   0.82130967464889339326e-2 + Q * zm2;
  846.             Q =  0.988504251128010129477e-1 + Q * zm2;
  847.             Q =   0.541391432071720958364e0 + Q * zm2;
  848.             Q =   0.148019669424231326694e1 + Q * zm2;
  849.             Q =   0.196202987197795200688e1 + Q * zm2;
  850.             Q =                       0.1e1 + Q * zm2;

  851.             final float Y = 0.158963680267333984375e0f;

  852.             final double r = zm2 * (z + 1);
  853.             final double R = P / Q;

  854.             result.addProduct(r, Y).addProduct(r, R);
  855.         } else {
  856.             //
  857.             // If z is less than 1 use recurrence to shift to
  858.             // z in the interval [1,2]:
  859.             //
  860.             if (z < 1) {
  861.                 result.add(-Math.log(z));
  862.                 zm2 = zm1;
  863.                 zm1 = z;
  864.                 z += 1;
  865.             }
  866.             //
  867.             // Two approximations, one for z in [1,1.5] and
  868.             // one for z in [1.5,2]:
  869.             //
  870.             if (z <= 1.5) {
  871.                 //
  872.                 // Use the following form:
  873.                 //
  874.                 // lgamma(z) = (z-1)(z-2)(Y + R(z-1
  875.                 //
  876.                 // where R(z-1) is a rational approximation optimised for
  877.                 // low absolute error - as long as its absolute error
  878.                 // is small compared to the constant Y - then any rounding
  879.                 // error in its computation will get wiped out.
  880.                 //
  881.                 // R(z-1) has the following properties:
  882.                 //
  883.                 // At double precision: Max error found:                1.230011e-17
  884.                 // At 80-bit long double precision:   Max error found:  5.631355e-21
  885.                 // Maximum Deviation Found:                             3.139e-021
  886.                 // Expected Error Term:                                 3.139e-021

  887.                 //
  888.                 final float Y = 0.52815341949462890625f;

  889.                 double P;
  890.                 P = -0.100346687696279557415e-2;
  891.                 P = -0.240149820648571559892e-1 + P * zm1;
  892.                 P =  -0.158413586390692192217e0 + P * zm1;
  893.                 P =  -0.406567124211938417342e0 + P * zm1;
  894.                 P =  -0.414983358359495381969e0 + P * zm1;
  895.                 P = -0.969117530159521214579e-1 + P * zm1;
  896.                 P =  0.490622454069039543534e-1 + P * zm1;
  897.                 double Q;
  898.                 Q = 0.195768102601107189171e-2;
  899.                 Q = 0.577039722690451849648e-1 + Q * zm1;
  900.                 Q =  0.507137738614363510846e0 + Q * zm1;
  901.                 Q =  0.191415588274426679201e1 + Q * zm1;
  902.                 Q =  0.348739585360723852576e1 + Q * zm1;
  903.                 Q =  0.302349829846463038743e1 + Q * zm1;
  904.                 Q =                      0.1e1 + Q * zm1;

  905.                 final double r = P / Q;
  906.                 final double prefix = zm1 * zm2;

  907.                 result.addProduct(prefix, Y).addProduct(prefix, r);
  908.             } else {
  909.                 //
  910.                 // Use the following form:
  911.                 //
  912.                 // lgamma(z) = (2-z)(1-z)(Y + R(2-z
  913.                 //
  914.                 // where R(2-z) is a rational approximation optimised for
  915.                 // low absolute error - as long as its absolute error
  916.                 // is small compared to the constant Y - then any rounding
  917.                 // error in its computation will get wiped out.
  918.                 //
  919.                 // R(2-z) has the following properties:
  920.                 //
  921.                 // At double precision, max error found:              1.797565e-17
  922.                 // At 80-bit long double precision, max error found:  9.306419e-21
  923.                 // Maximum Deviation Found:                           2.151e-021
  924.                 // Expected Error Term:                               2.150e-021
  925.                 //
  926.                 final float Y = 0.452017307281494140625f;

  927.                 final double mzm2 = -zm2;
  928.                 double P;
  929.                 P =  0.431171342679297331241e-3;
  930.                 P = -0.850535976868336437746e-2 + P * mzm2;
  931.                 P =  0.542809694055053558157e-1 + P * mzm2;
  932.                 P =  -0.142440390738631274135e0 + P * mzm2;
  933.                 P =   0.144216267757192309184e0 + P * mzm2;
  934.                 P = -0.292329721830270012337e-1 + P * mzm2;
  935.                 double Q;
  936.                 Q = -0.827193521891290553639e-6;
  937.                 Q = -0.100666795539143372762e-2 + Q * mzm2;
  938.                 Q =   0.25582797155975869989e-1 + Q * mzm2;
  939.                 Q =  -0.220095151814995745555e0 + Q * mzm2;
  940.                 Q =   0.846973248876495016101e0 + Q * mzm2;
  941.                 Q =  -0.150169356054485044494e1 + Q * mzm2;
  942.                 Q =                       0.1e1 + Q * mzm2;
  943.                 final double r = zm2 * zm1;
  944.                 final double R = P / Q;

  945.                 result.addProduct(r, Y).addProduct(r, R);
  946.             }
  947.         }
  948.         return result.getAsDouble();
  949.     }

  950.     /**
  951.      * Calculates tgamma(1+dz)-1.
  952.      *
  953.      * @param dz Argument
  954.      * @return tgamma(1+dz)-1
  955.      */
  956.     static double tgamma1pm1(double dz) {
  957.         //
  958.         // This helper calculates tgamma(1+dz)-1 without cancellation errors,
  959.         // used by the upper incomplete gamma with z < 1:
  960.         //
  961.         final double result;
  962.         if (dz < 0) {
  963.             if (dz < -0.5) {
  964.                 // Best method is simply to subtract 1 from tgamma:
  965.                 result = tgamma(1 + dz) - 1;
  966.             } else {
  967.                 // Use expm1 on lgamma:
  968.                 result = Math.expm1(-Math.log1p(dz) + lgammaSmall(dz + 2, dz + 1, dz));
  969.             }
  970.         } else {
  971.             if (dz < 2) {
  972.                 // Use expm1 on lgamma:
  973.                 result = Math.expm1(lgammaSmall(dz + 1, dz, dz - 1));
  974.             } else {
  975.                 // Best method is simply to subtract 1 from tgamma:
  976.                 result = tgamma(1 + dz) - 1;
  977.             }
  978.         }

  979.         return result;
  980.     }

  981.     /**
  982.      * Full upper incomplete gamma.
  983.      *
  984.      * @param a Argument a
  985.      * @param x Argument x
  986.      * @return upper gamma value
  987.      */
  988.     static double tgamma(double a, double x) {
  989.         return gammaIncompleteImp(a, x, false, true, Policy.getDefault());
  990.     }

  991.     /**
  992.      * Full upper incomplete gamma.
  993.      *
  994.      * @param a Argument a
  995.      * @param x Argument x
  996.      * @param policy Function evaluation policy
  997.      * @return upper gamma value
  998.      */
  999.     static double tgamma(double a, double x, Policy policy) {
  1000.         return gammaIncompleteImp(a, x, false, true, policy);
  1001.     }

  1002.     /**
  1003.      * Full lower incomplete gamma.
  1004.      *
  1005.      * @param a Argument a
  1006.      * @param x Argument x
  1007.      * @return lower gamma value
  1008.      */
  1009.     static double tgammaLower(double a, double x) {
  1010.         return gammaIncompleteImp(a, x, false, false, Policy.getDefault());
  1011.     }

  1012.     /**
  1013.      * Full lower incomplete gamma.
  1014.      *
  1015.      * @param a Argument a
  1016.      * @param x Argument x
  1017.      * @param policy Function evaluation policy
  1018.      * @return lower gamma value
  1019.      */
  1020.     static double tgammaLower(double a, double x, Policy policy) {
  1021.         return gammaIncompleteImp(a, x, false, false, policy);
  1022.     }

  1023.     /**
  1024.      * Regularised upper incomplete gamma.
  1025.      *
  1026.      * @param a Argument a
  1027.      * @param x Argument x
  1028.      * @return q
  1029.      */
  1030.     static double gammaQ(double a, double x) {
  1031.         return gammaIncompleteImp(a, x, true, true, Policy.getDefault());
  1032.     }

  1033.     /**
  1034.      * Regularised upper incomplete gamma.
  1035.      *
  1036.      * @param a Argument a
  1037.      * @param x Argument x
  1038.      * @param policy Function evaluation policy
  1039.      * @return q
  1040.      */
  1041.     static double gammaQ(double a, double x, Policy policy) {
  1042.         return gammaIncompleteImp(a, x, true, true, policy);
  1043.     }

  1044.     /**
  1045.      * Regularised lower incomplete gamma.
  1046.      *
  1047.      * @param a Argument a
  1048.      * @param x Argument x
  1049.      * @return p
  1050.      */
  1051.     static double gammaP(double a, double x) {
  1052.         return gammaIncompleteImp(a, x, true, false, Policy.getDefault());
  1053.     }

  1054.     /**
  1055.      * Regularised lower incomplete gamma.
  1056.      *
  1057.      * @param a Argument a
  1058.      * @param x Argument x
  1059.      * @param policy Function evaluation policy
  1060.      * @return p
  1061.      */
  1062.     static double gammaP(double a, double x, Policy policy) {
  1063.         return gammaIncompleteImp(a, x, true, false, policy);
  1064.     }

  1065.     /**
  1066.      * Derivative of the regularised lower incomplete gamma.
  1067.      * <p>\( \frac{e^{-x} x^{a-1}}{\Gamma(a)} \)
  1068.      *
  1069.      * <p>Adapted from {@code boost::math::detail::gamma_p_derivative_imp}
  1070.      *
  1071.      * @param a Argument a
  1072.      * @param x Argument x
  1073.      * @return p derivative
  1074.      */
  1075.     static double gammaPDerivative(double a, double x) {
  1076.         //
  1077.         // Usual error checks first:
  1078.         //
  1079.         if (Double.isNaN(a) || Double.isNaN(x) || a <= 0 || x < 0) {
  1080.             return Double.NaN;
  1081.         }
  1082.         //
  1083.         // Now special cases:
  1084.         //
  1085.         if (x == 0) {
  1086.             if (a > 1) {
  1087.                 return 0;
  1088.             }
  1089.             return (a == 1) ? 1 : Double.POSITIVE_INFINITY;
  1090.         }
  1091.         //
  1092.         // Normal case:
  1093.         //
  1094.         double f1 = regularisedGammaPrefix(a, x);
  1095.         if (f1 == 0) {
  1096.             // Underflow in calculation, use logs instead:
  1097.             f1 = a * Math.log(x) - x - lgamma(a) - Math.log(x);
  1098.             f1 = Math.exp(f1);
  1099.         } else {
  1100.             // Will overflow when (x < 1) && (Double.MAX_VALUE * x < f1).
  1101.             // There is no exception for this case so just return the result.
  1102.             f1 /= x;
  1103.         }

  1104.         return f1;
  1105.     }

  1106.     /**
  1107.      * Main incomplete gamma entry point, handles all four incomplete gammas.
  1108.      * Adapted from {@code boost::math::detail::gamma_incomplete_imp}.
  1109.      *
  1110.      * <p>The Boost code has a pointer {@code p_derivative} that can be set to the
  1111.      * value of the derivative. This is used for the inverse incomplete
  1112.      * gamma functions {@code gamma_(p|q)_inv_imp}. It is not required for the forward
  1113.      * evaluation functions.
  1114.      *
  1115.      * @param a Argument a
  1116.      * @param x Argument x
  1117.      * @param normalised true to compute the regularised value
  1118.      * @param invert true to compute the upper value Q (default is lower value P)
  1119.      * @param pol Function evaluation policy
  1120.      * @return gamma value
  1121.      */
  1122.     private static double gammaIncompleteImp(double a, double x,
  1123.             boolean normalised, boolean invert, Policy pol) {
  1124.         if (Double.isNaN(a) || Double.isNaN(x) || a <= 0 || x < 0) {
  1125.             return Double.NaN;
  1126.         }

  1127.         double result = 0;

  1128.         if (a >= MAX_FACTORIAL && !normalised) {
  1129.             //
  1130.             // When we're computing the non-normalized incomplete gamma
  1131.             // and a is large the result is rather hard to compute unless
  1132.             // we use logs. There are really two options - if x is a long
  1133.             // way from a in value then we can reliably use methods 2 and 4
  1134.             // below in logarithmic form and go straight to the result.
  1135.             // Otherwise we let the regularized gamma take the strain
  1136.             // (the result is unlikely to underflow in the central region anyway)
  1137.             // and combine with lgamma in the hopes that we get a finite result.
  1138.             //

  1139.             if (invert && (a * 4 < x)) {
  1140.                 // This is method 4 below, done in logs:
  1141.                 result = a * Math.log(x) - x;
  1142.                 result += Math.log(upperGammaFraction(a, x, pol));
  1143.             } else if (!invert && (a > 4 * x)) {
  1144.                 // This is method 2 below, done in logs:
  1145.                 result = a * Math.log(x) - x;
  1146.                 result += Math.log(lowerGammaSeries(a, x, 0, pol) / a);
  1147.             } else {
  1148.                 result = gammaIncompleteImp(a, x, true, invert, pol);
  1149.                 if (result == 0) {
  1150.                     if (invert) {
  1151.                         // Try http://functions.wolfram.com/06.06.06.0039.01
  1152.                         result = 1 + 1 / (12 * a) + 1 / (288 * a * a);
  1153.                         result = Math.log(result) - a + (a - 0.5f) * Math.log(a) + LOG_ROOT_TWO_PI;
  1154.                     } else {
  1155.                         // This is method 2 below, done in logs, we're really outside the
  1156.                         // range of this method, but since the result is almost certainly
  1157.                         // infinite, we should probably be OK:
  1158.                         result = a * Math.log(x) - x;
  1159.                         result += Math.log(lowerGammaSeries(a, x, 0, pol) / a);
  1160.                     }
  1161.                 } else {
  1162.                     result = Math.log(result) + lgamma(a);
  1163.                 }
  1164.             }
  1165.             // If result is > log(MAX_VALUE) the result will overflow.
  1166.             // There is no exception for this case so just return the result.
  1167.             return Math.exp(result);
  1168.         }

  1169.         final boolean isInt;
  1170.         final boolean isHalfInt;
  1171.         // Update. x must be safe for exp(-x). Change to -x > LOG_MIN_VALUE.
  1172.         final boolean isSmallA = (a < 30) && (a <= x + 1) && (-x > LOG_MIN_VALUE);
  1173.         if (isSmallA) {
  1174.             final double fa = Math.floor(a);
  1175.             isInt = fa == a;
  1176.             isHalfInt = !isInt && (Math.abs(fa - a) == 0.5f);
  1177.         } else {
  1178.             isInt = isHalfInt = false;
  1179.         }

  1180.         final int evalMethod;

  1181.         if (isInt && (x > 0.6)) {
  1182.             // calculate Q via finite sum:
  1183.             invert = !invert;
  1184.             evalMethod = 0;
  1185.         } else if (isHalfInt && (x > 0.2)) {
  1186.             // calculate Q via finite sum for half integer a:
  1187.             invert = !invert;
  1188.             evalMethod = 1;
  1189.         } else if ((x < ROOT_EPSILON) && (a > 1)) {
  1190.             evalMethod = 6;
  1191.         } else if ((x > 1000) && (a < x * 0.75f)) {
  1192.             // Note:
  1193.             // The branch is used in Boost when:
  1194.             // ((x > 1000) && ((a < x) || (Math.abs(a - 50) / x < 1)))
  1195.             //
  1196.             // This case was added after Boost 1_68_0.
  1197.             // See: https://github.com/boostorg/math/issues/168
  1198.             //
  1199.             // When using only double precision for the evaluation
  1200.             // it is a source of error when a ~ z as the asymptotic approximation
  1201.             // sums terms t_n+1 = t_n * (a - n - 1) / z starting from t_0 = 1.
  1202.             // These terms are close to 1 when a ~ z and the sum has many terms
  1203.             // with reduced precision.
  1204.             // This has been updated to allow only cases with fast convergence.
  1205.             // It will be used when x -> infinity and a << x.

  1206.             // calculate Q via asymptotic approximation:
  1207.             invert = !invert;
  1208.             evalMethod = 7;

  1209.         } else if (x < 0.5) {
  1210.             //
  1211.             // Changeover criterion chosen to give a changeover at Q ~ 0.33
  1212.             //
  1213.             if (-0.4 / Math.log(x) < a) {
  1214.                 // Compute P
  1215.                 evalMethod = 2;
  1216.             } else {
  1217.                 evalMethod = 3;
  1218.             }
  1219.         } else if (x < 1.1) {
  1220.             //
  1221.             // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
  1222.             //
  1223.             if (x * 0.75f < a) {
  1224.                 // Compute P
  1225.                 evalMethod = 2;
  1226.             } else {
  1227.                 evalMethod = 3;
  1228.             }
  1229.         } else {
  1230.             //
  1231.             // Begin by testing whether we're in the "bad" zone
  1232.             // where the result will be near 0.5 and the usual
  1233.             // series and continued fractions are slow to converge:
  1234.             //
  1235.             boolean useTemme = false;
  1236.             if (normalised && (a > 20)) {
  1237.                 final double sigma = Math.abs((x - a) / a);
  1238.                 if (a > 200) {
  1239.                     //
  1240.                     // This limit is chosen so that we use Temme's expansion
  1241.                     // only if the result would be larger than about 10^-6.
  1242.                     // Below that the regular series and continued fractions
  1243.                     // converge OK, and if we use Temme's method we get increasing
  1244.                     // errors from the dominant erfc term as its (inexact) argument
  1245.                     // increases in magnitude.
  1246.                     //
  1247.                     if (20 / a > sigma * sigma) {
  1248.                         useTemme = true;
  1249.                     }
  1250.                 } else {
  1251.                     // Note in this zone we can't use Temme's expansion for
  1252.                     // types longer than an 80-bit real:
  1253.                     // it would require too many terms in the polynomials.
  1254.                     if (sigma < 0.4) {
  1255.                         useTemme = true;
  1256.                     }
  1257.                 }
  1258.             }
  1259.             if (useTemme) {
  1260.                 evalMethod = 5;
  1261.             } else {
  1262.                 //
  1263.                 // Regular case where the result will not be too close to 0.5.
  1264.                 //
  1265.                 // Changeover here occurs at P ~ Q ~ 0.5
  1266.                 // Note that series computation of P is about x2 faster than continued fraction
  1267.                 // calculation of Q, so try and use the CF only when really necessary,
  1268.                 // especially for small x.
  1269.                 //
  1270.                 if (x - (1 / (3 * x)) < a) {
  1271.                     evalMethod = 2;
  1272.                 } else {
  1273.                     evalMethod = 4;
  1274.                     invert = !invert;
  1275.                 }
  1276.             }
  1277.         }

  1278.         switch (evalMethod) {
  1279.         case 0:
  1280.             result = finiteGammaQ(a, x);
  1281.             if (!normalised) {
  1282.                 result *= tgamma(a);
  1283.             }
  1284.             break;
  1285.         case 1:
  1286.             result = finiteHalfGammaQ(a, x);
  1287.             if (!normalised) {
  1288.                 result *= tgamma(a);
  1289.             }
  1290.             break;
  1291.         case 2:
  1292.             // Compute P:
  1293.             result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
  1294.             if (result != 0) {
  1295.                 //
  1296.                 // If we're going to be inverting the result then we can
  1297.                 // reduce the number of series evaluations by quite
  1298.                 // a few iterations if we set an initial value for the
  1299.                 // series sum based on what we'll end up subtracting it from
  1300.                 // at the end.
  1301.                 // Have to be careful though that this optimization doesn't
  1302.                 // lead to spurious numeric overflow. Note that the
  1303.                 // scary/expensive overflow checks below are more often
  1304.                 // than not bypassed in practice for "sensible" input
  1305.                 // values:
  1306.                 //

  1307.                 double initValue = 0;
  1308.                 boolean optimisedInvert = false;
  1309.                 if (invert) {
  1310.                     initValue = normalised ? 1 : tgamma(a);
  1311.                     if (normalised || (result >= 1) || (Double.MAX_VALUE * result > initValue)) {
  1312.                         initValue /= result;
  1313.                         if (normalised || (a < 1) || (Double.MAX_VALUE / a > initValue)) {
  1314.                             initValue *= -a;
  1315.                             optimisedInvert = true;
  1316.                         } else {
  1317.                             initValue = 0;
  1318.                         }
  1319.                     } else {
  1320.                         initValue = 0;
  1321.                     }
  1322.                 }
  1323.                 result *= lowerGammaSeries(a, x, initValue, pol) / a;
  1324.                 if (optimisedInvert) {
  1325.                     invert = false;
  1326.                     result = -result;
  1327.                 }
  1328.             }
  1329.             break;
  1330.         case 3:
  1331.             // Compute Q:
  1332.             invert = !invert;
  1333.             final double[] g = {0};
  1334.             result = tgammaSmallUpperPart(a, x, pol, g, invert);
  1335.             invert = false;
  1336.             if (normalised) {
  1337.                 // Addition to the Boost code:
  1338.                 if (g[0] == Double.POSITIVE_INFINITY) {
  1339.                     // Very small a will overflow gamma(a). Resort to logs.
  1340.                     // This method requires improvement as the error is very large.
  1341.                     // It is better than returning zero for a non-zero result.
  1342.                     result = Math.exp(Math.log(result) - lgamma(a));
  1343.                 } else {
  1344.                     result /= g[0];
  1345.                 }
  1346.             }
  1347.             break;
  1348.         case 4:
  1349.             // Compute Q:
  1350.             result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
  1351.             if (result != 0) {
  1352.                 result *= upperGammaFraction(a, x, pol);
  1353.             }
  1354.             break;
  1355.         case 5:
  1356.             // Call 53-bit version
  1357.             result = igammaTemmeLarge(a, x);
  1358.             if (x >= a) {
  1359.                 invert = !invert;
  1360.             }
  1361.             break;
  1362.         case 6:
  1363.             // x is so small that P is necessarily very small too, use
  1364.             // http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
  1365.             if (normalised) {
  1366.                 // If tgamma overflows then result = 0
  1367.                 result = Math.pow(x, a) / tgamma(a + 1);
  1368.             } else {
  1369.                 result = Math.pow(x, a) / a;
  1370.             }
  1371.             result *= 1 - a * x / (a + 1);
  1372.             break;
  1373.         case 7:
  1374.         default:
  1375.             // x is large,
  1376.             // Compute Q:
  1377.             result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
  1378.             result /= x;
  1379.             if (result != 0) {
  1380.                 result *= incompleteTgammaLargeX(a, x, pol);
  1381.             }
  1382.             break;
  1383.         }

  1384.         if (normalised && (result > 1)) {
  1385.             result = 1;
  1386.         }
  1387.         if (invert) {
  1388.             final double gam = normalised ? 1 : tgamma(a);
  1389.             result = gam - result;
  1390.         }

  1391.         return result;
  1392.     }

  1393.     /**
  1394.      * Upper gamma fraction.
  1395.      * Multiply result by z^a * e^-z to get the full
  1396.      * upper incomplete integral.  Divide by tgamma(z)
  1397.      * to normalise.
  1398.      *
  1399.      * @param a Argument a
  1400.      * @param z Argument z
  1401.      * @param pol Function evaluation policy
  1402.      * @return upper gamma fraction
  1403.      */
  1404.     // This is package-private for testing
  1405.     static double upperGammaFraction(double a, double z, Policy pol) {
  1406.         final double eps = pol.getEps();
  1407.         final int maxIterations = pol.getMaxIterations();

  1408.         // This is computing:
  1409.         //              1
  1410.         // ------------------------------
  1411.         // b0 + a1 / (b1 +     a2       )
  1412.         //                 -------------
  1413.         //                 b2 +    a3
  1414.         //                      --------
  1415.         //                      b3 + ...
  1416.         //
  1417.         // b0 = z - a + 1
  1418.         // a1 = a - 1
  1419.         //
  1420.         // It can be done several ways with variations in accuracy.
  1421.         // The current implementation has the best accuracy and matches the Boost code.

  1422.         final double zma1 = z - a + 1;

  1423.         final Supplier<Coefficient> gen = new Supplier<Coefficient>() {
  1424.             /** Iteration. */
  1425.             private int k;

  1426.             @Override
  1427.             public Coefficient get() {
  1428.                 ++k;
  1429.                 return Coefficient.of(k * (a - k), zma1 + 2.0 * k);
  1430.             }
  1431.         };

  1432.         return 1 / GeneralizedContinuedFraction.value(zma1, gen, eps, maxIterations);
  1433.     }

  1434.     /**
  1435.      * Upper gamma fraction for integer a.
  1436.      * Called when {@code a < 30} and {@code -x > LOG_MIN_VALUE}.
  1437.      *
  1438.      * @param a Argument a (assumed to be small)
  1439.      * @param x Argument x
  1440.      * @return upper gamma fraction
  1441.      */
  1442.     private static double finiteGammaQ(double a, double x) {
  1443.         //
  1444.         // Calculates normalised Q when a is an integer:
  1445.         //

  1446.         // Update:
  1447.         // Assume -x > log min value and no underflow to zero.

  1448.         double sum = Math.exp(-x);
  1449.         double term = sum;
  1450.         for (int n = 1; n < a; ++n) {
  1451.             term /= n;
  1452.             term *= x;
  1453.             sum += term;
  1454.         }
  1455.         return sum;
  1456.     }

  1457.     /**
  1458.      * Upper gamma fraction for half integer a.
  1459.      * Called when {@code a < 30} and {@code -x > LOG_MIN_VALUE}.
  1460.      *
  1461.      * @param a Argument a (assumed to be small)
  1462.      * @param x Argument x
  1463.      * @return upper gamma fraction
  1464.      */
  1465.     private static double finiteHalfGammaQ(double a, double x) {
  1466.         //
  1467.         // Calculates normalised Q when a is a half-integer:
  1468.         //

  1469.         // Update:
  1470.         // Assume -x > log min value:
  1471.         // erfc(sqrt(708)) = erfc(26.6) => erfc has a non-zero value

  1472.         double e = BoostErf.erfc(Math.sqrt(x));
  1473.         if (a > 1) {
  1474.             double term = Math.exp(-x) / Math.sqrt(Math.PI * x);
  1475.             term *= x;
  1476.             term /= 0.5;
  1477.             double sum = term;
  1478.             for (int n = 2; n < a; ++n) {
  1479.                 term /= n - 0.5;
  1480.                 term *= x;
  1481.                 sum += term;
  1482.             }
  1483.             e += sum;
  1484.         }
  1485.         return e;
  1486.     }

  1487.     /**
  1488.      * Lower gamma series.
  1489.      * Multiply result by ((z^a) * (e^-z) / a) to get the full
  1490.      * lower incomplete integral. Then divide by tgamma(a)
  1491.      * to get the normalised value.
  1492.      *
  1493.      * @param a Argument a
  1494.      * @param z Argument z
  1495.      * @param initValue Initial value
  1496.      * @param pol Function evaluation policy
  1497.      * @return lower gamma series
  1498.      */
  1499.     // This is package-private for testing
  1500.     static double lowerGammaSeries(double a, double z, double initValue, Policy pol) {
  1501.         final double eps = pol.getEps();
  1502.         final int maxIterations = pol.getMaxIterations();

  1503.         // Lower gamma series representation.
  1504.         final DoubleSupplier gen = new DoubleSupplier() {
  1505.             /** Next result. */
  1506.             private double result = 1;
  1507.             /** Iteration. */
  1508.             private int n;

  1509.             @Override
  1510.             public double getAsDouble() {
  1511.                 final double r = result;
  1512.                 n++;
  1513.                 result *= z / (a + n);
  1514.                 return r;
  1515.             }
  1516.         };

  1517.         return BoostTools.kahanSumSeries(gen, eps, maxIterations, initValue);
  1518.     }

  1519.     /**
  1520.      * Upper gamma fraction for very small a.
  1521.      *
  1522.      * @param a Argument a (assumed to be small)
  1523.      * @param x Argument x
  1524.      * @param pol Function evaluation policy
  1525.      * @param pgam set to value of gamma(a) on output
  1526.      * @param invert true to invert the result
  1527.      * @return upper gamma fraction
  1528.      */
  1529.     private static double tgammaSmallUpperPart(double a, double x, Policy pol, double[] pgam, boolean invert) {
  1530.         //
  1531.         // Compute the full upper fraction (Q) when a is very small:
  1532.         //
  1533.         double result;
  1534.         result = tgamma1pm1(a);

  1535.         // Note: Replacing this with tgamma(a) does not reduce error on current test data.

  1536.         // gamma(1+z) = z gamma(z)
  1537.         // pgam[0] == gamma(a)
  1538.         pgam[0] = (result + 1) / a;

  1539.         double p = BoostMath.powm1(x, a);
  1540.         result -= p;
  1541.         result /= a;
  1542.         // Removed subtraction of 10 from this value
  1543.         final int maxIter = pol.getMaxIterations();
  1544.         p += 1;
  1545.         final double initValue = invert ? pgam[0] : 0;

  1546.         // Series representation for upper fraction when z is small.
  1547.         final DoubleSupplier gen = new DoubleSupplier() {
  1548.             /** Result term. */
  1549.             private double result = -x;
  1550.             /** Argument x (this is negated on purpose). */
  1551.             private final double z = -x;
  1552.             /** Iteration. */
  1553.             private int n = 1;

  1554.             @Override
  1555.             public double getAsDouble() {
  1556.                 final double r = result / (a + n);
  1557.                 n++;
  1558.                 result = result * z / n;
  1559.                 return r;
  1560.             }
  1561.         };

  1562.         result = -p * BoostTools.kahanSumSeries(gen, pol.getEps(), maxIter, (initValue - result) / p);
  1563.         if (invert) {
  1564.             result = -result;
  1565.         }
  1566.         return result;
  1567.     }

  1568.     /**
  1569.      * Calculate power term prefix (z^a)(e^-z) used in the non-normalised
  1570.      * incomplete gammas.
  1571.      *
  1572.      * @param a Argument a
  1573.      * @param z Argument z
  1574.      * @return incomplete gamma prefix
  1575.      */
  1576.     static double fullIgammaPrefix(double a, double z) {
  1577.         if (z > Double.MAX_VALUE) {
  1578.             return 0;
  1579.         }
  1580.         final double alz = a * Math.log(z);
  1581.         final double prefix;

  1582.         if (z >= 1) {
  1583.             if ((alz < LOG_MAX_VALUE) && (-z > LOG_MIN_VALUE)) {
  1584.                 prefix = Math.pow(z, a) * Math.exp(-z);
  1585.             } else if (a >= 1) {
  1586.                 prefix = Math.pow(z / Math.exp(z / a), a);
  1587.             } else {
  1588.                 prefix = Math.exp(alz - z);
  1589.             }
  1590.         } else {
  1591.             if (alz > LOG_MIN_VALUE) {
  1592.                 prefix = Math.pow(z, a) * Math.exp(-z);
  1593.             } else {
  1594.                 // Updated to remove unreachable final branch using Math.exp(alz - z).
  1595.                 // This branch requires (z / a < LOG_MAX_VALUE) to avoid overflow in exp.
  1596.                 // At this point:
  1597.                 // 1. log(z) is negative;
  1598.                 // 2. a * log(z) <= -708 requires a > -708 / log(z).
  1599.                 // For any z < 1: -708 / log(z) is > z. Thus a is > z and
  1600.                 // z / a < LOG_MAX_VALUE is always true.
  1601.                 prefix = Math.pow(z / Math.exp(z / a), a);
  1602.             }
  1603.         }
  1604.         // Removed overflow check. Return infinity if it occurs.
  1605.         return prefix;
  1606.     }

  1607.     /**
  1608.      * Compute (z^a)(e^-z)/tgamma(a).
  1609.      * <p>Most of the error occurs in this function.
  1610.      *
  1611.      * @param a Argument a
  1612.      * @param z Argument z
  1613.      * @return regularized gamma prefix
  1614.      */
  1615.     // This is package-private for testing
  1616.     static double regularisedGammaPrefix(double a, double z) {
  1617.         if (z >= Double.MAX_VALUE) {
  1618.             return 0;
  1619.         }

  1620.         // Update this condition from: a < 1
  1621.         if (a <= 1) {
  1622.             //
  1623.             // We have to treat a < 1 as a special case because our Lanczos
  1624.             // approximations are optimised against the factorials with a > 1,
  1625.             // and for high precision types especially (128-bit reals for example)
  1626.             // very small values of a can give rather erroneous results for gamma
  1627.             // unless we do this:
  1628.             //
  1629.             // Boost todo: is this still required? Lanczos approx should be better now?
  1630.             //

  1631.             // Update this condition from: z <= LOG_MIN_VALUE
  1632.             // Require exp(-z) to not underflow:
  1633.             // -z > log(min_value)
  1634.             if (-z <= LOG_MIN_VALUE) {
  1635.                 // Oh dear, have to use logs, should be free of cancellation errors though:
  1636.                 return Math.exp(a * Math.log(z) - z - lgamma(a));
  1637.             }
  1638.             // direct calculation, no danger of overflow as gamma(a) < 1/a
  1639.             // for small a.
  1640.             return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
  1641.         }

  1642.         // Update to the Boost code.
  1643.         // Use some of the logic from fullIgammaPrefix(a, z) to use the direct
  1644.         // computation if it is valid. Assuming pow and exp are accurate to 1 ULP it
  1645.         // puts most of the error in evaluation of tgamma(a). This is accurate
  1646.         // enough that this reduces max error on the current test data.
  1647.         //
  1648.         // Overflow cases fall-through to the Lanczos approximation that incorporates
  1649.         // the pow and exp terms used in the tgamma(a) computation with the terms
  1650.         // z^a and e^-z into a single evaluation of pow and exp. See equation 15:
  1651.         // https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html
  1652.         if (a <= MAX_GAMMA_Z) {
  1653.             final double alz1 = a * Math.log(z);
  1654.             if (z >= 1) {
  1655.                 if ((alz1 < LOG_MAX_VALUE) && (-z > LOG_MIN_VALUE)) {
  1656.                     return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
  1657.                 }
  1658.             } else if (alz1 > LOG_MIN_VALUE) {
  1659.                 return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
  1660.             }
  1661.         }

  1662.         //
  1663.         // For smallish a and x combining the power terms with the Lanczos approximation
  1664.         // gives the greatest accuracy
  1665.         //

  1666.         final double agh = a + Lanczos.GMH;
  1667.         double prefix;

  1668.         final double factor = Math.sqrt(agh / Math.E) / Lanczos.lanczosSumExpGScaled(a);

  1669.         // Update to the Boost code.
  1670.         // Lower threshold for large a from 150 to 128 and compute d on demand.
  1671.         // See NUMBERS-179.
  1672.         if (a > 128) {
  1673.             final double d = ((z - a) - Lanczos.GMH) / agh;
  1674.             if (Math.abs(d * d * a) <= 100) {
  1675.                 // special case for large a and a ~ z.
  1676.                 // When a and x are large, we end up with a very large exponent with a base near one:
  1677.                 // this will not be computed accurately via the pow function, and taking logs simply
  1678.                 // leads to cancellation errors.
  1679.                 prefix = a * SpecialMath.log1pmx(d) + z * -Lanczos.GMH / agh;
  1680.                 prefix = Math.exp(prefix);
  1681.                 return prefix * factor;
  1682.             }
  1683.         }

  1684.         //
  1685.         // general case.
  1686.         // direct computation is most accurate, but use various fallbacks
  1687.         // for different parts of the problem domain:
  1688.         //

  1689.         final double alz = a * Math.log(z / agh);
  1690.         final double amz = a - z;
  1691.         if ((Math.min(alz, amz) <= LOG_MIN_VALUE) || (Math.max(alz, amz) >= LOG_MAX_VALUE)) {
  1692.             final double amza = amz / a;
  1693.             if ((Math.min(alz, amz) / 2 > LOG_MIN_VALUE) && (Math.max(alz, amz) / 2 < LOG_MAX_VALUE)) {
  1694.                 // compute square root of the result and then square it:
  1695.                 final double sq = Math.pow(z / agh, a / 2) * Math.exp(amz / 2);
  1696.                 prefix = sq * sq;
  1697.             } else if ((Math.min(alz, amz) / 4 > LOG_MIN_VALUE) &&
  1698.                     (Math.max(alz, amz) / 4 < LOG_MAX_VALUE) && (z > a)) {
  1699.                 // compute the 4th root of the result then square it twice:
  1700.                 final double sq = Math.pow(z / agh, a / 4) * Math.exp(amz / 4);
  1701.                 prefix = sq * sq;
  1702.                 prefix *= prefix;
  1703.             } else if ((amza > LOG_MIN_VALUE) && (amza < LOG_MAX_VALUE)) {
  1704.                 prefix = Math.pow((z * Math.exp(amza)) / agh, a);
  1705.             } else {
  1706.                 prefix = Math.exp(alz + amz);
  1707.             }
  1708.         } else {
  1709.             prefix = Math.pow(z / agh, a) * Math.exp(amz);
  1710.         }
  1711.         prefix *= factor;
  1712.         return prefix;
  1713.     }

  1714.     /**
  1715.      * Implements the asymptotic expansions of the incomplete
  1716.      * gamma functions P(a, x) and Q(a, x), used when a is large and
  1717.      * x ~ a.
  1718.      *
  1719.      * <p>The primary reference is:
  1720.      * <pre>
  1721.      * "The Asymptotic Expansion of the Incomplete Gamma Functions"
  1722.      * N. M. Temme.
  1723.      * Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
  1724.      * </pre>
  1725.      *
  1726.      * <p>A different way of evaluating these expansions,
  1727.      * plus a lot of very useful background information is in:
  1728.      * <pre>
  1729.      * "A Set of Algorithms For the Incomplete Gamma Functions."
  1730.      * N. M. Temme.
  1731.      * Probability in the Engineering and Informational Sciences,
  1732.      * 8, 1994, 291.
  1733.      * </pre>
  1734.      *
  1735.      * <p>An alternative implementation is in:
  1736.      * <pre>
  1737.      * "Computation of the Incomplete Gamma Function Ratios and their Inverse."
  1738.      * A. R. Didonato and A. H. Morris.
  1739.      * ACM TOMS, Vol 12, No 4, Dec 1986, p377.
  1740.      * </pre>
  1741.      *
  1742.      * <p>This is a port of the function accurate for 53-bit mantissas
  1743.      * (IEEE double precision or 10^-17). To understand the code, refer to Didonato
  1744.      * and Morris, from Eq 17 and 18 onwards.
  1745.      *
  1746.      * <p>The coefficients used here are not taken from Didonato and Morris:
  1747.      * the domain over which these expansions are used is slightly different
  1748.      * to theirs, and their constants are not quite accurate enough for
  1749.      * 128-bit long doubles.  Instead the coefficients were calculated
  1750.      * using the methods described by Temme p762 from Eq 3.8 onwards.
  1751.      * The values obtained agree with those obtained by Didonato and Morris
  1752.      * (at least to the first 30 digits that they provide).
  1753.      * At double precision the degrees of polynomial required for full
  1754.      * machine precision are close to those recommended to Didonato and Morris,
  1755.      * but of course many more terms are needed for larger types.
  1756.      *
  1757.      * <p>Adapted from {@code boost/math/special_functions/detail/igamma_large.hpp}.
  1758.      *
  1759.      * @param a the a
  1760.      * @param x the x
  1761.      * @return the double
  1762.      */
  1763.     // This is package-private for testing
  1764.     static double igammaTemmeLarge(double a, double x) {
  1765.         final double sigma = (x - a) / a;
  1766.         final double phi = -SpecialMath.log1pmx(sigma);
  1767.         final double y = a * phi;
  1768.         double z = Math.sqrt(2 * phi);
  1769.         if (x < a) {
  1770.             z = -z;
  1771.         }

  1772.         // The following polynomials are evaluated with a loop
  1773.         // with Horner's method. Variations exist using
  1774.         // a second order Horner's method with an unrolled loop.
  1775.         // These are chosen in Boost based on the C++ compiler.
  1776.         // For example:
  1777.         // boost/math/tools/detail/polynomial_horner1_15.hpp
  1778.         // boost/math/tools/detail/polynomial_horner2_15.hpp
  1779.         // boost/math/tools/detail/polynomial_horner3_15.hpp

  1780.         final double[] workspace = new double[10];

  1781.         final double[] C0 = {
  1782.             -0.33333333333333333,
  1783.             0.083333333333333333,
  1784.             -0.014814814814814815,
  1785.             0.0011574074074074074,
  1786.             0.0003527336860670194,
  1787.             -0.00017875514403292181,
  1788.             0.39192631785224378e-4,
  1789.             -0.21854485106799922e-5,
  1790.             -0.185406221071516e-5,
  1791.             0.8296711340953086e-6,
  1792.             -0.17665952736826079e-6,
  1793.             0.67078535434014986e-8,
  1794.             0.10261809784240308e-7,
  1795.             -0.43820360184533532e-8,
  1796.             0.91476995822367902e-9,
  1797.         };
  1798.         workspace[0] = BoostTools.evaluatePolynomial(C0, z);

  1799.         final double[] C1 = {
  1800.             -0.0018518518518518519,
  1801.             -0.0034722222222222222,
  1802.             0.0026455026455026455,
  1803.             -0.00099022633744855967,
  1804.             0.00020576131687242798,
  1805.             -0.40187757201646091e-6,
  1806.             -0.18098550334489978e-4,
  1807.             0.76491609160811101e-5,
  1808.             -0.16120900894563446e-5,
  1809.             0.46471278028074343e-8,
  1810.             0.1378633446915721e-6,
  1811.             -0.5752545603517705e-7,
  1812.             0.11951628599778147e-7,
  1813.         };
  1814.         workspace[1] = BoostTools.evaluatePolynomial(C1, z);

  1815.         final double[] C2 = {
  1816.             0.0041335978835978836,
  1817.             -0.0026813271604938272,
  1818.             0.00077160493827160494,
  1819.             0.20093878600823045e-5,
  1820.             -0.00010736653226365161,
  1821.             0.52923448829120125e-4,
  1822.             -0.12760635188618728e-4,
  1823.             0.34235787340961381e-7,
  1824.             0.13721957309062933e-5,
  1825.             -0.6298992138380055e-6,
  1826.             0.14280614206064242e-6,
  1827.         };
  1828.         workspace[2] = BoostTools.evaluatePolynomial(C2, z);

  1829.         final double[] C3 = {
  1830.             0.00064943415637860082,
  1831.             0.00022947209362139918,
  1832.             -0.00046918949439525571,
  1833.             0.00026772063206283885,
  1834.             -0.75618016718839764e-4,
  1835.             -0.23965051138672967e-6,
  1836.             0.11082654115347302e-4,
  1837.             -0.56749528269915966e-5,
  1838.             0.14230900732435884e-5,
  1839.         };
  1840.         workspace[3] = BoostTools.evaluatePolynomial(C3, z);

  1841.         final double[] C4 = {
  1842.             -0.0008618882909167117,
  1843.             0.00078403922172006663,
  1844.             -0.00029907248030319018,
  1845.             -0.14638452578843418e-5,
  1846.             0.66414982154651222e-4,
  1847.             -0.39683650471794347e-4,
  1848.             0.11375726970678419e-4,
  1849.         };
  1850.         workspace[4] = BoostTools.evaluatePolynomial(C4, z);

  1851.         final double[] C5 = {
  1852.             -0.00033679855336635815,
  1853.             -0.69728137583658578e-4,
  1854.             0.00027727532449593921,
  1855.             -0.00019932570516188848,
  1856.             0.67977804779372078e-4,
  1857.             0.1419062920643967e-6,
  1858.             -0.13594048189768693e-4,
  1859.             0.80184702563342015e-5,
  1860.             -0.22914811765080952e-5,
  1861.         };
  1862.         workspace[5] = BoostTools.evaluatePolynomial(C5, z);

  1863.         final double[] C6 = {
  1864.             0.00053130793646399222,
  1865.             -0.00059216643735369388,
  1866.             0.00027087820967180448,
  1867.             0.79023532326603279e-6,
  1868.             -0.81539693675619688e-4,
  1869.             0.56116827531062497e-4,
  1870.             -0.18329116582843376e-4,
  1871.         };
  1872.         workspace[6] = BoostTools.evaluatePolynomial(C6, z);

  1873.         final double[] C7 = {
  1874.             0.00034436760689237767,
  1875.             0.51717909082605922e-4,
  1876.             -0.00033493161081142236,
  1877.             0.0002812695154763237,
  1878.             -0.00010976582244684731,
  1879.         };
  1880.         workspace[7] = BoostTools.evaluatePolynomial(C7, z);

  1881.         final double[] C8 = {
  1882.             -0.00065262391859530942,
  1883.             0.00083949872067208728,
  1884.             -0.00043829709854172101,
  1885.         };
  1886.         workspace[8] = BoostTools.evaluatePolynomial(C8, z);
  1887.         workspace[9] = -0.00059676129019274625;

  1888.         double result = BoostTools.evaluatePolynomial(workspace, 1 / a);
  1889.         result *= Math.exp(-y) / Math.sqrt(2 * Math.PI * a);
  1890.         if (x < a) {
  1891.             result = -result;
  1892.         }

  1893.         result += BoostErf.erfc(Math.sqrt(y)) / 2;

  1894.         return result;
  1895.     }

  1896.     /**
  1897.      * Incomplete tgamma for large X.
  1898.      *
  1899.      * <p>This summation is a source of error as the series starts at 1 and descends to zero.
  1900.      * It can have thousands of iterations when a and z are large and close in value.
  1901.      *
  1902.      * @param a Argument a
  1903.      * @param x Argument x
  1904.      * @param pol Function evaluation policy
  1905.      * @return incomplete tgamma
  1906.      */
  1907.     // This is package-private for testing
  1908.     static double incompleteTgammaLargeX(double a, double x, Policy pol) {
  1909.         final double eps = pol.getEps();
  1910.         final int maxIterations = pol.getMaxIterations();

  1911.         // Asymptotic approximation for large argument, see: https://dlmf.nist.gov/8.11#E2.
  1912.         final DoubleSupplier gen = new DoubleSupplier() {
  1913.             /** Result term. */
  1914.             private double term = 1;
  1915.             /** Iteration. */
  1916.             private int n;

  1917.             @Override
  1918.             public double getAsDouble() {
  1919.                 final double result = term;
  1920.                 n++;
  1921.                 term *= (a - n) / x;
  1922.                 return result;
  1923.             }
  1924.         };

  1925.         return BoostTools.kahanSumSeries(gen, eps, maxIterations);
  1926.     }

  1927.     /**
  1928.      * Return true if the array is not null and has non-zero length.
  1929.      *
  1930.      * @param array Array
  1931.      * @return true if a non-zero length array
  1932.      */
  1933.     private static boolean nonZeroLength(int[] array) {
  1934.         return array != null && array.length != 0;
  1935.     }

  1936.     /**
  1937.      * Ratio of gamma functions.
  1938.      *
  1939.      * <p>\[ tgamma_ratio(z, delta) = \frac{\Gamma(z)}{\Gamma(z + delta)} \]
  1940.      *
  1941.      * <p>Adapted from {@code tgamma_delta_ratio_imp}. The use of
  1942.      * {@code max_factorial<double>::value == 170} has been replaced with
  1943.      * {@code MAX_GAMMA_Z == 171}. This threshold is used when it is possible
  1944.      * to call the gamma function without overflow.
  1945.      *
  1946.      * @param z Argument z
  1947.      * @param delta The difference
  1948.      * @return gamma ratio
  1949.      */
  1950.     static double tgammaDeltaRatio(double z, double delta) {
  1951.         final double zDelta = z + delta;
  1952.         if (Double.isNaN(zDelta)) {
  1953.             // One or both arguments are NaN
  1954.             return Double.NaN;
  1955.         }
  1956.         if (z <= 0 || zDelta <= 0) {
  1957.             // This isn't very sophisticated, or accurate, but it does work:
  1958.             return tgamma(z) / tgamma(zDelta);
  1959.         }

  1960.         // Note: Directly calling tgamma(z) / tgamma(z + delta) if possible
  1961.         // without overflow is not more accurate

  1962.         if (Math.rint(delta) == delta) {
  1963.             if (delta == 0) {
  1964.                 return 1;
  1965.             }
  1966.             //
  1967.             // If both z and delta are integers, see if we can just use table lookup
  1968.             // of the factorials to get the result:
  1969.             //
  1970.             if (Math.rint(z) == z &&
  1971.                 z <= MAX_GAMMA_Z && zDelta <= MAX_GAMMA_Z) {
  1972.                 return FACTORIAL[(int) z - 1] / FACTORIAL[(int) zDelta - 1];
  1973.             }
  1974.             if (Math.abs(delta) < 20) {
  1975.                 //
  1976.                 // delta is a small integer, we can use a finite product:
  1977.                 //
  1978.                 if (delta < 0) {
  1979.                     z -= 1;
  1980.                     double result = z;
  1981.                     for (int d = (int) (delta + 1); d != 0; d++) {
  1982.                         z -= 1;
  1983.                         result *= z;
  1984.                     }
  1985.                     return result;
  1986.                 }
  1987.                 double result = 1 / z;
  1988.                 for (int d = (int) (delta - 1); d != 0; d--) {
  1989.                     z += 1;
  1990.                     result /= z;
  1991.                 }
  1992.                 return result;
  1993.             }
  1994.         }
  1995.         return tgammaDeltaRatioImpLanczos(z, delta);
  1996.     }

  1997.     /**
  1998.      * Ratio of gamma functions using Lanczos support.
  1999.      *
  2000.      * <p>\[ tgamma_delta_ratio(z, delta) = \frac{\Gamma(z)}{\Gamma(z + delta)}
  2001.      *
  2002.      * <p>Adapted from {@code tgamma_delta_ratio_imp_lanczos}. The use of
  2003.      * {@code max_factorial<double>::value == 170} has been replaced with
  2004.      * {@code MAX_GAMMA_Z == 171}. This threshold is used when it is possible
  2005.      * to use the precomputed factorial table.
  2006.      *
  2007.      * @param z Argument z
  2008.      * @param delta The difference
  2009.      * @return gamma ratio
  2010.      */
  2011.     private static double tgammaDeltaRatioImpLanczos(double z, double delta) {
  2012.         if (z < EPSILON) {
  2013.             //
  2014.             // We get spurious numeric overflow unless we're very careful, this
  2015.             // can occur either inside Lanczos::lanczos_sum(z) or in the
  2016.             // final combination of terms, to avoid this, split the product up
  2017.             // into 2 (or 3) parts:
  2018.             //
  2019.             // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta
  2020.             // z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial
  2021.             //
  2022.             if (MAX_GAMMA_Z < delta) {
  2023.                 double ratio = tgammaDeltaRatioImpLanczos(delta, MAX_GAMMA_Z - delta);
  2024.                 ratio *= z;
  2025.                 ratio *= FACTORIAL[MAX_FACTORIAL];
  2026.                 return 1 / ratio;
  2027.             }
  2028.             return 1 / (z * tgamma(z + delta));
  2029.         }
  2030.         final double zgh = z + Lanczos.GMH;
  2031.         double result;
  2032.         if (z + delta == z) {
  2033.             // Here:
  2034.             // lanczosSum(z) / lanczosSum(z + delta) == 1

  2035.             // Update to the Boost code to remove unreachable code:
  2036.             // Given z + delta == z then |delta / z| < EPSILON; and z < zgh
  2037.             // assume (abs(delta / zgh) < EPSILON) and remove unreachable
  2038.             // else branch which sets result = 1

  2039.             // We have:
  2040.             // result = exp((0.5 - z) * log1p(delta / zgh));
  2041.             // 0.5 - z == -z
  2042.             // log1p(delta / zgh) = delta / zgh = delta / z
  2043.             // multiplying we get -delta.

  2044.             // Note:
  2045.             // This can be different from
  2046.             // exp((0.5 - z) * log1p(delta / zgh)) when z is small.
  2047.             // In this case the result is exp(small) and the result
  2048.             // is within 1 ULP of 1.0. This is left as the original
  2049.             // Boost method using exp(-delta).

  2050.             result = Math.exp(-delta);
  2051.         } else {
  2052.             if (Math.abs(delta) < 10) {
  2053.                 result = Math.exp((0.5 - z) * Math.log1p(delta / zgh));
  2054.             } else {
  2055.                 result = Math.pow(zgh / (zgh + delta), z - 0.5);
  2056.             }
  2057.             // Split the calculation up to avoid spurious overflow:
  2058.             result *= Lanczos.lanczosSum(z) / Lanczos.lanczosSum(z + delta);
  2059.         }
  2060.         result *= Math.pow(Math.E / (zgh + delta), delta);
  2061.         return result;
  2062.     }

  2063.     /**
  2064.      * Ratio of gamma functions.
  2065.      *
  2066.      * <p>\[ tgamma_ratio(x, y) = \frac{\Gamma(x)}{\Gamma(y)}
  2067.      *
  2068.      * <p>Adapted from {@code tgamma_ratio_imp}. The use of
  2069.      * {@code max_factorial<double>::value == 170} has been replaced with
  2070.      * {@code MAX_GAMMA_Z == 171}. This threshold is used when it is possible
  2071.      * to call the gamma function without overflow.
  2072.      *
  2073.      * @param x Argument x (must be positive finite)
  2074.      * @param y Argument y (must be positive finite)
  2075.      * @return gamma ratio (or nan)
  2076.      */
  2077.     static double tgammaRatio(double x, double y) {
  2078.         if (x <= 0 || !Double.isFinite(x) || y <= 0 || !Double.isFinite(y)) {
  2079.             return Double.NaN;
  2080.         }
  2081.         if (x <= Double.MIN_NORMAL) {
  2082.             // Special case for denorms...Ugh.
  2083.             return TWO_POW_53 * tgammaRatio(x * TWO_POW_53, y);
  2084.         }

  2085.         if (x <= MAX_GAMMA_Z && y <= MAX_GAMMA_Z) {
  2086.             // Rather than subtracting values, lets just call the gamma functions directly:
  2087.             return tgamma(x) / tgamma(y);
  2088.         }
  2089.         double prefix = 1;
  2090.         if (x < 1) {
  2091.             if (y < 2 * MAX_GAMMA_Z) {
  2092.                 // We need to sidestep on x as well, otherwise we'll underflow
  2093.                 // before we get to factor in the prefix term:
  2094.                 prefix /= x;
  2095.                 x += 1;
  2096.                 while (y >= MAX_GAMMA_Z) {
  2097.                     y -= 1;
  2098.                     prefix /= y;
  2099.                 }
  2100.                 return prefix * tgamma(x) / tgamma(y);
  2101.             }
  2102.             //
  2103.             // result is almost certainly going to underflow to zero, try logs just in case:
  2104.             //
  2105.             return Math.exp(lgamma(x) - lgamma(y));
  2106.         }
  2107.         if (y < 1) {
  2108.             if (x < 2 * MAX_GAMMA_Z) {
  2109.                 // We need to sidestep on y as well, otherwise we'll overflow
  2110.                 // before we get to factor in the prefix term:
  2111.                 prefix *= y;
  2112.                 y += 1;
  2113.                 while (x >= MAX_GAMMA_Z) {
  2114.                     x -= 1;
  2115.                     prefix *= x;
  2116.                 }
  2117.                 return prefix * tgamma(x) / tgamma(y);
  2118.             }
  2119.             //
  2120.             // Result will almost certainly overflow, try logs just in case:
  2121.             //
  2122.             return Math.exp(lgamma(x) - lgamma(y));
  2123.         }
  2124.         //
  2125.         // Regular case, x and y both large and similar in magnitude:
  2126.         //
  2127.         return tgammaDeltaRatio(x, y - x);
  2128.     }
  2129. }