BoostGamma.java
- /*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
- // Copyright John Maddock 2006-7, 2013-20.
- // Copyright Paul A. Bristow 2007, 2013-14.
- // Copyright Nikhar Agrawal 2013-14
- // Copyright Christopher Kormanyos 2013-14, 2020
- // Use, modification and distribution are subject to the
- // Boost Software License, Version 1.0. (See accompanying file
- // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
- package org.apache.commons.numbers.gamma;
- import java.util.function.DoubleSupplier;
- import java.util.function.Supplier;
- import org.apache.commons.numbers.core.Sum;
- import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction;
- import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction.Coefficient;
- /**
- * Implementation of the
- * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">Regularized Gamma functions</a> and
- * <a href="https://mathworld.wolfram.com/IncompleteGammaFunction.html">Incomplete Gamma functions</a>.
- *
- * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
- * {@code c++} implementation {@code <boost/math/special_functions/gamma.hpp>}.
- * All work is copyright to the original authors and subject to the Boost Software License.
- *
- * @see
- * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma.html">
- * Boost C++ Gamma functions</a>
- */
- final class BoostGamma {
- //
- // Code ported from Boost 1.77.0
- //
- // boost/math/special_functions/gamma.hpp
- // boost/math/special_functions/detail/igamma_large.hpp
- // boost/math/special_functions/lanczos.hpp
- //
- // Original code comments are preserved.
- //
- // Changes to the Boost implementation:
- // - Update method names to replace underscores with camel case
- // - Explicitly inline the polynomial function evaluation
- // using Horner's method (https://en.wikipedia.org/wiki/Horner%27s_method)
- // - Remove checks for under/overflow. In this implementation no error is raised
- // for overflow (infinity is returned) or underflow (sub-normal or zero is returned).
- // This follows the conventions in java.lang.Math for the same conditions.
- // - Removed the pointer p_derivative in the gammaIncompleteImp. This is used
- // in the Boost code for the gamma_(p|q)_inv functions for a derivative
- // based inverse function. This is currently not supported.
- // - Added extended precision arithmetic for some series summations or other computations.
- // The Boost default policy is to evaluate in long double for a double result. Extended
- // precision is not possible for the entire computation but has been used where
- // possible for some terms to reduce errors. Error reduction verified on the test data.
- // - Altered the tgamma(x) function to use the double-precision NSWC Library of Mathematics
- // Subroutines when the error is lower. This is for the non Lanczos code.
- // - Altered the condition used for the asymptotic approximation method to avoid
- // loss of precision in the series summation when a ~ z.
- // - Altered series generators to use integer counters added to the double term
- // replacing directly incrementing a double term. When the term is large it cannot
- // be incremented: 1e16 + 1 == 1e16.
- // - Removed unreachable code branch in tgammaDeltaRatioImpLanczos when z + delta == z.
- //
- // Note:
- // The major source of error is in the function regularisedGammaPrefix when computing
- // (z^a)(e^-z)/tgamma(a) with extreme input to the power and exponential terms.
- // An extended precision pow and exp function returning a quad length result would
- // be required to reduce error for these arguments. Tests using the Dfp class
- // from o.a.c.math4.legacy.core have been demonstrated to effectively eliminate the
- // errors from the power terms and improve accuracy on the current test data.
- // In the interest of performance the Dfp class is not used in this version.
- /** Default epsilon value for relative error.
- * This is equal to the Boost constant {@code boost::math::tools::epsilon<double>()}. */
- private static final double EPSILON = 0x1.0p-52;
- /** Value for the sqrt of the epsilon for relative error.
- * This is equal to the Boost constant {@code boost::math::tools::root_epsilon<double>()}. */
- private static final double ROOT_EPSILON = 1.4901161193847656E-8;
- /** Approximate value for ln(Double.MAX_VALUE).
- * This is equal to the Boost constant {@code boost::math::tools::log_max_value<double>()}.
- * No term {@code x} should be used in {@code exp(x)} if {@code x > LOG_MAX_VALUE} to avoid
- * overflow. */
- private static final int LOG_MAX_VALUE = 709;
- /** Approximate value for ln(Double.MIN_VALUE).
- * This is equal to the Boost constant {@code boost::math::tools::log_min_value<double>()}.
- * No term {@code x} should be used in {@code exp(x)} if {@code x < LOG_MIN_VALUE} to avoid
- * underflow to sub-normal or zero. */
- private static final int LOG_MIN_VALUE = -708;
- /** The largest factorial that can be represented as a double.
- * This is equal to the Boost constant {@code boost::math::max_factorial<double>::value}. */
- private static final int MAX_FACTORIAL = 170;
- /** The largest integer value for gamma(z) that can be represented as a double. */
- private static final int MAX_GAMMA_Z = MAX_FACTORIAL + 1;
- /** ln(sqrt(2 pi)). Computed to 25-digits precision. */
- private static final double LOG_ROOT_TWO_PI = 0.9189385332046727417803297;
- /** ln(pi). Computed to 25-digits precision. */
- private static final double LOG_PI = 1.144729885849400174143427;
- /** Euler's constant. */
- private static final double EULER = 0.5772156649015328606065120900824024310;
- /** The threshold value for choosing the Lanczos approximation. */
- private static final int LANCZOS_THRESHOLD = 20;
- /** 2^53. */
- private static final double TWO_POW_53 = 0x1.0p53;
- /** All factorials that can be represented as a double. Size = 171. */
- private static final double[] FACTORIAL = {
- 1,
- 1,
- 2,
- 6,
- 24,
- 120,
- 720,
- 5040,
- 40320,
- 362880.0,
- 3628800.0,
- 39916800.0,
- 479001600.0,
- 6227020800.0,
- 87178291200.0,
- 1307674368000.0,
- 20922789888000.0,
- 355687428096000.0,
- 6402373705728000.0,
- 121645100408832000.0,
- 0.243290200817664e19,
- 0.5109094217170944e20,
- 0.112400072777760768e22,
- 0.2585201673888497664e23,
- 0.62044840173323943936e24,
- 0.15511210043330985984e26,
- 0.403291461126605635584e27,
- 0.10888869450418352160768e29,
- 0.304888344611713860501504e30,
- 0.8841761993739701954543616e31,
- 0.26525285981219105863630848e33,
- 0.822283865417792281772556288e34,
- 0.26313083693369353016721801216e36,
- 0.868331761881188649551819440128e37,
- 0.29523279903960414084761860964352e39,
- 0.103331479663861449296666513375232e41,
- 0.3719933267899012174679994481508352e42,
- 0.137637530912263450463159795815809024e44,
- 0.5230226174666011117600072241000742912e45,
- 0.203978820811974433586402817399028973568e47,
- 0.815915283247897734345611269596115894272e48,
- 0.3345252661316380710817006205344075166515e50,
- 0.1405006117752879898543142606244511569936e52,
- 0.6041526306337383563735513206851399750726e53,
- 0.265827157478844876804362581101461589032e55,
- 0.1196222208654801945619631614956577150644e57,
- 0.5502622159812088949850305428800254892962e58,
- 0.2586232415111681806429643551536119799692e60,
- 0.1241391559253607267086228904737337503852e62,
- 0.6082818640342675608722521633212953768876e63,
- 0.3041409320171337804361260816606476884438e65,
- 0.1551118753287382280224243016469303211063e67,
- 0.8065817517094387857166063685640376697529e68,
- 0.427488328406002556429801375338939964969e70,
- 0.2308436973392413804720927426830275810833e72,
- 0.1269640335365827592596510084756651695958e74,
- 0.7109985878048634518540456474637249497365e75,
- 0.4052691950487721675568060190543232213498e77,
- 0.2350561331282878571829474910515074683829e79,
- 0.1386831185456898357379390197203894063459e81,
- 0.8320987112741390144276341183223364380754e82,
- 0.507580213877224798800856812176625227226e84,
- 0.3146997326038793752565312235495076408801e86,
- 0.1982608315404440064116146708361898137545e88,
- 0.1268869321858841641034333893351614808029e90,
- 0.8247650592082470666723170306785496252186e91,
- 0.5443449390774430640037292402478427526443e93,
- 0.3647111091818868528824985909660546442717e95,
- 0.2480035542436830599600990418569171581047e97,
- 0.1711224524281413113724683388812728390923e99,
- 0.1197857166996989179607278372168909873646e101,
- 0.8504785885678623175211676442399260102886e102,
- 0.6123445837688608686152407038527467274078e104,
- 0.4470115461512684340891257138125051110077e106,
- 0.3307885441519386412259530282212537821457e108,
- 0.2480914081139539809194647711659403366093e110,
- 0.188549470166605025498793226086114655823e112,
- 0.1451830920282858696340707840863082849837e114,
- 0.1132428117820629783145752115873204622873e116,
- 0.8946182130782975286851441715398316520698e117,
- 0.7156945704626380229481153372318653216558e119,
- 0.5797126020747367985879734231578109105412e121,
- 0.4753643337012841748421382069894049466438e123,
- 0.3945523969720658651189747118012061057144e125,
- 0.3314240134565353266999387579130131288001e127,
- 0.2817104114380550276949479442260611594801e129,
- 0.2422709538367273238176552320344125971528e131,
- 0.210775729837952771721360051869938959523e133,
- 0.1854826422573984391147968456455462843802e135,
- 0.1650795516090846108121691926245361930984e137,
- 0.1485715964481761497309522733620825737886e139,
- 0.1352001527678402962551665687594951421476e141,
- 0.1243841405464130725547532432587355307758e143,
- 0.1156772507081641574759205162306240436215e145,
- 0.1087366156656743080273652852567866010042e147,
- 0.103299784882390592625997020993947270954e149,
- 0.9916779348709496892095714015418938011582e150,
- 0.9619275968248211985332842594956369871234e152,
- 0.942689044888324774562618574305724247381e154,
- 0.9332621544394415268169923885626670049072e156,
- 0.9332621544394415268169923885626670049072e158,
- 0.9425947759838359420851623124482936749562e160,
- 0.9614466715035126609268655586972595484554e162,
- 0.990290071648618040754671525458177334909e164,
- 0.1029901674514562762384858386476504428305e167,
- 0.1081396758240290900504101305800329649721e169,
- 0.1146280563734708354534347384148349428704e171,
- 0.1226520203196137939351751701038733888713e173,
- 0.132464181945182897449989183712183259981e175,
- 0.1443859583202493582204882102462797533793e177,
- 0.1588245541522742940425370312709077287172e179,
- 0.1762952551090244663872161047107075788761e181,
- 0.1974506857221074023536820372759924883413e183,
- 0.2231192748659813646596607021218715118256e185,
- 0.2543559733472187557120132004189335234812e187,
- 0.2925093693493015690688151804817735520034e189,
- 0.339310868445189820119825609358857320324e191,
- 0.396993716080872089540195962949863064779e193,
- 0.4684525849754290656574312362808384164393e195,
- 0.5574585761207605881323431711741977155627e197,
- 0.6689502913449127057588118054090372586753e199,
- 0.8094298525273443739681622845449350829971e201,
- 0.9875044200833601362411579871448208012564e203,
- 0.1214630436702532967576624324188129585545e206,
- 0.1506141741511140879795014161993280686076e208,
- 0.1882677176888926099743767702491600857595e210,
- 0.237217324288004688567714730513941708057e212,
- 0.3012660018457659544809977077527059692324e214,
- 0.3856204823625804217356770659234636406175e216,
- 0.4974504222477287440390234150412680963966e218,
- 0.6466855489220473672507304395536485253155e220,
- 0.8471580690878820510984568758152795681634e222,
- 0.1118248651196004307449963076076169029976e225,
- 0.1487270706090685728908450891181304809868e227,
- 0.1992942746161518876737324194182948445223e229,
- 0.269047270731805048359538766214698040105e231,
- 0.3659042881952548657689727220519893345429e233,
- 0.5012888748274991661034926292112253883237e235,
- 0.6917786472619488492228198283114910358867e237,
- 0.9615723196941089004197195613529725398826e239,
- 0.1346201247571752460587607385894161555836e242,
- 0.1898143759076170969428526414110767793728e244,
- 0.2695364137888162776588507508037290267094e246,
- 0.3854370717180072770521565736493325081944e248,
- 0.5550293832739304789551054660550388118e250,
- 0.80479260574719919448490292577980627711e252,
- 0.1174997204390910823947958271638517164581e255,
- 0.1727245890454638911203498659308620231933e257,
- 0.2556323917872865588581178015776757943262e259,
- 0.380892263763056972698595524350736933546e261,
- 0.571338395644585459047893286526105400319e263,
- 0.8627209774233240431623188626544191544816e265,
- 0.1311335885683452545606724671234717114812e268,
- 0.2006343905095682394778288746989117185662e270,
- 0.308976961384735088795856467036324046592e272,
- 0.4789142901463393876335775239063022722176e274,
- 0.7471062926282894447083809372938315446595e276,
- 0.1172956879426414428192158071551315525115e279,
- 0.1853271869493734796543609753051078529682e281,
- 0.2946702272495038326504339507351214862195e283,
- 0.4714723635992061322406943211761943779512e285,
- 0.7590705053947218729075178570936729485014e287,
- 0.1229694218739449434110178928491750176572e290,
- 0.2004401576545302577599591653441552787813e292,
- 0.3287218585534296227263330311644146572013e294,
- 0.5423910666131588774984495014212841843822e296,
- 0.9003691705778437366474261723593317460744e298,
- 0.1503616514864999040201201707840084015944e301,
- 0.2526075744973198387538018869171341146786e303,
- 0.4269068009004705274939251888899566538069e305,
- 0.7257415615307998967396728211129263114717e307,
- };
- /**
- * 53-bit precision implementation of the Lanczos approximation.
- *
- * <p>This implementation is in partial fraction form with the leading constant
- * of \( \sqrt{2\pi} \) absorbed into the sum.
- *
- * <p>It is related to the Gamma function by the following equation
- * \[
- * \Gamma(z) = \frac{(z + g - 0.5)^{z - 0.5}}{e^{z + g - 0.5}} \mathrm{lanczos}(z)
- * \]
- * where \( g \) is the Lanczos constant.
- *
- * <h2>Warning</h2>
- *
- * <p>This is not a substitute for {@link LanczosApproximation}. The approximation is
- * written in partial fraction form with the leading constants absorbed by the
- * coefficients in the sum.
- *
- * @see <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/lanczos.html">
- * Boost Lanczos Approximation</a>
- */
- static final class Lanczos {
- // Optimal values for G for each N are taken from
- // http://web.mala.bc.ca/pughg/phdThesis/phdThesis.pdf,
- // as are the theoretical error bounds.
- //
- // Constants calculated using the method described by Godfrey
- // http://my.fit.edu/~gabdo/gamma.txt and elaborated by Toth at
- // http://www.rskey.org/gamma.htm using NTL::RR at 1000 bit precision.
- //
- // Lanczos Coefficients for N=13 G=6.024680040776729583740234375
- // Max experimental error (with arbitrary precision arithmetic) 1.196214e-17
- // Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
- //
- /**
- * Lanczos constant G.
- */
- static final double G = 6.024680040776729583740234375;
- /**
- * Lanczos constant G - half.
- *
- * <p>Note: The form {@code (g - 0.5)} is used when computing the gamma function.
- */
- static final double GMH = 5.524680040776729583740234375;
- /** Common denominator used for the rational evaluation. */
- private static final int[] DENOM = {
- 0,
- 39916800,
- 120543840,
- 150917976,
- 105258076,
- 45995730,
- 13339535,
- 2637558,
- 357423,
- 32670,
- 1925,
- 66,
- 1
- };
- /** Private constructor. */
- private Lanczos() {
- // intentionally empty.
- }
- /**
- * Computes the Lanczos approximation.
- *
- * @param z Argument.
- * @return the Lanczos approximation.
- */
- static double lanczosSum(double z) {
- final double[] num = {
- 23531376880.41075968857200767445163675473,
- 42919803642.64909876895789904700198885093,
- 35711959237.35566804944018545154716670596,
- 17921034426.03720969991975575445893111267,
- 6039542586.35202800506429164430729792107,
- 1439720407.311721673663223072794912393972,
- 248874557.8620541565114603864132294232163,
- 31426415.58540019438061423162831820536287,
- 2876370.628935372441225409051620849613599,
- 186056.2653952234950402949897160456992822,
- 8071.672002365816210638002902272250613822,
- 210.8242777515793458725097339207133627117,
- 2.506628274631000270164908177133837338626
- };
- return evaluateRational(num, DENOM, z);
- }
- /**
- * Computes the Lanczos approximation scaled by {@code exp(g)}.
- *
- * @param z Argument.
- * @return the scaled Lanczos approximation.
- */
- static double lanczosSumExpGScaled(double z) {
- // As above with numerator divided by exp(g) = 413.509...
- final double[] num = {
- 56906521.91347156388090791033559122686859,
- 103794043.1163445451906271053616070238554,
- 86363131.28813859145546927288977868422342,
- 43338889.32467613834773723740590533316085,
- 14605578.08768506808414169982791359218571,
- 3481712.15498064590882071018964774556468,
- 601859.6171681098786670226533699352302507,
- 75999.29304014542649875303443598909137092,
- 6955.999602515376140356310115515198987526,
- 449.9445569063168119446858607650988409623,
- 19.51992788247617482847860966235652136208,
- 0.5098416655656676188125178644804694509993,
- 0.006061842346248906525783753964555936883222
- };
- return evaluateRational(num, DENOM, z);
- }
- /**
- * Evaluate the rational number as two polynomials.
- *
- * <p>Adapted from {@code boost/math/tools/detail/rational_horner3_13.hpp}.
- * Note: There are 3 variations of the unrolled rational evaluation.
- * These methods change the order based on the sign of x. This
- * should be used for the Lanczos code as this comment in
- * {@code boost/math/tools/rational.hpp} notes:
- *
- * <blockquote>
- * However, there
- * are some tricks we can use to prevent overflow that might otherwise
- * occur in polynomial evaluation, if z is large. This is important
- * in our Lanczos code for example.
- * </blockquote>
- *
- * @param a Coefficients of the numerator polynomial
- * @param b Coefficients of the denominator polynomial
- * @param x Value
- * @return the rational number
- */
- private static double evaluateRational(double[] a, int[] b, double x) {
- // The choice of algorithm in Boost is based on the compiler
- // to suite the available optimisations.
- //
- // Tests against rational_horner1_13.hpp which uses a first order
- // Horner method (no x*x term) show only minor variations in
- // error. rational_horner2_13.hpp has the same second order Horner
- // method with different code layout of the same sum.
- // rational_horner3_13.hpp
- if (x <= 1) {
- final double x2 = x * x;
- double t0 = a[12] * x2 + a[10];
- double t1 = a[11] * x2 + a[9];
- double t2 = b[12] * x2 + b[10];
- double t3 = b[11] * x2 + b[9];
- t0 *= x2;
- t1 *= x2;
- t2 *= x2;
- t3 *= x2;
- t0 += a[8];
- t1 += a[7];
- t2 += b[8];
- t3 += b[7];
- t0 *= x2;
- t1 *= x2;
- t2 *= x2;
- t3 *= x2;
- t0 += a[6];
- t1 += a[5];
- t2 += b[6];
- t3 += b[5];
- t0 *= x2;
- t1 *= x2;
- t2 *= x2;
- t3 *= x2;
- t0 += a[4];
- t1 += a[3];
- t2 += b[4];
- t3 += b[3];
- t0 *= x2;
- t1 *= x2;
- t2 *= x2;
- t3 *= x2;
- t0 += a[2];
- t1 += a[1];
- t2 += b[2];
- t3 += b[1];
- t0 *= x2;
- t2 *= x2;
- t0 += a[0];
- t2 += b[0];
- t1 *= x;
- t3 *= x;
- return (t0 + t1) / (t2 + t3);
- }
- final double z = 1 / x;
- final double z2 = 1 / (x * x);
- double t0 = a[0] * z2 + a[2];
- double t1 = a[1] * z2 + a[3];
- double t2 = b[0] * z2 + b[2];
- double t3 = b[1] * z2 + b[3];
- t0 *= z2;
- t1 *= z2;
- t2 *= z2;
- t3 *= z2;
- t0 += a[4];
- t1 += a[5];
- t2 += b[4];
- t3 += b[5];
- t0 *= z2;
- t1 *= z2;
- t2 *= z2;
- t3 *= z2;
- t0 += a[6];
- t1 += a[7];
- t2 += b[6];
- t3 += b[7];
- t0 *= z2;
- t1 *= z2;
- t2 *= z2;
- t3 *= z2;
- t0 += a[8];
- t1 += a[9];
- t2 += b[8];
- t3 += b[9];
- t0 *= z2;
- t1 *= z2;
- t2 *= z2;
- t3 *= z2;
- t0 += a[10];
- t1 += a[11];
- t2 += b[10];
- t3 += b[11];
- t0 *= z2;
- t2 *= z2;
- t0 += a[12];
- t2 += b[12];
- t1 *= z;
- t3 *= z;
- return (t0 + t1) / (t2 + t3);
- }
- // Not implemented:
- // lanczos_sum_near_1
- // lanczos_sum_near_2
- }
- /** Private constructor. */
- private BoostGamma() {
- // intentionally empty.
- }
- /**
- * All factorials that are representable as a double.
- * This data is exposed for testing.
- *
- * @return factorials
- */
- static double[] getFactorials() {
- return FACTORIAL.clone();
- }
- /**
- * Returns the factorial of n.
- * This is unchecked as an index out of bound exception will occur
- * if the value n is not within the range [0, 170].
- * This function is exposed for use in {@link BoostBeta}.
- *
- * @param n Argument n (must be in [0, 170])
- * @return n!
- */
- static double uncheckedFactorial(int n) {
- return FACTORIAL[n];
- }
- /**
- * Gamma function.
- *
- * <p>For small {@code z} this is based on the <em>NSWC Library of Mathematics
- * Subroutines</em> double precision implementation, {@code DGAMMA}.
- *
- * <p>For large {@code z} this is an implementation of the Boost C++ tgamma
- * function with Lanczos support.
- *
- * <p>Integers are handled using a look-up table of factorials.
- *
- * <p>Note: The Boost C++ implementation uses the Lanczos sum for all {@code z}.
- * When promotion of double to long double is not available this has larger
- * errors than the double precision specific NSWC implementation. For larger
- * {@code z} the Boost C++ Lanczos implementation incorporates the sqrt(2 pi)
- * factor and has lower error than the implementation using the
- * {@link LanczosApproximation} class.
- *
- * @param z Argument z
- * @return gamma value
- */
- static double tgamma(double z) {
- // Handle integers
- if (Math.rint(z) == z) {
- if (z <= 0) {
- // Pole error
- return Double.NaN;
- }
- if (z <= MAX_GAMMA_Z) {
- // Gamma(n) = (n-1)!
- return FACTORIAL[(int) z - 1];
- }
- // Overflow
- return Double.POSITIVE_INFINITY;
- }
- if (Math.abs(z) <= LANCZOS_THRESHOLD) {
- // Small z
- // NSWC Library of Mathematics Subroutines
- // Note:
- // This does not benefit from using extended precision to track the sum (t).
- // Extended precision on the product reduces the error but the majority
- // of error remains in InvGamma1pm1.
- if (z >= 1) {
- /*
- * From the recurrence relation
- * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
- * then
- * Gamma(t) = 1 / [1 + InvGamma1pm1.value(t - 1)],
- * where t = x - n. This means that t must satisfy
- * -0.5 <= t - 1 <= 1.5.
- */
- double prod = 1;
- double t = z;
- while (t > 2.5) {
- t -= 1;
- prod *= t;
- }
- return prod / (1 + InvGamma1pm1.value(t - 1));
- }
- /*
- * From the recurrence relation
- * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
- * then
- * Gamma(x + n + 1) = 1 / [1 + InvGamma1pm1.value(x + n)],
- * which requires -0.5 <= x + n <= 1.5.
- */
- double prod = z;
- double t = z;
- while (t < -0.5) {
- t += 1;
- prod *= t;
- }
- return 1 / (prod * (1 + InvGamma1pm1.value(t)));
- }
- // Large non-integer z
- // Boost C++ tgamma implementation
- if (z < 0) {
- /*
- * From the reflection formula
- * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
- * and the recurrence relation
- * Gamma(1 - x) = -x * Gamma(-x),
- * it is found
- * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
- */
- return -Math.PI / (sinpx(z) * tgamma(-z));
- } else if (z > MAX_GAMMA_Z + 1) {
- // Addition to the Boost code: Simple overflow detection
- return Double.POSITIVE_INFINITY;
- }
- double result = Lanczos.lanczosSum(z);
- final double zgh = z + Lanczos.GMH;
- final double lzgh = Math.log(zgh);
- if (z * lzgh > LOG_MAX_VALUE) {
- // we're going to overflow unless this is done with care:
- // Updated
- // Check for overflow removed:
- // if (lzgh * z / 2 > LOG_MAX_VALUE) ... overflow
- // This is replaced by checking z > MAX_FACTORIAL + 2
- final double hp = Math.pow(zgh, (z / 2) - 0.25);
- result *= hp / Math.exp(zgh);
- // Check for overflow has been removed:
- // if (Double.MAX_VALUE / hp < result) ... overflow
- result *= hp;
- } else {
- result *= Math.pow(zgh, z - 0.5) / Math.exp(zgh);
- }
- return result;
- }
- /**
- * Ad hoc function calculates x * sin(pi * x), taking extra care near when x is
- * near a whole number.
- *
- * @param x Value (assumed to be negative)
- * @return x * sin(pi * x)
- */
- static double sinpx(double x) {
- int sign = 1;
- // This is always called with a negative
- // if (x < 0)
- x = -x;
- double fl = Math.floor(x);
- double dist;
- if (isOdd(fl)) {
- fl += 1;
- dist = fl - x;
- sign = -sign;
- } else {
- dist = x - fl;
- }
- if (dist > 0.5f) {
- dist = 1 - dist;
- }
- final double result = Math.sin(dist * Math.PI);
- return sign * x * result;
- }
- /**
- * Checks if the value is odd.
- *
- * @param v Value (assumed to be positive and an integer)
- * @return true if odd
- */
- private static boolean isOdd(double v) {
- // Note:
- // Any value larger than 2^53 should be even.
- // If the input is positive then truncation of extreme doubles (>2^63)
- // to the primitive long creates an odd value: 2^63-1.
- // This is corrected by inverting the sign of v and the extreme is even: -2^63.
- // This function is never called when the argument is this large
- // as this is a pole error in tgamma so the effect is never observed.
- // However the isOdd function is correct for all positive finite v.
- return (((long) -v) & 0x1) == 1;
- }
- /**
- * Log Gamma function.
- * Defined as the natural logarithm of the absolute value of tgamma(z).
- *
- * @param z Argument z
- * @return log gamma value
- */
- static double lgamma(double z) {
- return lgamma(z, null);
- }
- /**
- * Log Gamma function.
- * Defined as the natural logarithm of the absolute value of tgamma(z).
- *
- * @param z Argument z
- * @param sign If a non-zero length array the first index is set on output to the sign of tgamma(z)
- * @return log gamma value
- */
- static double lgamma(double z, int[] sign) {
- double result;
- int sresult = 1;
- if (z <= -ROOT_EPSILON) {
- // reflection formula:
- if (Math.rint(z) == z) {
- // Pole error
- return Double.NaN;
- }
- double t = sinpx(z);
- z = -z;
- if (t < 0) {
- t = -t;
- } else {
- sresult = -sresult;
- }
- // This summation can have large magnitudes with opposite signs.
- // Use an extended precision sum to reduce cancellation.
- result = Sum.of(-lgamma(z)).add(-Math.log(t)).add(LOG_PI).getAsDouble();
- } else if (z < ROOT_EPSILON) {
- if (z == 0) {
- // Pole error
- return Double.NaN;
- }
- if (4 * Math.abs(z) < EPSILON) {
- result = -Math.log(Math.abs(z));
- } else {
- result = Math.log(Math.abs(1 / z - EULER));
- }
- if (z < 0) {
- sresult = -1;
- }
- } else if (z < 15) {
- result = lgammaSmall(z, z - 1, z - 2);
- // The z > 3 condition is always true
- //} else if (z > 3 && z < 100) {
- } else if (z < 100) {
- // taking the log of tgamma reduces the error, no danger of overflow here:
- result = Math.log(tgamma(z));
- } else {
- // regular evaluation:
- final double zgh = z + Lanczos.GMH;
- result = Math.log(zgh) - 1;
- result *= z - 0.5f;
- //
- // Only add on the lanczos sum part if we're going to need it:
- //
- if (result * EPSILON < 20) {
- result += Math.log(Lanczos.lanczosSumExpGScaled(z));
- }
- }
- if (nonZeroLength(sign)) {
- sign[0] = sresult;
- }
- return result;
- }
- /**
- * Log Gamma function for small z.
- *
- * @param z Argument z
- * @param zm1 {@code z - 1}
- * @param zm2 {@code z - 2}
- * @return log gamma value
- */
- private static double lgammaSmall(double z, double zm1, double zm2) {
- // This version uses rational approximations for small
- // values of z accurate enough for 64-bit mantissas
- // (80-bit long doubles), works well for 53-bit doubles as well.
- // Updated to use an extended precision sum
- final Sum result = Sum.create();
- // Note:
- // Removed z < EPSILON branch.
- // The function is called
- // from lgamma:
- // ROOT_EPSILON <= z < 15
- // from tgamma1pm1:
- // 1.5 <= z < 2
- // 1 <= z < 3
- if ((zm1 == 0) || (zm2 == 0)) {
- // nothing to do, result is zero....
- return 0;
- } else if (z > 2) {
- //
- // Begin by performing argument reduction until
- // z is in [2,3):
- //
- if (z >= 3) {
- do {
- z -= 1;
- result.add(Math.log(z));
- } while (z >= 3);
- // Update zm2, we need it below:
- zm2 = z - 2;
- }
- //
- // Use the following form:
- //
- // lgamma(z) = (z-2)(z+1)(Y + R(z-2))
- //
- // where R(z-2) is a rational approximation optimised for
- // low absolute error - as long as its absolute error
- // is small compared to the constant Y - then any rounding
- // error in its computation will get wiped out.
- //
- // R(z-2) has the following properties:
- //
- // At double: Max error found: 4.231e-18
- // At long double: Max error found: 1.987e-21
- // Maximum Deviation Found (approximation error): 5.900e-24
- //
- double P;
- P = -0.324588649825948492091e-4;
- P = -0.541009869215204396339e-3 + P * zm2;
- P = -0.259453563205438108893e-3 + P * zm2;
- P = 0.172491608709613993966e-1 + P * zm2;
- P = 0.494103151567532234274e-1 + P * zm2;
- P = 0.25126649619989678683e-1 + P * zm2;
- P = -0.180355685678449379109e-1 + P * zm2;
- double Q;
- Q = -0.223352763208617092964e-6;
- Q = 0.224936291922115757597e-3 + Q * zm2;
- Q = 0.82130967464889339326e-2 + Q * zm2;
- Q = 0.988504251128010129477e-1 + Q * zm2;
- Q = 0.541391432071720958364e0 + Q * zm2;
- Q = 0.148019669424231326694e1 + Q * zm2;
- Q = 0.196202987197795200688e1 + Q * zm2;
- Q = 0.1e1 + Q * zm2;
- final float Y = 0.158963680267333984375e0f;
- final double r = zm2 * (z + 1);
- final double R = P / Q;
- result.addProduct(r, Y).addProduct(r, R);
- } else {
- //
- // If z is less than 1 use recurrence to shift to
- // z in the interval [1,2]:
- //
- if (z < 1) {
- result.add(-Math.log(z));
- zm2 = zm1;
- zm1 = z;
- z += 1;
- }
- //
- // Two approximations, one for z in [1,1.5] and
- // one for z in [1.5,2]:
- //
- if (z <= 1.5) {
- //
- // Use the following form:
- //
- // lgamma(z) = (z-1)(z-2)(Y + R(z-1
- //
- // where R(z-1) is a rational approximation optimised for
- // low absolute error - as long as its absolute error
- // is small compared to the constant Y - then any rounding
- // error in its computation will get wiped out.
- //
- // R(z-1) has the following properties:
- //
- // At double precision: Max error found: 1.230011e-17
- // At 80-bit long double precision: Max error found: 5.631355e-21
- // Maximum Deviation Found: 3.139e-021
- // Expected Error Term: 3.139e-021
- //
- final float Y = 0.52815341949462890625f;
- double P;
- P = -0.100346687696279557415e-2;
- P = -0.240149820648571559892e-1 + P * zm1;
- P = -0.158413586390692192217e0 + P * zm1;
- P = -0.406567124211938417342e0 + P * zm1;
- P = -0.414983358359495381969e0 + P * zm1;
- P = -0.969117530159521214579e-1 + P * zm1;
- P = 0.490622454069039543534e-1 + P * zm1;
- double Q;
- Q = 0.195768102601107189171e-2;
- Q = 0.577039722690451849648e-1 + Q * zm1;
- Q = 0.507137738614363510846e0 + Q * zm1;
- Q = 0.191415588274426679201e1 + Q * zm1;
- Q = 0.348739585360723852576e1 + Q * zm1;
- Q = 0.302349829846463038743e1 + Q * zm1;
- Q = 0.1e1 + Q * zm1;
- final double r = P / Q;
- final double prefix = zm1 * zm2;
- result.addProduct(prefix, Y).addProduct(prefix, r);
- } else {
- //
- // Use the following form:
- //
- // lgamma(z) = (2-z)(1-z)(Y + R(2-z
- //
- // where R(2-z) is a rational approximation optimised for
- // low absolute error - as long as its absolute error
- // is small compared to the constant Y - then any rounding
- // error in its computation will get wiped out.
- //
- // R(2-z) has the following properties:
- //
- // At double precision, max error found: 1.797565e-17
- // At 80-bit long double precision, max error found: 9.306419e-21
- // Maximum Deviation Found: 2.151e-021
- // Expected Error Term: 2.150e-021
- //
- final float Y = 0.452017307281494140625f;
- final double mzm2 = -zm2;
- double P;
- P = 0.431171342679297331241e-3;
- P = -0.850535976868336437746e-2 + P * mzm2;
- P = 0.542809694055053558157e-1 + P * mzm2;
- P = -0.142440390738631274135e0 + P * mzm2;
- P = 0.144216267757192309184e0 + P * mzm2;
- P = -0.292329721830270012337e-1 + P * mzm2;
- double Q;
- Q = -0.827193521891290553639e-6;
- Q = -0.100666795539143372762e-2 + Q * mzm2;
- Q = 0.25582797155975869989e-1 + Q * mzm2;
- Q = -0.220095151814995745555e0 + Q * mzm2;
- Q = 0.846973248876495016101e0 + Q * mzm2;
- Q = -0.150169356054485044494e1 + Q * mzm2;
- Q = 0.1e1 + Q * mzm2;
- final double r = zm2 * zm1;
- final double R = P / Q;
- result.addProduct(r, Y).addProduct(r, R);
- }
- }
- return result.getAsDouble();
- }
- /**
- * Calculates tgamma(1+dz)-1.
- *
- * @param dz Argument
- * @return tgamma(1+dz)-1
- */
- static double tgamma1pm1(double dz) {
- //
- // This helper calculates tgamma(1+dz)-1 without cancellation errors,
- // used by the upper incomplete gamma with z < 1:
- //
- final double result;
- if (dz < 0) {
- if (dz < -0.5) {
- // Best method is simply to subtract 1 from tgamma:
- result = tgamma(1 + dz) - 1;
- } else {
- // Use expm1 on lgamma:
- result = Math.expm1(-Math.log1p(dz) + lgammaSmall(dz + 2, dz + 1, dz));
- }
- } else {
- if (dz < 2) {
- // Use expm1 on lgamma:
- result = Math.expm1(lgammaSmall(dz + 1, dz, dz - 1));
- } else {
- // Best method is simply to subtract 1 from tgamma:
- result = tgamma(1 + dz) - 1;
- }
- }
- return result;
- }
- /**
- * Full upper incomplete gamma.
- *
- * @param a Argument a
- * @param x Argument x
- * @return upper gamma value
- */
- static double tgamma(double a, double x) {
- return gammaIncompleteImp(a, x, false, true, Policy.getDefault());
- }
- /**
- * Full upper incomplete gamma.
- *
- * @param a Argument a
- * @param x Argument x
- * @param policy Function evaluation policy
- * @return upper gamma value
- */
- static double tgamma(double a, double x, Policy policy) {
- return gammaIncompleteImp(a, x, false, true, policy);
- }
- /**
- * Full lower incomplete gamma.
- *
- * @param a Argument a
- * @param x Argument x
- * @return lower gamma value
- */
- static double tgammaLower(double a, double x) {
- return gammaIncompleteImp(a, x, false, false, Policy.getDefault());
- }
- /**
- * Full lower incomplete gamma.
- *
- * @param a Argument a
- * @param x Argument x
- * @param policy Function evaluation policy
- * @return lower gamma value
- */
- static double tgammaLower(double a, double x, Policy policy) {
- return gammaIncompleteImp(a, x, false, false, policy);
- }
- /**
- * Regularised upper incomplete gamma.
- *
- * @param a Argument a
- * @param x Argument x
- * @return q
- */
- static double gammaQ(double a, double x) {
- return gammaIncompleteImp(a, x, true, true, Policy.getDefault());
- }
- /**
- * Regularised upper incomplete gamma.
- *
- * @param a Argument a
- * @param x Argument x
- * @param policy Function evaluation policy
- * @return q
- */
- static double gammaQ(double a, double x, Policy policy) {
- return gammaIncompleteImp(a, x, true, true, policy);
- }
- /**
- * Regularised lower incomplete gamma.
- *
- * @param a Argument a
- * @param x Argument x
- * @return p
- */
- static double gammaP(double a, double x) {
- return gammaIncompleteImp(a, x, true, false, Policy.getDefault());
- }
- /**
- * Regularised lower incomplete gamma.
- *
- * @param a Argument a
- * @param x Argument x
- * @param policy Function evaluation policy
- * @return p
- */
- static double gammaP(double a, double x, Policy policy) {
- return gammaIncompleteImp(a, x, true, false, policy);
- }
- /**
- * Derivative of the regularised lower incomplete gamma.
- * <p>\( \frac{e^{-x} x^{a-1}}{\Gamma(a)} \)
- *
- * <p>Adapted from {@code boost::math::detail::gamma_p_derivative_imp}
- *
- * @param a Argument a
- * @param x Argument x
- * @return p derivative
- */
- static double gammaPDerivative(double a, double x) {
- //
- // Usual error checks first:
- //
- if (Double.isNaN(a) || Double.isNaN(x) || a <= 0 || x < 0) {
- return Double.NaN;
- }
- //
- // Now special cases:
- //
- if (x == 0) {
- if (a > 1) {
- return 0;
- }
- return (a == 1) ? 1 : Double.POSITIVE_INFINITY;
- }
- //
- // Normal case:
- //
- double f1 = regularisedGammaPrefix(a, x);
- if (f1 == 0) {
- // Underflow in calculation, use logs instead:
- f1 = a * Math.log(x) - x - lgamma(a) - Math.log(x);
- f1 = Math.exp(f1);
- } else {
- // Will overflow when (x < 1) && (Double.MAX_VALUE * x < f1).
- // There is no exception for this case so just return the result.
- f1 /= x;
- }
- return f1;
- }
- /**
- * Main incomplete gamma entry point, handles all four incomplete gammas.
- * Adapted from {@code boost::math::detail::gamma_incomplete_imp}.
- *
- * <p>The Boost code has a pointer {@code p_derivative} that can be set to the
- * value of the derivative. This is used for the inverse incomplete
- * gamma functions {@code gamma_(p|q)_inv_imp}. It is not required for the forward
- * evaluation functions.
- *
- * @param a Argument a
- * @param x Argument x
- * @param normalised true to compute the regularised value
- * @param invert true to compute the upper value Q (default is lower value P)
- * @param pol Function evaluation policy
- * @return gamma value
- */
- private static double gammaIncompleteImp(double a, double x,
- boolean normalised, boolean invert, Policy pol) {
- if (Double.isNaN(a) || Double.isNaN(x) || a <= 0 || x < 0) {
- return Double.NaN;
- }
- double result = 0;
- if (a >= MAX_FACTORIAL && !normalised) {
- //
- // When we're computing the non-normalized incomplete gamma
- // and a is large the result is rather hard to compute unless
- // we use logs. There are really two options - if x is a long
- // way from a in value then we can reliably use methods 2 and 4
- // below in logarithmic form and go straight to the result.
- // Otherwise we let the regularized gamma take the strain
- // (the result is unlikely to underflow in the central region anyway)
- // and combine with lgamma in the hopes that we get a finite result.
- //
- if (invert && (a * 4 < x)) {
- // This is method 4 below, done in logs:
- result = a * Math.log(x) - x;
- result += Math.log(upperGammaFraction(a, x, pol));
- } else if (!invert && (a > 4 * x)) {
- // This is method 2 below, done in logs:
- result = a * Math.log(x) - x;
- result += Math.log(lowerGammaSeries(a, x, 0, pol) / a);
- } else {
- result = gammaIncompleteImp(a, x, true, invert, pol);
- if (result == 0) {
- if (invert) {
- // Try http://functions.wolfram.com/06.06.06.0039.01
- result = 1 + 1 / (12 * a) + 1 / (288 * a * a);
- result = Math.log(result) - a + (a - 0.5f) * Math.log(a) + LOG_ROOT_TWO_PI;
- } else {
- // This is method 2 below, done in logs, we're really outside the
- // range of this method, but since the result is almost certainly
- // infinite, we should probably be OK:
- result = a * Math.log(x) - x;
- result += Math.log(lowerGammaSeries(a, x, 0, pol) / a);
- }
- } else {
- result = Math.log(result) + lgamma(a);
- }
- }
- // If result is > log(MAX_VALUE) the result will overflow.
- // There is no exception for this case so just return the result.
- return Math.exp(result);
- }
- final boolean isInt;
- final boolean isHalfInt;
- // Update. x must be safe for exp(-x). Change to -x > LOG_MIN_VALUE.
- final boolean isSmallA = (a < 30) && (a <= x + 1) && (-x > LOG_MIN_VALUE);
- if (isSmallA) {
- final double fa = Math.floor(a);
- isInt = fa == a;
- isHalfInt = !isInt && (Math.abs(fa - a) == 0.5f);
- } else {
- isInt = isHalfInt = false;
- }
- final int evalMethod;
- if (isInt && (x > 0.6)) {
- // calculate Q via finite sum:
- invert = !invert;
- evalMethod = 0;
- } else if (isHalfInt && (x > 0.2)) {
- // calculate Q via finite sum for half integer a:
- invert = !invert;
- evalMethod = 1;
- } else if ((x < ROOT_EPSILON) && (a > 1)) {
- evalMethod = 6;
- } else if ((x > 1000) && (a < x * 0.75f)) {
- // Note:
- // The branch is used in Boost when:
- // ((x > 1000) && ((a < x) || (Math.abs(a - 50) / x < 1)))
- //
- // This case was added after Boost 1_68_0.
- // See: https://github.com/boostorg/math/issues/168
- //
- // When using only double precision for the evaluation
- // it is a source of error when a ~ z as the asymptotic approximation
- // sums terms t_n+1 = t_n * (a - n - 1) / z starting from t_0 = 1.
- // These terms are close to 1 when a ~ z and the sum has many terms
- // with reduced precision.
- // This has been updated to allow only cases with fast convergence.
- // It will be used when x -> infinity and a << x.
- // calculate Q via asymptotic approximation:
- invert = !invert;
- evalMethod = 7;
- } else if (x < 0.5) {
- //
- // Changeover criterion chosen to give a changeover at Q ~ 0.33
- //
- if (-0.4 / Math.log(x) < a) {
- // Compute P
- evalMethod = 2;
- } else {
- evalMethod = 3;
- }
- } else if (x < 1.1) {
- //
- // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
- //
- if (x * 0.75f < a) {
- // Compute P
- evalMethod = 2;
- } else {
- evalMethod = 3;
- }
- } else {
- //
- // Begin by testing whether we're in the "bad" zone
- // where the result will be near 0.5 and the usual
- // series and continued fractions are slow to converge:
- //
- boolean useTemme = false;
- if (normalised && (a > 20)) {
- final double sigma = Math.abs((x - a) / a);
- if (a > 200) {
- //
- // This limit is chosen so that we use Temme's expansion
- // only if the result would be larger than about 10^-6.
- // Below that the regular series and continued fractions
- // converge OK, and if we use Temme's method we get increasing
- // errors from the dominant erfc term as its (inexact) argument
- // increases in magnitude.
- //
- if (20 / a > sigma * sigma) {
- useTemme = true;
- }
- } else {
- // Note in this zone we can't use Temme's expansion for
- // types longer than an 80-bit real:
- // it would require too many terms in the polynomials.
- if (sigma < 0.4) {
- useTemme = true;
- }
- }
- }
- if (useTemme) {
- evalMethod = 5;
- } else {
- //
- // Regular case where the result will not be too close to 0.5.
- //
- // Changeover here occurs at P ~ Q ~ 0.5
- // Note that series computation of P is about x2 faster than continued fraction
- // calculation of Q, so try and use the CF only when really necessary,
- // especially for small x.
- //
- if (x - (1 / (3 * x)) < a) {
- evalMethod = 2;
- } else {
- evalMethod = 4;
- invert = !invert;
- }
- }
- }
- switch (evalMethod) {
- case 0:
- result = finiteGammaQ(a, x);
- if (!normalised) {
- result *= tgamma(a);
- }
- break;
- case 1:
- result = finiteHalfGammaQ(a, x);
- if (!normalised) {
- result *= tgamma(a);
- }
- break;
- case 2:
- // Compute P:
- result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
- if (result != 0) {
- //
- // If we're going to be inverting the result then we can
- // reduce the number of series evaluations by quite
- // a few iterations if we set an initial value for the
- // series sum based on what we'll end up subtracting it from
- // at the end.
- // Have to be careful though that this optimization doesn't
- // lead to spurious numeric overflow. Note that the
- // scary/expensive overflow checks below are more often
- // than not bypassed in practice for "sensible" input
- // values:
- //
- double initValue = 0;
- boolean optimisedInvert = false;
- if (invert) {
- initValue = normalised ? 1 : tgamma(a);
- if (normalised || (result >= 1) || (Double.MAX_VALUE * result > initValue)) {
- initValue /= result;
- if (normalised || (a < 1) || (Double.MAX_VALUE / a > initValue)) {
- initValue *= -a;
- optimisedInvert = true;
- } else {
- initValue = 0;
- }
- } else {
- initValue = 0;
- }
- }
- result *= lowerGammaSeries(a, x, initValue, pol) / a;
- if (optimisedInvert) {
- invert = false;
- result = -result;
- }
- }
- break;
- case 3:
- // Compute Q:
- invert = !invert;
- final double[] g = {0};
- result = tgammaSmallUpperPart(a, x, pol, g, invert);
- invert = false;
- if (normalised) {
- // Addition to the Boost code:
- if (g[0] == Double.POSITIVE_INFINITY) {
- // Very small a will overflow gamma(a). Resort to logs.
- // This method requires improvement as the error is very large.
- // It is better than returning zero for a non-zero result.
- result = Math.exp(Math.log(result) - lgamma(a));
- } else {
- result /= g[0];
- }
- }
- break;
- case 4:
- // Compute Q:
- result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
- if (result != 0) {
- result *= upperGammaFraction(a, x, pol);
- }
- break;
- case 5:
- // Call 53-bit version
- result = igammaTemmeLarge(a, x);
- if (x >= a) {
- invert = !invert;
- }
- break;
- case 6:
- // x is so small that P is necessarily very small too, use
- // http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
- if (normalised) {
- // If tgamma overflows then result = 0
- result = Math.pow(x, a) / tgamma(a + 1);
- } else {
- result = Math.pow(x, a) / a;
- }
- result *= 1 - a * x / (a + 1);
- break;
- case 7:
- default:
- // x is large,
- // Compute Q:
- result = normalised ? regularisedGammaPrefix(a, x) : fullIgammaPrefix(a, x);
- result /= x;
- if (result != 0) {
- result *= incompleteTgammaLargeX(a, x, pol);
- }
- break;
- }
- if (normalised && (result > 1)) {
- result = 1;
- }
- if (invert) {
- final double gam = normalised ? 1 : tgamma(a);
- result = gam - result;
- }
- return result;
- }
- /**
- * Upper gamma fraction.
- * Multiply result by z^a * e^-z to get the full
- * upper incomplete integral. Divide by tgamma(z)
- * to normalise.
- *
- * @param a Argument a
- * @param z Argument z
- * @param pol Function evaluation policy
- * @return upper gamma fraction
- */
- // This is package-private for testing
- static double upperGammaFraction(double a, double z, Policy pol) {
- final double eps = pol.getEps();
- final int maxIterations = pol.getMaxIterations();
- // This is computing:
- // 1
- // ------------------------------
- // b0 + a1 / (b1 + a2 )
- // -------------
- // b2 + a3
- // --------
- // b3 + ...
- //
- // b0 = z - a + 1
- // a1 = a - 1
- //
- // It can be done several ways with variations in accuracy.
- // The current implementation has the best accuracy and matches the Boost code.
- final double zma1 = z - a + 1;
- final Supplier<Coefficient> gen = new Supplier<Coefficient>() {
- /** Iteration. */
- private int k;
- @Override
- public Coefficient get() {
- ++k;
- return Coefficient.of(k * (a - k), zma1 + 2.0 * k);
- }
- };
- return 1 / GeneralizedContinuedFraction.value(zma1, gen, eps, maxIterations);
- }
- /**
- * Upper gamma fraction for integer a.
- * Called when {@code a < 30} and {@code -x > LOG_MIN_VALUE}.
- *
- * @param a Argument a (assumed to be small)
- * @param x Argument x
- * @return upper gamma fraction
- */
- private static double finiteGammaQ(double a, double x) {
- //
- // Calculates normalised Q when a is an integer:
- //
- // Update:
- // Assume -x > log min value and no underflow to zero.
- double sum = Math.exp(-x);
- double term = sum;
- for (int n = 1; n < a; ++n) {
- term /= n;
- term *= x;
- sum += term;
- }
- return sum;
- }
- /**
- * Upper gamma fraction for half integer a.
- * Called when {@code a < 30} and {@code -x > LOG_MIN_VALUE}.
- *
- * @param a Argument a (assumed to be small)
- * @param x Argument x
- * @return upper gamma fraction
- */
- private static double finiteHalfGammaQ(double a, double x) {
- //
- // Calculates normalised Q when a is a half-integer:
- //
- // Update:
- // Assume -x > log min value:
- // erfc(sqrt(708)) = erfc(26.6) => erfc has a non-zero value
- double e = BoostErf.erfc(Math.sqrt(x));
- if (a > 1) {
- double term = Math.exp(-x) / Math.sqrt(Math.PI * x);
- term *= x;
- term /= 0.5;
- double sum = term;
- for (int n = 2; n < a; ++n) {
- term /= n - 0.5;
- term *= x;
- sum += term;
- }
- e += sum;
- }
- return e;
- }
- /**
- * Lower gamma series.
- * Multiply result by ((z^a) * (e^-z) / a) to get the full
- * lower incomplete integral. Then divide by tgamma(a)
- * to get the normalised value.
- *
- * @param a Argument a
- * @param z Argument z
- * @param initValue Initial value
- * @param pol Function evaluation policy
- * @return lower gamma series
- */
- // This is package-private for testing
- static double lowerGammaSeries(double a, double z, double initValue, Policy pol) {
- final double eps = pol.getEps();
- final int maxIterations = pol.getMaxIterations();
- // Lower gamma series representation.
- final DoubleSupplier gen = new DoubleSupplier() {
- /** Next result. */
- private double result = 1;
- /** Iteration. */
- private int n;
- @Override
- public double getAsDouble() {
- final double r = result;
- n++;
- result *= z / (a + n);
- return r;
- }
- };
- return BoostTools.kahanSumSeries(gen, eps, maxIterations, initValue);
- }
- /**
- * Upper gamma fraction for very small a.
- *
- * @param a Argument a (assumed to be small)
- * @param x Argument x
- * @param pol Function evaluation policy
- * @param pgam set to value of gamma(a) on output
- * @param invert true to invert the result
- * @return upper gamma fraction
- */
- private static double tgammaSmallUpperPart(double a, double x, Policy pol, double[] pgam, boolean invert) {
- //
- // Compute the full upper fraction (Q) when a is very small:
- //
- double result;
- result = tgamma1pm1(a);
- // Note: Replacing this with tgamma(a) does not reduce error on current test data.
- // gamma(1+z) = z gamma(z)
- // pgam[0] == gamma(a)
- pgam[0] = (result + 1) / a;
- double p = BoostMath.powm1(x, a);
- result -= p;
- result /= a;
- // Removed subtraction of 10 from this value
- final int maxIter = pol.getMaxIterations();
- p += 1;
- final double initValue = invert ? pgam[0] : 0;
- // Series representation for upper fraction when z is small.
- final DoubleSupplier gen = new DoubleSupplier() {
- /** Result term. */
- private double result = -x;
- /** Argument x (this is negated on purpose). */
- private final double z = -x;
- /** Iteration. */
- private int n = 1;
- @Override
- public double getAsDouble() {
- final double r = result / (a + n);
- n++;
- result = result * z / n;
- return r;
- }
- };
- result = -p * BoostTools.kahanSumSeries(gen, pol.getEps(), maxIter, (initValue - result) / p);
- if (invert) {
- result = -result;
- }
- return result;
- }
- /**
- * Calculate power term prefix (z^a)(e^-z) used in the non-normalised
- * incomplete gammas.
- *
- * @param a Argument a
- * @param z Argument z
- * @return incomplete gamma prefix
- */
- static double fullIgammaPrefix(double a, double z) {
- if (z > Double.MAX_VALUE) {
- return 0;
- }
- final double alz = a * Math.log(z);
- final double prefix;
- if (z >= 1) {
- if ((alz < LOG_MAX_VALUE) && (-z > LOG_MIN_VALUE)) {
- prefix = Math.pow(z, a) * Math.exp(-z);
- } else if (a >= 1) {
- prefix = Math.pow(z / Math.exp(z / a), a);
- } else {
- prefix = Math.exp(alz - z);
- }
- } else {
- if (alz > LOG_MIN_VALUE) {
- prefix = Math.pow(z, a) * Math.exp(-z);
- } else {
- // Updated to remove unreachable final branch using Math.exp(alz - z).
- // This branch requires (z / a < LOG_MAX_VALUE) to avoid overflow in exp.
- // At this point:
- // 1. log(z) is negative;
- // 2. a * log(z) <= -708 requires a > -708 / log(z).
- // For any z < 1: -708 / log(z) is > z. Thus a is > z and
- // z / a < LOG_MAX_VALUE is always true.
- prefix = Math.pow(z / Math.exp(z / a), a);
- }
- }
- // Removed overflow check. Return infinity if it occurs.
- return prefix;
- }
- /**
- * Compute (z^a)(e^-z)/tgamma(a).
- * <p>Most of the error occurs in this function.
- *
- * @param a Argument a
- * @param z Argument z
- * @return regularized gamma prefix
- */
- // This is package-private for testing
- static double regularisedGammaPrefix(double a, double z) {
- if (z >= Double.MAX_VALUE) {
- return 0;
- }
- // Update this condition from: a < 1
- if (a <= 1) {
- //
- // We have to treat a < 1 as a special case because our Lanczos
- // approximations are optimised against the factorials with a > 1,
- // and for high precision types especially (128-bit reals for example)
- // very small values of a can give rather erroneous results for gamma
- // unless we do this:
- //
- // Boost todo: is this still required? Lanczos approx should be better now?
- //
- // Update this condition from: z <= LOG_MIN_VALUE
- // Require exp(-z) to not underflow:
- // -z > log(min_value)
- if (-z <= LOG_MIN_VALUE) {
- // Oh dear, have to use logs, should be free of cancellation errors though:
- return Math.exp(a * Math.log(z) - z - lgamma(a));
- }
- // direct calculation, no danger of overflow as gamma(a) < 1/a
- // for small a.
- return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
- }
- // Update to the Boost code.
- // Use some of the logic from fullIgammaPrefix(a, z) to use the direct
- // computation if it is valid. Assuming pow and exp are accurate to 1 ULP it
- // puts most of the error in evaluation of tgamma(a). This is accurate
- // enough that this reduces max error on the current test data.
- //
- // Overflow cases fall-through to the Lanczos approximation that incorporates
- // the pow and exp terms used in the tgamma(a) computation with the terms
- // z^a and e^-z into a single evaluation of pow and exp. See equation 15:
- // https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html
- if (a <= MAX_GAMMA_Z) {
- final double alz1 = a * Math.log(z);
- if (z >= 1) {
- if ((alz1 < LOG_MAX_VALUE) && (-z > LOG_MIN_VALUE)) {
- return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
- }
- } else if (alz1 > LOG_MIN_VALUE) {
- return Math.pow(z, a) * Math.exp(-z) / tgamma(a);
- }
- }
- //
- // For smallish a and x combining the power terms with the Lanczos approximation
- // gives the greatest accuracy
- //
- final double agh = a + Lanczos.GMH;
- double prefix;
- final double factor = Math.sqrt(agh / Math.E) / Lanczos.lanczosSumExpGScaled(a);
- // Update to the Boost code.
- // Lower threshold for large a from 150 to 128 and compute d on demand.
- // See NUMBERS-179.
- if (a > 128) {
- final double d = ((z - a) - Lanczos.GMH) / agh;
- if (Math.abs(d * d * a) <= 100) {
- // special case for large a and a ~ z.
- // When a and x are large, we end up with a very large exponent with a base near one:
- // this will not be computed accurately via the pow function, and taking logs simply
- // leads to cancellation errors.
- prefix = a * SpecialMath.log1pmx(d) + z * -Lanczos.GMH / agh;
- prefix = Math.exp(prefix);
- return prefix * factor;
- }
- }
- //
- // general case.
- // direct computation is most accurate, but use various fallbacks
- // for different parts of the problem domain:
- //
- final double alz = a * Math.log(z / agh);
- final double amz = a - z;
- if ((Math.min(alz, amz) <= LOG_MIN_VALUE) || (Math.max(alz, amz) >= LOG_MAX_VALUE)) {
- final double amza = amz / a;
- if ((Math.min(alz, amz) / 2 > LOG_MIN_VALUE) && (Math.max(alz, amz) / 2 < LOG_MAX_VALUE)) {
- // compute square root of the result and then square it:
- final double sq = Math.pow(z / agh, a / 2) * Math.exp(amz / 2);
- prefix = sq * sq;
- } else if ((Math.min(alz, amz) / 4 > LOG_MIN_VALUE) &&
- (Math.max(alz, amz) / 4 < LOG_MAX_VALUE) && (z > a)) {
- // compute the 4th root of the result then square it twice:
- final double sq = Math.pow(z / agh, a / 4) * Math.exp(amz / 4);
- prefix = sq * sq;
- prefix *= prefix;
- } else if ((amza > LOG_MIN_VALUE) && (amza < LOG_MAX_VALUE)) {
- prefix = Math.pow((z * Math.exp(amza)) / agh, a);
- } else {
- prefix = Math.exp(alz + amz);
- }
- } else {
- prefix = Math.pow(z / agh, a) * Math.exp(amz);
- }
- prefix *= factor;
- return prefix;
- }
- /**
- * Implements the asymptotic expansions of the incomplete
- * gamma functions P(a, x) and Q(a, x), used when a is large and
- * x ~ a.
- *
- * <p>The primary reference is:
- * <pre>
- * "The Asymptotic Expansion of the Incomplete Gamma Functions"
- * N. M. Temme.
- * Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
- * </pre>
- *
- * <p>A different way of evaluating these expansions,
- * plus a lot of very useful background information is in:
- * <pre>
- * "A Set of Algorithms For the Incomplete Gamma Functions."
- * N. M. Temme.
- * Probability in the Engineering and Informational Sciences,
- * 8, 1994, 291.
- * </pre>
- *
- * <p>An alternative implementation is in:
- * <pre>
- * "Computation of the Incomplete Gamma Function Ratios and their Inverse."
- * A. R. Didonato and A. H. Morris.
- * ACM TOMS, Vol 12, No 4, Dec 1986, p377.
- * </pre>
- *
- * <p>This is a port of the function accurate for 53-bit mantissas
- * (IEEE double precision or 10^-17). To understand the code, refer to Didonato
- * and Morris, from Eq 17 and 18 onwards.
- *
- * <p>The coefficients used here are not taken from Didonato and Morris:
- * the domain over which these expansions are used is slightly different
- * to theirs, and their constants are not quite accurate enough for
- * 128-bit long doubles. Instead the coefficients were calculated
- * using the methods described by Temme p762 from Eq 3.8 onwards.
- * The values obtained agree with those obtained by Didonato and Morris
- * (at least to the first 30 digits that they provide).
- * At double precision the degrees of polynomial required for full
- * machine precision are close to those recommended to Didonato and Morris,
- * but of course many more terms are needed for larger types.
- *
- * <p>Adapted from {@code boost/math/special_functions/detail/igamma_large.hpp}.
- *
- * @param a the a
- * @param x the x
- * @return the double
- */
- // This is package-private for testing
- static double igammaTemmeLarge(double a, double x) {
- final double sigma = (x - a) / a;
- final double phi = -SpecialMath.log1pmx(sigma);
- final double y = a * phi;
- double z = Math.sqrt(2 * phi);
- if (x < a) {
- z = -z;
- }
- // The following polynomials are evaluated with a loop
- // with Horner's method. Variations exist using
- // a second order Horner's method with an unrolled loop.
- // These are chosen in Boost based on the C++ compiler.
- // For example:
- // boost/math/tools/detail/polynomial_horner1_15.hpp
- // boost/math/tools/detail/polynomial_horner2_15.hpp
- // boost/math/tools/detail/polynomial_horner3_15.hpp
- final double[] workspace = new double[10];
- final double[] C0 = {
- -0.33333333333333333,
- 0.083333333333333333,
- -0.014814814814814815,
- 0.0011574074074074074,
- 0.0003527336860670194,
- -0.00017875514403292181,
- 0.39192631785224378e-4,
- -0.21854485106799922e-5,
- -0.185406221071516e-5,
- 0.8296711340953086e-6,
- -0.17665952736826079e-6,
- 0.67078535434014986e-8,
- 0.10261809784240308e-7,
- -0.43820360184533532e-8,
- 0.91476995822367902e-9,
- };
- workspace[0] = BoostTools.evaluatePolynomial(C0, z);
- final double[] C1 = {
- -0.0018518518518518519,
- -0.0034722222222222222,
- 0.0026455026455026455,
- -0.00099022633744855967,
- 0.00020576131687242798,
- -0.40187757201646091e-6,
- -0.18098550334489978e-4,
- 0.76491609160811101e-5,
- -0.16120900894563446e-5,
- 0.46471278028074343e-8,
- 0.1378633446915721e-6,
- -0.5752545603517705e-7,
- 0.11951628599778147e-7,
- };
- workspace[1] = BoostTools.evaluatePolynomial(C1, z);
- final double[] C2 = {
- 0.0041335978835978836,
- -0.0026813271604938272,
- 0.00077160493827160494,
- 0.20093878600823045e-5,
- -0.00010736653226365161,
- 0.52923448829120125e-4,
- -0.12760635188618728e-4,
- 0.34235787340961381e-7,
- 0.13721957309062933e-5,
- -0.6298992138380055e-6,
- 0.14280614206064242e-6,
- };
- workspace[2] = BoostTools.evaluatePolynomial(C2, z);
- final double[] C3 = {
- 0.00064943415637860082,
- 0.00022947209362139918,
- -0.00046918949439525571,
- 0.00026772063206283885,
- -0.75618016718839764e-4,
- -0.23965051138672967e-6,
- 0.11082654115347302e-4,
- -0.56749528269915966e-5,
- 0.14230900732435884e-5,
- };
- workspace[3] = BoostTools.evaluatePolynomial(C3, z);
- final double[] C4 = {
- -0.0008618882909167117,
- 0.00078403922172006663,
- -0.00029907248030319018,
- -0.14638452578843418e-5,
- 0.66414982154651222e-4,
- -0.39683650471794347e-4,
- 0.11375726970678419e-4,
- };
- workspace[4] = BoostTools.evaluatePolynomial(C4, z);
- final double[] C5 = {
- -0.00033679855336635815,
- -0.69728137583658578e-4,
- 0.00027727532449593921,
- -0.00019932570516188848,
- 0.67977804779372078e-4,
- 0.1419062920643967e-6,
- -0.13594048189768693e-4,
- 0.80184702563342015e-5,
- -0.22914811765080952e-5,
- };
- workspace[5] = BoostTools.evaluatePolynomial(C5, z);
- final double[] C6 = {
- 0.00053130793646399222,
- -0.00059216643735369388,
- 0.00027087820967180448,
- 0.79023532326603279e-6,
- -0.81539693675619688e-4,
- 0.56116827531062497e-4,
- -0.18329116582843376e-4,
- };
- workspace[6] = BoostTools.evaluatePolynomial(C6, z);
- final double[] C7 = {
- 0.00034436760689237767,
- 0.51717909082605922e-4,
- -0.00033493161081142236,
- 0.0002812695154763237,
- -0.00010976582244684731,
- };
- workspace[7] = BoostTools.evaluatePolynomial(C7, z);
- final double[] C8 = {
- -0.00065262391859530942,
- 0.00083949872067208728,
- -0.00043829709854172101,
- };
- workspace[8] = BoostTools.evaluatePolynomial(C8, z);
- workspace[9] = -0.00059676129019274625;
- double result = BoostTools.evaluatePolynomial(workspace, 1 / a);
- result *= Math.exp(-y) / Math.sqrt(2 * Math.PI * a);
- if (x < a) {
- result = -result;
- }
- result += BoostErf.erfc(Math.sqrt(y)) / 2;
- return result;
- }
- /**
- * Incomplete tgamma for large X.
- *
- * <p>This summation is a source of error as the series starts at 1 and descends to zero.
- * It can have thousands of iterations when a and z are large and close in value.
- *
- * @param a Argument a
- * @param x Argument x
- * @param pol Function evaluation policy
- * @return incomplete tgamma
- */
- // This is package-private for testing
- static double incompleteTgammaLargeX(double a, double x, Policy pol) {
- final double eps = pol.getEps();
- final int maxIterations = pol.getMaxIterations();
- // Asymptotic approximation for large argument, see: https://dlmf.nist.gov/8.11#E2.
- final DoubleSupplier gen = new DoubleSupplier() {
- /** Result term. */
- private double term = 1;
- /** Iteration. */
- private int n;
- @Override
- public double getAsDouble() {
- final double result = term;
- n++;
- term *= (a - n) / x;
- return result;
- }
- };
- return BoostTools.kahanSumSeries(gen, eps, maxIterations);
- }
- /**
- * Return true if the array is not null and has non-zero length.
- *
- * @param array Array
- * @return true if a non-zero length array
- */
- private static boolean nonZeroLength(int[] array) {
- return array != null && array.length != 0;
- }
- /**
- * Ratio of gamma functions.
- *
- * <p>\[ tgamma_ratio(z, delta) = \frac{\Gamma(z)}{\Gamma(z + delta)} \]
- *
- * <p>Adapted from {@code tgamma_delta_ratio_imp}. The use of
- * {@code max_factorial<double>::value == 170} has been replaced with
- * {@code MAX_GAMMA_Z == 171}. This threshold is used when it is possible
- * to call the gamma function without overflow.
- *
- * @param z Argument z
- * @param delta The difference
- * @return gamma ratio
- */
- static double tgammaDeltaRatio(double z, double delta) {
- final double zDelta = z + delta;
- if (Double.isNaN(zDelta)) {
- // One or both arguments are NaN
- return Double.NaN;
- }
- if (z <= 0 || zDelta <= 0) {
- // This isn't very sophisticated, or accurate, but it does work:
- return tgamma(z) / tgamma(zDelta);
- }
- // Note: Directly calling tgamma(z) / tgamma(z + delta) if possible
- // without overflow is not more accurate
- if (Math.rint(delta) == delta) {
- if (delta == 0) {
- return 1;
- }
- //
- // If both z and delta are integers, see if we can just use table lookup
- // of the factorials to get the result:
- //
- if (Math.rint(z) == z &&
- z <= MAX_GAMMA_Z && zDelta <= MAX_GAMMA_Z) {
- return FACTORIAL[(int) z - 1] / FACTORIAL[(int) zDelta - 1];
- }
- if (Math.abs(delta) < 20) {
- //
- // delta is a small integer, we can use a finite product:
- //
- if (delta < 0) {
- z -= 1;
- double result = z;
- for (int d = (int) (delta + 1); d != 0; d++) {
- z -= 1;
- result *= z;
- }
- return result;
- }
- double result = 1 / z;
- for (int d = (int) (delta - 1); d != 0; d--) {
- z += 1;
- result /= z;
- }
- return result;
- }
- }
- return tgammaDeltaRatioImpLanczos(z, delta);
- }
- /**
- * Ratio of gamma functions using Lanczos support.
- *
- * <p>\[ tgamma_delta_ratio(z, delta) = \frac{\Gamma(z)}{\Gamma(z + delta)}
- *
- * <p>Adapted from {@code tgamma_delta_ratio_imp_lanczos}. The use of
- * {@code max_factorial<double>::value == 170} has been replaced with
- * {@code MAX_GAMMA_Z == 171}. This threshold is used when it is possible
- * to use the precomputed factorial table.
- *
- * @param z Argument z
- * @param delta The difference
- * @return gamma ratio
- */
- private static double tgammaDeltaRatioImpLanczos(double z, double delta) {
- if (z < EPSILON) {
- //
- // We get spurious numeric overflow unless we're very careful, this
- // can occur either inside Lanczos::lanczos_sum(z) or in the
- // final combination of terms, to avoid this, split the product up
- // into 2 (or 3) parts:
- //
- // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta
- // z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial
- //
- if (MAX_GAMMA_Z < delta) {
- double ratio = tgammaDeltaRatioImpLanczos(delta, MAX_GAMMA_Z - delta);
- ratio *= z;
- ratio *= FACTORIAL[MAX_FACTORIAL];
- return 1 / ratio;
- }
- return 1 / (z * tgamma(z + delta));
- }
- final double zgh = z + Lanczos.GMH;
- double result;
- if (z + delta == z) {
- // Here:
- // lanczosSum(z) / lanczosSum(z + delta) == 1
- // Update to the Boost code to remove unreachable code:
- // Given z + delta == z then |delta / z| < EPSILON; and z < zgh
- // assume (abs(delta / zgh) < EPSILON) and remove unreachable
- // else branch which sets result = 1
- // We have:
- // result = exp((0.5 - z) * log1p(delta / zgh));
- // 0.5 - z == -z
- // log1p(delta / zgh) = delta / zgh = delta / z
- // multiplying we get -delta.
- // Note:
- // This can be different from
- // exp((0.5 - z) * log1p(delta / zgh)) when z is small.
- // In this case the result is exp(small) and the result
- // is within 1 ULP of 1.0. This is left as the original
- // Boost method using exp(-delta).
- result = Math.exp(-delta);
- } else {
- if (Math.abs(delta) < 10) {
- result = Math.exp((0.5 - z) * Math.log1p(delta / zgh));
- } else {
- result = Math.pow(zgh / (zgh + delta), z - 0.5);
- }
- // Split the calculation up to avoid spurious overflow:
- result *= Lanczos.lanczosSum(z) / Lanczos.lanczosSum(z + delta);
- }
- result *= Math.pow(Math.E / (zgh + delta), delta);
- return result;
- }
- /**
- * Ratio of gamma functions.
- *
- * <p>\[ tgamma_ratio(x, y) = \frac{\Gamma(x)}{\Gamma(y)}
- *
- * <p>Adapted from {@code tgamma_ratio_imp}. The use of
- * {@code max_factorial<double>::value == 170} has been replaced with
- * {@code MAX_GAMMA_Z == 171}. This threshold is used when it is possible
- * to call the gamma function without overflow.
- *
- * @param x Argument x (must be positive finite)
- * @param y Argument y (must be positive finite)
- * @return gamma ratio (or nan)
- */
- static double tgammaRatio(double x, double y) {
- if (x <= 0 || !Double.isFinite(x) || y <= 0 || !Double.isFinite(y)) {
- return Double.NaN;
- }
- if (x <= Double.MIN_NORMAL) {
- // Special case for denorms...Ugh.
- return TWO_POW_53 * tgammaRatio(x * TWO_POW_53, y);
- }
- if (x <= MAX_GAMMA_Z && y <= MAX_GAMMA_Z) {
- // Rather than subtracting values, lets just call the gamma functions directly:
- return tgamma(x) / tgamma(y);
- }
- double prefix = 1;
- if (x < 1) {
- if (y < 2 * MAX_GAMMA_Z) {
- // We need to sidestep on x as well, otherwise we'll underflow
- // before we get to factor in the prefix term:
- prefix /= x;
- x += 1;
- while (y >= MAX_GAMMA_Z) {
- y -= 1;
- prefix /= y;
- }
- return prefix * tgamma(x) / tgamma(y);
- }
- //
- // result is almost certainly going to underflow to zero, try logs just in case:
- //
- return Math.exp(lgamma(x) - lgamma(y));
- }
- if (y < 1) {
- if (x < 2 * MAX_GAMMA_Z) {
- // We need to sidestep on y as well, otherwise we'll overflow
- // before we get to factor in the prefix term:
- prefix *= y;
- y += 1;
- while (x >= MAX_GAMMA_Z) {
- x -= 1;
- prefix *= x;
- }
- return prefix * tgamma(x) / tgamma(y);
- }
- //
- // Result will almost certainly overflow, try logs just in case:
- //
- return Math.exp(lgamma(x) - lgamma(y));
- }
- //
- // Regular case, x and y both large and similar in magnitude:
- //
- return tgammaDeltaRatio(x, y - x);
- }
- }