Gamma.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.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  17. /*
  18.  * This is not the original file distributed by the Apache Software Foundation
  19.  * It has been modified by the Hipparchus project
  20.  */
  21. package org.hipparchus.special;

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.Field;
  24. import org.hipparchus.exception.LocalizedCoreFormats;
  25. import org.hipparchus.exception.MathIllegalArgumentException;
  26. import org.hipparchus.exception.MathIllegalStateException;
  27. import org.hipparchus.util.ContinuedFraction;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.FieldContinuedFraction;

  30. /**
  31.  * <p>
  32.  * This is a utility class that provides computation methods related to the
  33.  * &Gamma; (Gamma) family of functions.
  34.  * </p>
  35.  * <p>
  36.  * Implementation of {@link #invGamma1pm1(double)} and
  37.  * {@link #logGamma1p(double)} is based on the algorithms described in
  38.  * </p>
  39.  * <ul>
  40.  * <li><a href="http://dx.doi.org/10.1145/22721.23109">Didonato and Morris
  41.  * (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios and
  42.  *     their Inverse</em>, TOMS 12(4), 377-393,</li>
  43.  * <li><a href="http://dx.doi.org/10.1145/131766.131776">Didonato and Morris
  44.  * (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the
  45.  *     Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li>
  46.  * </ul>
  47.  * <p>
  48.  * and implemented in the
  49.  * <a href="http://www.dtic.mil/docs/citations/ADA476840">NSWC Library of Mathematical Functions</a>,
  50.  * available
  51.  * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html">here</a>.
  52.  * This library is "approved for public release", and the
  53.  * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf">Copyright guidance</a>
  54.  * indicates that unless otherwise stated in the code, all FORTRAN functions in
  55.  * this library are license free. Since no such notice appears in the code these
  56.  * functions can safely be ported to Hipparchus.
  57.  * </p>
  58.  *
  59.  */
  60. public class Gamma {
  61.     /**
  62.      * <a href="http://en.wikipedia.org/wiki/Euler-Mascheroni_constant">Euler-Mascheroni constant</a>
  63.      */
  64.     public static final double GAMMA = 0.577215664901532860606512090082; // NOPMD - the fact the function and the constant have the same name is intentional and comes from mathematics conventions

  65.     /**
  66.      * The value of the {@code g} constant in the Lanczos approximation, see
  67.      * {@link #lanczos(double)}.
  68.      */
  69.     public static final double LANCZOS_G = 607.0 / 128.0;

  70.     /** Maximum allowed numerical error. */
  71.     private static final double DEFAULT_EPSILON = 10e-15;

  72.     /** Lanczos coefficients */
  73.     private static final double[] LANCZOS = {
  74.         0.99999999999999709182,
  75.         57.156235665862923517,
  76.         -59.597960355475491248,
  77.         14.136097974741747174,
  78.         -0.49191381609762019978,
  79.         .33994649984811888699e-4,
  80.         .46523628927048575665e-4,
  81.         -.98374475304879564677e-4,
  82.         .15808870322491248884e-3,
  83.         -.21026444172410488319e-3,
  84.         .21743961811521264320e-3,
  85.         -.16431810653676389022e-3,
  86.         .84418223983852743293e-4,
  87.         -.26190838401581408670e-4,
  88.         .36899182659531622704e-5,
  89.     };

  90.     /** Avoid repeated computation of log of 2 PI in logGamma */
  91.     private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(2.0 * FastMath.PI);

  92.     /** The constant value of &radic;(2&pi;). */
  93.     private static final double SQRT_TWO_PI = 2.506628274631000502;

  94.     // limits for switching algorithm in digamma
  95.     /** C limit. */
  96.     private static final double C_LIMIT = 49;

  97.     /** S limit. */
  98.     private static final double S_LIMIT = 1e-8;

  99.     /*
  100.      * Constants for the computation of double invGamma1pm1(double).
  101.      * Copied from DGAM1 in the NSWC library.
  102.      */

  103.     /** The constant {@code A0} defined in {@code DGAM1}. */
  104.     private static final double INV_GAMMA1P_M1_A0 = .611609510448141581788E-08;

  105.     /** The constant {@code A1} defined in {@code DGAM1}. */
  106.     private static final double INV_GAMMA1P_M1_A1 = .624730830116465516210E-08;

  107.     /** The constant {@code B1} defined in {@code DGAM1}. */
  108.     private static final double INV_GAMMA1P_M1_B1 = .203610414066806987300E+00;

  109.     /** The constant {@code B2} defined in {@code DGAM1}. */
  110.     private static final double INV_GAMMA1P_M1_B2 = .266205348428949217746E-01;

  111.     /** The constant {@code B3} defined in {@code DGAM1}. */
  112.     private static final double INV_GAMMA1P_M1_B3 = .493944979382446875238E-03;

  113.     /** The constant {@code B4} defined in {@code DGAM1}. */
  114.     private static final double INV_GAMMA1P_M1_B4 = -.851419432440314906588E-05;

  115.     /** The constant {@code B5} defined in {@code DGAM1}. */
  116.     private static final double INV_GAMMA1P_M1_B5 = -.643045481779353022248E-05;

  117.     /** The constant {@code B6} defined in {@code DGAM1}. */
  118.     private static final double INV_GAMMA1P_M1_B6 = .992641840672773722196E-06;

  119.     /** The constant {@code B7} defined in {@code DGAM1}. */
  120.     private static final double INV_GAMMA1P_M1_B7 = -.607761895722825260739E-07;

  121.     /** The constant {@code B8} defined in {@code DGAM1}. */
  122.     private static final double INV_GAMMA1P_M1_B8 = .195755836614639731882E-09;

  123.     /** The constant {@code P0} defined in {@code DGAM1}. */
  124.     private static final double INV_GAMMA1P_M1_P0 = .6116095104481415817861E-08;

  125.     /** The constant {@code P1} defined in {@code DGAM1}. */
  126.     private static final double INV_GAMMA1P_M1_P1 = .6871674113067198736152E-08;

  127.     /** The constant {@code P2} defined in {@code DGAM1}. */
  128.     private static final double INV_GAMMA1P_M1_P2 = .6820161668496170657918E-09;

  129.     /** The constant {@code P3} defined in {@code DGAM1}. */
  130.     private static final double INV_GAMMA1P_M1_P3 = .4686843322948848031080E-10;

  131.     /** The constant {@code P4} defined in {@code DGAM1}. */
  132.     private static final double INV_GAMMA1P_M1_P4 = .1572833027710446286995E-11;

  133.     /** The constant {@code P5} defined in {@code DGAM1}. */
  134.     private static final double INV_GAMMA1P_M1_P5 = -.1249441572276366213222E-12;

  135.     /** The constant {@code P6} defined in {@code DGAM1}. */
  136.     private static final double INV_GAMMA1P_M1_P6 = .4343529937408594255178E-14;

  137.     /** The constant {@code Q1} defined in {@code DGAM1}. */
  138.     private static final double INV_GAMMA1P_M1_Q1 = .3056961078365221025009E+00;

  139.     /** The constant {@code Q2} defined in {@code DGAM1}. */
  140.     private static final double INV_GAMMA1P_M1_Q2 = .5464213086042296536016E-01;

  141.     /** The constant {@code Q3} defined in {@code DGAM1}. */
  142.     private static final double INV_GAMMA1P_M1_Q3 = .4956830093825887312020E-02;

  143.     /** The constant {@code Q4} defined in {@code DGAM1}. */
  144.     private static final double INV_GAMMA1P_M1_Q4 = .2692369466186361192876E-03;

  145.     /** The constant {@code C} defined in {@code DGAM1}. */
  146.     private static final double INV_GAMMA1P_M1_C = -.422784335098467139393487909917598E+00;

  147.     /** The constant {@code C0} defined in {@code DGAM1}. */
  148.     private static final double INV_GAMMA1P_M1_C0 = .577215664901532860606512090082402E+00;

  149.     /** The constant {@code C1} defined in {@code DGAM1}. */
  150.     private static final double INV_GAMMA1P_M1_C1 = -.655878071520253881077019515145390E+00;

  151.     /** The constant {@code C2} defined in {@code DGAM1}. */
  152.     private static final double INV_GAMMA1P_M1_C2 = -.420026350340952355290039348754298E-01;

  153.     /** The constant {@code C3} defined in {@code DGAM1}. */
  154.     private static final double INV_GAMMA1P_M1_C3 = .166538611382291489501700795102105E+00;

  155.     /** The constant {@code C4} defined in {@code DGAM1}. */
  156.     private static final double INV_GAMMA1P_M1_C4 = -.421977345555443367482083012891874E-01;

  157.     /** The constant {@code C5} defined in {@code DGAM1}. */
  158.     private static final double INV_GAMMA1P_M1_C5 = -.962197152787697356211492167234820E-02;

  159.     /** The constant {@code C6} defined in {@code DGAM1}. */
  160.     private static final double INV_GAMMA1P_M1_C6 = .721894324666309954239501034044657E-02;

  161.     /** The constant {@code C7} defined in {@code DGAM1}. */
  162.     private static final double INV_GAMMA1P_M1_C7 = -.116516759185906511211397108401839E-02;

  163.     /** The constant {@code C8} defined in {@code DGAM1}. */
  164.     private static final double INV_GAMMA1P_M1_C8 = -.215241674114950972815729963053648E-03;

  165.     /** The constant {@code C9} defined in {@code DGAM1}. */
  166.     private static final double INV_GAMMA1P_M1_C9 = .128050282388116186153198626328164E-03;

  167.     /** The constant {@code C10} defined in {@code DGAM1}. */
  168.     private static final double INV_GAMMA1P_M1_C10 = -.201348547807882386556893914210218E-04;

  169.     /** The constant {@code C11} defined in {@code DGAM1}. */
  170.     private static final double INV_GAMMA1P_M1_C11 = -.125049348214267065734535947383309E-05;

  171.     /** The constant {@code C12} defined in {@code DGAM1}. */
  172.     private static final double INV_GAMMA1P_M1_C12 = .113302723198169588237412962033074E-05;

  173.     /** The constant {@code C13} defined in {@code DGAM1}. */
  174.     private static final double INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06;

  175.     /**
  176.      * Default constructor.  Prohibit instantiation.
  177.      */
  178.     private Gamma() {}

  179.     /**
  180.      * <p>
  181.      * Returns the value of log&nbsp;&Gamma;(x) for x&nbsp;&gt;&nbsp;0.
  182.      * </p>
  183.      * <p>
  184.      * For x &le; 8, the implementation is based on the double precision
  185.      * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
  186.      * {@code DGAMLN}. For x &gt; 8, the implementation is based on
  187.      * </p>
  188.      * <ul>
  189.      * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
  190.      *     Function</a>, equation (28).</li>
  191.      * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
  192.      *     Lanczos Approximation</a>, equations (1) through (5).</li>
  193.      * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
  194.      *     the computation of the convergent Lanczos complex Gamma
  195.      *     approximation</a></li>
  196.      * </ul>
  197.      *
  198.      * @param x Argument.
  199.      * @return the value of {@code log(Gamma(x))}, {@code Double.NaN} if
  200.      * {@code x <= 0.0}.
  201.      */
  202.     public static double logGamma(double x) {
  203.         double ret;

  204.         if (Double.isNaN(x) || (x <= 0.0)) {
  205.             ret = Double.NaN;
  206.         } else if (x < 0.5) {
  207.             return logGamma1p(x) - FastMath.log(x);
  208.         } else if (x <= 2.5) {
  209.             return logGamma1p((x - 0.5) - 0.5);
  210.         } else if (x <= 8.0) {
  211.             final int n = (int) FastMath.floor(x - 1.5);
  212.             double prod = 1.0;
  213.             for (int i = 1; i <= n; i++) {
  214.                 prod *= x - i;
  215.             }
  216.             return logGamma1p(x - (n + 1)) + FastMath.log(prod);
  217.         } else {
  218.             double sum = lanczos(x);
  219.             double tmp = x + LANCZOS_G + .5;
  220.             ret = ((x + .5) * FastMath.log(tmp)) - tmp +
  221.                 HALF_LOG_2_PI + FastMath.log(sum / x);
  222.         }

  223.         return ret;
  224.     }

  225.     /**
  226.      * <p>
  227.      * Returns the value of log&nbsp;&Gamma;(x) for x&nbsp;&gt;&nbsp;0.
  228.      * </p>
  229.      * <p>
  230.      * For x &le; 8, the implementation is based on the double precision
  231.      * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
  232.      * {@code DGAMLN}. For x &gt; 8, the implementation is based on
  233.      * </p>
  234.      * <ul>
  235.      * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
  236.      *     Function</a>, equation (28).</li>
  237.      * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
  238.      *     Lanczos Approximation</a>, equations (1) through (5).</li>
  239.      * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
  240.      *     the computation of the convergent Lanczos complex Gamma
  241.      *     approximation</a></li>
  242.      * </ul>
  243.      *
  244.      * @param x Argument.
  245.      * @param <T> Type of the field elements.
  246.      * @return the value of {@code log(Gamma(x))}, {@code Double.NaN} if
  247.      * {@code x <= 0.0}.
  248.      */
  249.     public static <T extends CalculusFieldElement<T>> T logGamma(T x) {
  250.         final Field<T> field = x.getField();
  251.         T              ret;

  252.         if (x.isNaN() || (x.getReal() <= 0.0)) {
  253.             ret = field.getOne().multiply(Double.NaN);
  254.         }
  255.         else if (x.getReal() < 0.5) {
  256.             return logGamma1p(x).subtract(x.log());
  257.         }
  258.         else if (x.getReal() <= 2.5) {
  259.             return logGamma1p(x.subtract(1));
  260.         }
  261.         else if (x.getReal() <= 8.0) {
  262.             final int n    = (int) x.subtract(1.5).floor().getReal();
  263.             T         prod = field.getOne();
  264.             for (int i = 1; i <= n; i++) {
  265.                 prod = prod.multiply(x.subtract(i));
  266.             }
  267.             return logGamma1p(x.subtract(n + 1)).add(prod.log());
  268.         }
  269.         else {
  270.             T sum = lanczos(x);
  271.             T tmp = x.add(LANCZOS_G + .5);
  272.             ret = x.add(.5).multiply(tmp.log()).subtract(tmp).add(HALF_LOG_2_PI).add(sum.divide(x).log());
  273.         }

  274.         return ret;
  275.     }

  276.     /**
  277.      * Returns the regularized gamma function P(a, x).
  278.      *
  279.      * @param a Parameter.
  280.      * @param x Value.
  281.      * @return the regularized gamma function P(a, x).
  282.      * @throws MathIllegalStateException if the algorithm fails to converge.
  283.      */
  284.     public static double regularizedGammaP(double a, double x) {
  285.         return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
  286.     }

  287.     /**
  288.      * Returns the regularized gamma function P(a, x).
  289.      *
  290.      * @param a Parameter.
  291.      * @param x Value.
  292.      * @param <T> Type of the field elements.
  293.      * @return the regularized gamma function P(a, x).
  294.      * @throws MathIllegalStateException if the algorithm fails to converge.
  295.      */
  296.     public static <T extends CalculusFieldElement<T>> T regularizedGammaP(T a, T x) {
  297.         return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
  298.     }

  299.     /**
  300.      * Returns the regularized gamma function P(a, x).
  301.      * <p>
  302.      * The implementation of this method is based on:
  303.      * <ul>
  304.      *  <li>
  305.      *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
  306.      *   Regularized Gamma Function</a>, equation (1)
  307.      *  </li>
  308.      *  <li>
  309.      *   <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
  310.      *   Incomplete Gamma Function</a>, equation (4).
  311.      *  </li>
  312.      *  <li>
  313.      *   <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html">
  314.      *   Confluent Hypergeometric Function of the First Kind</a>, equation (1).
  315.      *  </li>
  316.      * </ul>
  317.      *
  318.      * @param a the a parameter.
  319.      * @param x the value.
  320.      * @param epsilon When the absolute value of the nth item in the
  321.      * series is less than epsilon the approximation ceases to calculate
  322.      * further elements in the series.
  323.      * @param maxIterations Maximum number of "iterations" to complete.
  324.      * @return the regularized gamma function P(a, x)
  325.      * @throws MathIllegalStateException if the algorithm fails to converge.
  326.      */
  327.     public static double regularizedGammaP(double a,
  328.                                            double x,
  329.                                            double epsilon,
  330.                                            int maxIterations) {
  331.         double ret;

  332.         if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
  333.             ret = Double.NaN;
  334.         } else if (x == 0.0) {
  335.             ret = 0.0;
  336.         } else if (x >= a + 1) {
  337.             // use regularizedGammaQ because it should converge faster in this
  338.             // case.
  339.             ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
  340.         } else {
  341.             // calculate series
  342.             double n = 0.0; // current element index
  343.             double an = 1.0 / a; // n-th element in the series
  344.             double sum = an; // partial sum
  345.             while (FastMath.abs(an/sum) > epsilon &&
  346.                    n < maxIterations &&
  347.                    sum < Double.POSITIVE_INFINITY) {
  348.                 // compute next element in the series
  349.                 n += 1.0;
  350.                 an *= x / (a + n);

  351.                 // update partial sum
  352.                 sum += an;
  353.             }
  354.             if (n >= maxIterations) {
  355.                 throw new MathIllegalStateException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, maxIterations);
  356.             } else if (Double.isInfinite(sum)) {
  357.                 ret = 1.0;
  358.             } else {
  359.                 ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * sum;
  360.             }
  361.         }

  362.         return ret;
  363.     }

  364.     /**
  365.      * Returns the regularized gamma function P(a, x).
  366.      * <p>
  367.      * The implementation of this method is based on:
  368.      * <ul>
  369.      *  <li>
  370.      *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
  371.      *   Regularized Gamma Function</a>, equation (1)
  372.      *  </li>
  373.      *  <li>
  374.      *   <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
  375.      *   Incomplete Gamma Function</a>, equation (4).
  376.      *  </li>
  377.      *  <li>
  378.      *   <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html">
  379.      *   Confluent Hypergeometric Function of the First Kind</a>, equation (1).
  380.      *  </li>
  381.      * </ul>
  382.      *
  383.      * @param a the a parameter.
  384.      * @param x the value.
  385.      * @param epsilon When the absolute value of the nth item in the
  386.      * series is less than epsilon the approximation ceases to calculate
  387.      * further elements in the series.
  388.      * @param maxIterations Maximum number of "iterations" to complete.
  389.      * @param <T> Type of the field elements.
  390.      * @return the regularized gamma function P(a, x)
  391.      * @throws MathIllegalStateException if the algorithm fails to converge.
  392.      */
  393.     public static <T extends CalculusFieldElement<T>> T regularizedGammaP(T a,
  394.                                            T x,
  395.                                            double epsilon,
  396.                                            int maxIterations) {
  397.         final Field<T> field = x.getField();
  398.         final T        zero  = field.getZero();
  399.         final T        one   = field.getOne();

  400.         T ret;

  401.         if (a.isNaN() || x.isNaN() || (a.getReal() <= 0.0) || (x.getReal() < 0.0)) {
  402.             ret = one.multiply(Double.NaN);
  403.         }
  404.         else if (x.getReal() == 0.0) {
  405.             ret = zero;
  406.         }
  407.         else if (x.getReal() >= a.add(1).getReal()) {
  408.             // use regularizedGammaQ because it should converge faster in this
  409.             // case.
  410.             ret = one.subtract(regularizedGammaQ(a, x, epsilon, maxIterations));
  411.         }
  412.         else {
  413.             // calculate series
  414.             double n   = 0.0; // current element index
  415.             T      an  = one.divide(a); // n-th element in the series
  416.             T      sum = an; // partial sum
  417.             while (an.divide(sum).abs().getReal() > epsilon &&
  418.                     n < maxIterations &&
  419.                     sum.getReal() < Double.POSITIVE_INFINITY) {
  420.                 // compute next element in the series
  421.                 n += 1.0;
  422.                 an = an.multiply(x.divide(a.add(n)));

  423.                 // update partial sum
  424.                 sum = sum.add(an);
  425.             }
  426.             if (n >= maxIterations) {
  427.                 throw new MathIllegalStateException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, maxIterations);
  428.             }
  429.             else if (sum.isInfinite()) {
  430.                 ret = one;
  431.             }
  432.             else {
  433.                 ret = a.multiply(x.log()).subtract(logGamma(a)).subtract(x).exp().multiply(sum);
  434.             }
  435.         }

  436.         return ret;
  437.     }

  438.     /**
  439.      * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
  440.      *
  441.      * @param a the a parameter.
  442.      * @param x the value.
  443.      * @return the regularized gamma function Q(a, x)
  444.      * @throws MathIllegalStateException if the algorithm fails to converge.
  445.      */
  446.     public static double regularizedGammaQ(double a, double x) {
  447.         return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
  448.     }

  449.     /**
  450.      * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
  451.      *
  452.      * @param a the a parameter.
  453.      * @param x the value.
  454.      * @param <T> Type of the field elements.
  455.      * @return the regularized gamma function Q(a, x)
  456.      * @throws MathIllegalStateException if the algorithm fails to converge.
  457.      */
  458.     public static <T extends CalculusFieldElement<T>> T regularizedGammaQ(T a, T x) {
  459.         return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
  460.     }

  461.     /**
  462.      * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
  463.      * <p>
  464.      * The implementation of this method is based on:
  465.      * <ul>
  466.      *  <li>
  467.      *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
  468.      *   Regularized Gamma Function</a>, equation (1).
  469.      *  </li>
  470.      *  <li>
  471.      *   <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
  472.      *   Regularized incomplete gamma function: Continued fraction representations
  473.      *   (formula 06.08.10.0003)</a>
  474.      *  </li>
  475.      * </ul>
  476.      *
  477.      * @param a the a parameter.
  478.      * @param x the value.
  479.      * @param epsilon When the absolute value of the nth item in the
  480.      * series is less than epsilon the approximation ceases to calculate
  481.      * further elements in the series.
  482.      * @param maxIterations Maximum number of "iterations" to complete.
  483.      * @return the regularized gamma function P(a, x)
  484.      * @throws MathIllegalStateException if the algorithm fails to converge.
  485.      */
  486.     public static double regularizedGammaQ(final double a,
  487.                                            double x,
  488.                                            double epsilon,
  489.                                            int maxIterations) {
  490.         double ret;

  491.         if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
  492.             ret = Double.NaN;
  493.         } else if (x == 0.0) {
  494.             ret = 1.0;
  495.         } else if (x < a + 1.0) {
  496.             // use regularizedGammaP because it should converge faster in this
  497.             // case.
  498.             ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations);
  499.         } else {
  500.             // create continued fraction
  501.             ContinuedFraction cf = new ContinuedFraction() {

  502.                 /** {@inheritDoc} */
  503.                 @Override
  504.                 protected double getA(int n, double x) {
  505.                     return ((2.0 * n) + 1.0) - a + x;
  506.                 }

  507.                 /** {@inheritDoc} */
  508.                 @Override
  509.                 protected double getB(int n, double x) {
  510.                     return n * (a - n);
  511.                 }
  512.             };

  513.             ret = 1.0 / cf.evaluate(x, epsilon, maxIterations);
  514.             ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * ret;
  515.         }

  516.         return ret;
  517.     }

  518.     /**
  519.      * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
  520.      * <p>
  521.      * The implementation of this method is based on:
  522.      * <ul>
  523.      *  <li>
  524.      *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
  525.      *   Regularized Gamma Function</a>, equation (1).
  526.      *  </li>
  527.      *  <li>
  528.      *   <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
  529.      *   Regularized incomplete gamma function: Continued fraction representations
  530.      *   (formula 06.08.10.0003)</a>
  531.      *  </li>
  532.      * </ul>
  533.      *
  534.      * @param a the a parameter.
  535.      * @param x the value.
  536.      * @param epsilon When the absolute value of the nth item in the
  537.      * series is less than epsilon the approximation ceases to calculate
  538.      * further elements in the series.
  539.      * @param maxIterations Maximum number of "iterations" to complete.
  540.      * @param <T> Type fo the field elements.
  541.      * @return the regularized gamma function P(a, x)
  542.      * @throws MathIllegalStateException if the algorithm fails to converge.
  543.      */
  544.     public static <T extends CalculusFieldElement<T>> T regularizedGammaQ(final T a,
  545.                                                                           T x,
  546.                                                                           double epsilon,
  547.                                                                           int maxIterations) {
  548.         final Field<T> field = x.getField();
  549.         final T        one   = field.getOne();

  550.         T ret;

  551.         if (a.isNaN() || x.isNaN() || a.getReal() <= 0.0 || x.getReal() < 0.0) {
  552.             ret = field.getOne().multiply(Double.NaN);
  553.         }
  554.         else if (x.getReal() == 0.0) {
  555.             ret = one;
  556.         }
  557.         else if (x.getReal() < a.add(1.0).getReal()) {
  558.             // use regularizedGammaP because it should converge faster in this
  559.             // case.
  560.             ret = one.subtract(regularizedGammaP(a, x, epsilon, maxIterations));
  561.         }
  562.         else {
  563.             // create continued fraction
  564.             FieldContinuedFraction cf = new FieldContinuedFraction() {

  565.                 /** {@inheritDoc} */
  566.                 @Override
  567.                 @SuppressWarnings("unchecked")
  568.                 public <C extends CalculusFieldElement<C>> C getA(final int n, final C x) {
  569.                     return x.subtract((C) a).add((2.0 * n) + 1.0);
  570.                 }

  571.                 /** {@inheritDoc} */
  572.                 @Override
  573.                 @SuppressWarnings("unchecked")
  574.                 public <C extends CalculusFieldElement<C>> C getB(final int n, final C x) {
  575.                     return (C) a.subtract(n).multiply(n);
  576.                 }
  577.             };

  578.             ret = one.divide(cf.evaluate(x, epsilon, maxIterations));
  579.             ret = a.multiply(x.log()).subtract(logGamma(a)).subtract(x).exp().multiply(ret);
  580.         }

  581.         return ret;
  582.     }

  583.     /**
  584.      * <p>Computes the digamma function of x.</p>
  585.      *
  586.      * <p>This is an independently written implementation of the algorithm described in
  587.      * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p>
  588.      *
  589.      * <p>Some of the constants have been changed to increase accuracy at the moderate expense
  590.      * of run-time.  The result should be accurate to within 10^-8 absolute tolerance for
  591.      * x &gt;= 10^-5 and within 10^-8 relative tolerance for x &gt; 0.</p>
  592.      *
  593.      * <p>Performance for large negative values of x will be quite expensive (proportional to
  594.      * |x|).  Accuracy for negative values of x should be about 10^-8 absolute for results
  595.      * less than 10^5 and 10^-8 relative for results larger than that.</p>
  596.      *
  597.      * @param x Argument.
  598.      * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller.
  599.      * @see <a href="http://en.wikipedia.org/wiki/Digamma_function">Digamma</a>
  600.      * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf">Bernardo&#39;s original article </a>
  601.      */
  602.     public static double digamma(double x) {
  603.         if (Double.isNaN(x) || Double.isInfinite(x)) {
  604.             return x;
  605.         }

  606.         if (x > 0 && x <= S_LIMIT) {
  607.             // use method 5 from Bernardo AS103
  608.             // accurate to O(x)
  609.             return -GAMMA - 1 / x;
  610.         }
  611.         if (x >= C_LIMIT) {
  612.             // use method 8 (accurate to O(1/x^8))
  613.             double inv = 1 / (x * x);
  614.             //            1       1        1         1         1         5           691         1
  615.             // log(x) -  --- - ------ + ------- - ------- + ------- - ------- +  ---------- - -------
  616.             //           2 x   12 x^2   120 x^4   252 x^6   240 x^8   660 x^10   32760 x^12   12 x^14
  617.             return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv * (1.0 / 252 + inv *
  618.                     (1.0 / 240 - inv * (5.0 / 660 + inv * (691.0 / 32760 - inv / 12))))));
  619.         }

  620.         return digamma(x + 1) - 1 / x;
  621.     }

  622.     /**
  623.      * <p>Computes the digamma function of x.</p>
  624.      *
  625.      * <p>This is an independently written implementation of the algorithm described in
  626.      * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p>
  627.      *
  628.      * <p>Some of the constants have been changed to increase accuracy at the moderate expense
  629.      * of run-time.  The result should be accurate to within 10^-8 absolute tolerance for
  630.      * x &gt;= 10^-5 and within 10^-8 relative tolerance for x &gt; 0.</p>
  631.      *
  632.      * <p>Performance for large negative values of x will be quite expensive (proportional to
  633.      * |x|).  Accuracy for negative values of x should be about 10^-8 absolute for results
  634.      * less than 10^5 and 10^-8 relative for results larger than that.</p>
  635.      *
  636.      * @param x Argument.
  637.      * @param <T> Type of the field elements.
  638.      * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller.
  639.      * @see <a href="http://en.wikipedia.org/wiki/Digamma_function">Digamma</a>
  640.      * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf">Bernardo&#39;s original article </a>
  641.      */
  642.     public static <T extends CalculusFieldElement<T>> T digamma(T x) {
  643.         if (x.isNaN() || x.isInfinite()) {
  644.             return x;
  645.         }

  646.         if (x.getReal() > 0 && x.getReal() <= S_LIMIT) {
  647.             // use method 5 from Bernardo AS103
  648.             // accurate to O(x)
  649.             return x.pow(-1).negate().subtract(GAMMA);
  650.         }

  651.         if (x.getReal() >= C_LIMIT) {
  652.             // use method 8 (accurate to O(1/x^8))
  653.             T inv = x.square().reciprocal();
  654.             //            1       1        1         1         1         5           691         1
  655.             // log(x) -  --- - ------ + ------- - ------- + ------- - ------- +  ---------- - -------
  656.             //           2 x   12 x^2   120 x^4   252 x^6   240 x^8   660 x^10   32760 x^12   12 x^14
  657.             return x.log().subtract(x.pow(-1).multiply(0.5)).add(
  658.                     inv.multiply(
  659.                             inv.multiply(
  660.                                     inv.multiply(
  661.                                             inv.multiply(
  662.                                                     inv.multiply(
  663.                                                             inv.multiply(inv.divide(-12.)
  664.                                                                             .add(691. / 32760))
  665.                                                                .subtract(5. / 660))
  666.                                                        .add(1.0 / 240))
  667.                                                .subtract(1.0 / 252))
  668.                                        .add(1.0 / 120))
  669.                                .subtract(1.0 / 12)));
  670.         }

  671.         return digamma(x.add(1.)).subtract(x.pow(-1));
  672.     }

  673.     /**
  674.      * Computes the trigamma function of x.
  675.      * This function is derived by taking the derivative of the implementation
  676.      * of digamma.
  677.      *
  678.      * @param x Argument.
  679.      * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller
  680.      * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma</a>
  681.      * @see Gamma#digamma(double)
  682.      */
  683.     public static double trigamma(double x) {
  684.         if (Double.isNaN(x) || Double.isInfinite(x)) {
  685.             return x;
  686.         }

  687.         if (x > 0 && x <= S_LIMIT) {
  688.             return 1 / (x * x);
  689.         }

  690.         if (x >= C_LIMIT) {
  691.             double inv = 1 / (x * x);
  692.             //  1    1      1       1       1      1         5        691        7
  693.             //  - + ---- + ---- - ----- + ----- - ----- + ------- - -------- + ------
  694.             //  x      2      3       5       7       9        11         13       15
  695.             //      2 x    6 x    30 x    42 x    30 x    66 x      2730 x      6 x
  696.             return 1 / x + inv * 0.5 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv * (1.0 / 42 - inv * (1.0 / 30 + inv *
  697.                     (5.0 / 66 - inv * (691. / 2730 + inv * 7. / 15))))));
  698.         }

  699.         return trigamma(x + 1) + 1 / (x * x);
  700.     }

  701.     /**
  702.      * Computes the trigamma function of x.
  703.      * This function is derived by taking the derivative of the implementation
  704.      * of digamma.
  705.      *
  706.      * @param x Argument.
  707.      * @param <T> Type of the field elements.
  708.      * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller
  709.      * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma</a>
  710.      * @see Gamma#digamma(double)
  711.      */
  712.     public static <T extends CalculusFieldElement<T>> T trigamma(T x) {
  713.         if (x.isNaN() || x.isInfinite()) {
  714.             return x;
  715.         }

  716.         if (x.getReal() > 0 && x.getReal() <= S_LIMIT) {
  717.             // use method 5 from Bernardo AS103
  718.             // accurate to O(x)
  719.             return x.square().reciprocal();
  720.         }

  721.         if (x.getReal() >= C_LIMIT) {
  722.             // use method 4 (accurate to O(1/x^8)
  723.             T inv    = x.square().reciprocal();
  724.             T invCub = inv.multiply(x.reciprocal());
  725.             //  1    1      1       1       1      1         5        691        7
  726.             //  - + ---- + ---- - ----- + ----- + ----- + ------- - -------- + ------
  727.             //  x      2      3       5       7       9        11         13       15
  728.             //      2 x    6 x    30 x    42 x    30 x    66 x      2730 x      6 x
  729.             return x.pow(-1).add(
  730.                     inv.multiply(0.5)).add(
  731.                             invCub.multiply(
  732.                                     inv.multiply(
  733.                                             inv.multiply(
  734.                                                     inv.multiply(
  735.                                                             inv.multiply(
  736.                                                                     inv.multiply(inv.multiply(7. / 6)
  737.                                                                                     .subtract(691. / 2730))
  738.                                                                        .add(5. / 66))
  739.                                                                .subtract(1.0 / 30))
  740.                                                        .add(1.0 / 42))
  741.                                                .subtract(1.0 / 30))
  742.                                        .add(1.0 / 6)));
  743.         }

  744.         return trigamma(x.add(1.)).add(x.square().reciprocal());
  745.     }

  746.     /**
  747.      * <p>
  748.      * Returns the Lanczos approximation used to compute the gamma function.
  749.      * The Lanczos approximation is related to the Gamma function by the
  750.      * following equation
  751.      * \[
  752.      * \Gamma(x) = \frac{\sqrt{2\pi}}{x} \times (x + g + \frac{1}{2}) ^ (x + \frac{1}{2})
  753.      *                   \times e^{-x - g - 0.5} \times \mathrm{lanczos}(x)
  754.      * \]
  755.      * where {@code g} is the Lanczos constant.
  756.      * </p>
  757.      *
  758.      * @param x Argument.
  759.      * @return The Lanczos approximation.
  760.      * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos Approximation</a>
  761.      * equations (1) through (5), and Paul Godfrey's
  762.      * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation
  763.      * of the convergent Lanczos complex Gamma approximation</a>
  764.      */
  765.     public static double lanczos(final double x) {
  766.         double sum = 0.0;
  767.         for (int i = LANCZOS.length - 1; i > 0; --i) {
  768.             sum += LANCZOS[i] / (x + i);
  769.         }
  770.         return sum + LANCZOS[0];
  771.     }

  772.     /**
  773.      * <p>
  774.      * Returns the Lanczos approximation used to compute the gamma function.
  775.      * The Lanczos approximation is related to the Gamma function by the
  776.      * following equation
  777.      * \[
  778.      * \Gamma(x) = \frac{\sqrt{2\pi}}{x} \times (x + g + \frac{1}{2}) ^ (x + \frac{1}{2})
  779.      *                   \times e^{-x - g - 0.5} \times \mathrm{lanczos}(x)
  780.      * \]
  781.      * where {@code g} is the Lanczos constant.
  782.      * </p>
  783.      *
  784.      * @param x Argument.
  785.      * @param <T> Type of the field elements.
  786.      * @return The Lanczos approximation.
  787.      * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos Approximation</a>
  788.      * equations (1) through (5), and Paul Godfrey's
  789.      * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation
  790.      * of the convergent Lanczos complex Gamma approximation</a>
  791.      */
  792.     public static <T extends CalculusFieldElement<T>> T lanczos(final T x) {
  793.         final Field<T> field = x.getField();
  794.         T              sum   = field.getZero();
  795.         for (int i = LANCZOS.length - 1; i > 0; --i) {
  796.             sum = sum.add(x.add(i).pow(-1.).multiply(LANCZOS[i]));
  797.         }
  798.         return sum.add(LANCZOS[0]);
  799.     }

  800.     /**
  801.      * Returns the value of 1 / &Gamma;(1 + x) - 1 for -0&#46;5 &le; x &le;
  802.      * 1&#46;5. This implementation is based on the double precision
  803.      * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
  804.      * {@code DGAM1}.
  805.      *
  806.      * @param x Argument.
  807.      * @return The value of {@code 1.0 / Gamma(1.0 + x) - 1.0}.
  808.      * @throws MathIllegalArgumentException if {@code x < -0.5}
  809.      * @throws MathIllegalArgumentException if {@code x > 1.5}
  810.      */
  811.     public static double invGamma1pm1(final double x) {

  812.         if (x < -0.5) {
  813.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
  814.                                                    x, -0.5);
  815.         }
  816.         if (x > 1.5) {
  817.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
  818.                                                    x, 1.5);
  819.         }

  820.         final double ret;
  821.         final double t = x <= 0.5 ? x : (x - 0.5) - 0.5;
  822.         if (t < 0.0) {
  823.             final double a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1;
  824.             double b = INV_GAMMA1P_M1_B8;
  825.             b = INV_GAMMA1P_M1_B7 + t * b;
  826.             b = INV_GAMMA1P_M1_B6 + t * b;
  827.             b = INV_GAMMA1P_M1_B5 + t * b;
  828.             b = INV_GAMMA1P_M1_B4 + t * b;
  829.             b = INV_GAMMA1P_M1_B3 + t * b;
  830.             b = INV_GAMMA1P_M1_B2 + t * b;
  831.             b = INV_GAMMA1P_M1_B1 + t * b;
  832.             b = 1.0 + t * b;

  833.             double c = INV_GAMMA1P_M1_C13 + t * (a / b);
  834.             c = INV_GAMMA1P_M1_C12 + t * c;
  835.             c = INV_GAMMA1P_M1_C11 + t * c;
  836.             c = INV_GAMMA1P_M1_C10 + t * c;
  837.             c = INV_GAMMA1P_M1_C9 + t * c;
  838.             c = INV_GAMMA1P_M1_C8 + t * c;
  839.             c = INV_GAMMA1P_M1_C7 + t * c;
  840.             c = INV_GAMMA1P_M1_C6 + t * c;
  841.             c = INV_GAMMA1P_M1_C5 + t * c;
  842.             c = INV_GAMMA1P_M1_C4 + t * c;
  843.             c = INV_GAMMA1P_M1_C3 + t * c;
  844.             c = INV_GAMMA1P_M1_C2 + t * c;
  845.             c = INV_GAMMA1P_M1_C1 + t * c;
  846.             c = INV_GAMMA1P_M1_C + t * c;
  847.             if (x > 0.5) {
  848.                 ret = t * c / x;
  849.             } else {
  850.                 ret = x * ((c + 0.5) + 0.5);
  851.             }
  852.         } else {
  853.             double p = INV_GAMMA1P_M1_P6;
  854.             p = INV_GAMMA1P_M1_P5 + t * p;
  855.             p = INV_GAMMA1P_M1_P4 + t * p;
  856.             p = INV_GAMMA1P_M1_P3 + t * p;
  857.             p = INV_GAMMA1P_M1_P2 + t * p;
  858.             p = INV_GAMMA1P_M1_P1 + t * p;
  859.             p = INV_GAMMA1P_M1_P0 + t * p;

  860.             double q = INV_GAMMA1P_M1_Q4;
  861.             q = INV_GAMMA1P_M1_Q3 + t * q;
  862.             q = INV_GAMMA1P_M1_Q2 + t * q;
  863.             q = INV_GAMMA1P_M1_Q1 + t * q;
  864.             q = 1.0 + t * q;

  865.             double c = INV_GAMMA1P_M1_C13 + (p / q) * t;
  866.             c = INV_GAMMA1P_M1_C12 + t * c;
  867.             c = INV_GAMMA1P_M1_C11 + t * c;
  868.             c = INV_GAMMA1P_M1_C10 + t * c;
  869.             c = INV_GAMMA1P_M1_C9 + t * c;
  870.             c = INV_GAMMA1P_M1_C8 + t * c;
  871.             c = INV_GAMMA1P_M1_C7 + t * c;
  872.             c = INV_GAMMA1P_M1_C6 + t * c;
  873.             c = INV_GAMMA1P_M1_C5 + t * c;
  874.             c = INV_GAMMA1P_M1_C4 + t * c;
  875.             c = INV_GAMMA1P_M1_C3 + t * c;
  876.             c = INV_GAMMA1P_M1_C2 + t * c;
  877.             c = INV_GAMMA1P_M1_C1 + t * c;
  878.             c = INV_GAMMA1P_M1_C0 + t * c;

  879.             if (x > 0.5) {
  880.                 ret = (t / x) * ((c - 0.5) - 0.5);
  881.             } else {
  882.                 ret = x * c;
  883.             }
  884.         }

  885.         return ret;
  886.     }

  887.     /**
  888.      * Returns the value of 1 / &Gamma;(1 + x) - 1 for -0&#46;5 &le; x &le;
  889.      * 1&#46;5. This implementation is based on the double precision
  890.      * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
  891.      * {@code DGAM1}.
  892.      *
  893.      * @param x Argument.
  894.      * @param <T> Type of the field elements.
  895.      * @return The value of {@code 1.0 / Gamma(1.0 + x) - 1.0}.
  896.      * @throws MathIllegalArgumentException if {@code x < -0.5}
  897.      * @throws MathIllegalArgumentException if {@code x > 1.5}
  898.      */
  899.     public static <T extends CalculusFieldElement<T>> T invGamma1pm1(final T x) {
  900.         final T one = x.getField().getOne();

  901.         if (x.getReal() < -0.5) {
  902.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
  903.                                                    x, -0.5);
  904.         }
  905.         if (x.getReal() > 1.5) {
  906.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
  907.                                                    x, 1.5);
  908.         }

  909.         final T ret;
  910.         final T t = x.getReal() <= 0.5 ? x : x.subtract(1);
  911.         if (t.getReal() < 0.0) {
  912.             final T a = one.newInstance(INV_GAMMA1P_M1_A0).add(t.multiply(INV_GAMMA1P_M1_A1));
  913.             T       b = one.newInstance(INV_GAMMA1P_M1_B8);
  914.             b = t.multiply(b).add(INV_GAMMA1P_M1_B7);
  915.             b = t.multiply(b).add(INV_GAMMA1P_M1_B6);
  916.             b = t.multiply(b).add(INV_GAMMA1P_M1_B5);
  917.             b = t.multiply(b).add(INV_GAMMA1P_M1_B4);
  918.             b = t.multiply(b).add(INV_GAMMA1P_M1_B3);
  919.             b = t.multiply(b).add(INV_GAMMA1P_M1_B2);
  920.             b = t.multiply(b).add(INV_GAMMA1P_M1_B1);
  921.             b = t.multiply(b).add(1.);

  922.             T c = one.newInstance(INV_GAMMA1P_M1_C13).add(t.multiply(a.divide(b)));
  923.             c = t.multiply(c).add(INV_GAMMA1P_M1_C12);
  924.             c = t.multiply(c).add(INV_GAMMA1P_M1_C11);
  925.             c = t.multiply(c).add(INV_GAMMA1P_M1_C10);
  926.             c = t.multiply(c).add(INV_GAMMA1P_M1_C9);
  927.             c = t.multiply(c).add(INV_GAMMA1P_M1_C8);
  928.             c = t.multiply(c).add(INV_GAMMA1P_M1_C7);
  929.             c = t.multiply(c).add(INV_GAMMA1P_M1_C6);
  930.             c = t.multiply(c).add(INV_GAMMA1P_M1_C5);
  931.             c = t.multiply(c).add(INV_GAMMA1P_M1_C4);
  932.             c = t.multiply(c).add(INV_GAMMA1P_M1_C3);
  933.             c = t.multiply(c).add(INV_GAMMA1P_M1_C2);
  934.             c = t.multiply(c).add(INV_GAMMA1P_M1_C1);
  935.             c = t.multiply(c).add(INV_GAMMA1P_M1_C);
  936.             if (x.getReal() > 0.5) {
  937.                 ret = t.multiply(c).divide(x);
  938.             }
  939.             else {
  940.                 ret = x.multiply(c.add(1));
  941.             }
  942.         }
  943.         else {
  944.             T p = one.newInstance(INV_GAMMA1P_M1_P6);
  945.             p = t.multiply(p).add(INV_GAMMA1P_M1_P5);
  946.             p = t.multiply(p).add(INV_GAMMA1P_M1_P4);
  947.             p = t.multiply(p).add(INV_GAMMA1P_M1_P3);
  948.             p = t.multiply(p).add(INV_GAMMA1P_M1_P2);
  949.             p = t.multiply(p).add(INV_GAMMA1P_M1_P1);
  950.             p = t.multiply(p).add(INV_GAMMA1P_M1_P0);

  951.             T q = one.newInstance(INV_GAMMA1P_M1_Q4);
  952.             q = t.multiply(q).add(INV_GAMMA1P_M1_Q3);
  953.             q = t.multiply(q).add(INV_GAMMA1P_M1_Q2);
  954.             q = t.multiply(q).add(INV_GAMMA1P_M1_Q1);
  955.             q = t.multiply(q).add(1.);

  956.             T c = one.newInstance(INV_GAMMA1P_M1_C13).add(t.multiply(p.divide(q)));
  957.             c = t.multiply(c).add(INV_GAMMA1P_M1_C12);
  958.             c = t.multiply(c).add(INV_GAMMA1P_M1_C11);
  959.             c = t.multiply(c).add(INV_GAMMA1P_M1_C10);
  960.             c = t.multiply(c).add(INV_GAMMA1P_M1_C9);
  961.             c = t.multiply(c).add(INV_GAMMA1P_M1_C8);
  962.             c = t.multiply(c).add(INV_GAMMA1P_M1_C7);
  963.             c = t.multiply(c).add(INV_GAMMA1P_M1_C6);
  964.             c = t.multiply(c).add(INV_GAMMA1P_M1_C5);
  965.             c = t.multiply(c).add(INV_GAMMA1P_M1_C4);
  966.             c = t.multiply(c).add(INV_GAMMA1P_M1_C3);
  967.             c = t.multiply(c).add(INV_GAMMA1P_M1_C2);
  968.             c = t.multiply(c).add(INV_GAMMA1P_M1_C1);
  969.             c = t.multiply(c).add(INV_GAMMA1P_M1_C0);

  970.             if (x.getReal() > 0.5) {
  971.                 ret = t.divide(x).multiply(c.subtract(1));
  972.             }
  973.             else {
  974.                 ret = x.multiply(c);
  975.             }
  976.         }

  977.         return ret;
  978.     }

  979.     /**
  980.      * Returns the value of log &Gamma;(1 + x) for -0&#46;5 &le; x &le; 1&#46;5.
  981.      * This implementation is based on the double precision implementation in
  982.      * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}.
  983.      *
  984.      * @param x Argument.
  985.      * @return The value of {@code log(Gamma(1 + x))}.
  986.      * @throws MathIllegalArgumentException if {@code x < -0.5}.
  987.      * @throws MathIllegalArgumentException if {@code x > 1.5}.
  988.      */
  989.     public static double logGamma1p(final double x)
  990.         throws MathIllegalArgumentException {

  991.         if (x < -0.5) {
  992.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
  993.                                                    x, -0.5);
  994.         }
  995.         if (x > 1.5) {
  996.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
  997.                                                    x, 1.5);
  998.         }

  999.         return -FastMath.log1p(invGamma1pm1(x));
  1000.     }

  1001.     /**
  1002.      * Returns the value of log &Gamma;(1 + x) for -0&#46;5 &le; x &le; 1&#46;5.
  1003.      * This implementation is based on the double precision implementation in
  1004.      * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}.
  1005.      *
  1006.      * @param x Argument.
  1007.      * @param <T> Type of the field elements.
  1008.      * @return The value of {@code log(Gamma(1 + x))}.
  1009.      * @throws MathIllegalArgumentException if {@code x < -0.5}.
  1010.      * @throws MathIllegalArgumentException if {@code x > 1.5}.
  1011.      */
  1012.     public static <T extends CalculusFieldElement<T>> T logGamma1p(final T x)
  1013.             throws MathIllegalArgumentException {

  1014.         if (x.getReal() < -0.5) {
  1015.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
  1016.                                                    x, -0.5);
  1017.         }
  1018.         if (x.getReal() > 1.5) {
  1019.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
  1020.                                                    x, 1.5);
  1021.         }

  1022.         return invGamma1pm1(x).log1p().negate();
  1023.     }


  1024.     /**
  1025.      * Returns the value of Γ(x). Based on the <em>NSWC Library of
  1026.      * Mathematics Subroutines</em> double precision implementation,
  1027.      * {@code DGAMMA}.
  1028.      *
  1029.      * @param x Argument.
  1030.      * @return the value of {@code Gamma(x)}.
  1031.      */
  1032.     public static double gamma(final double x) {

  1033.         if ((x == FastMath.rint(x)) && (x <= 0.0)) {
  1034.             return Double.NaN;
  1035.         }

  1036.         final double ret;
  1037.         final double absX = FastMath.abs(x);
  1038.         if (absX <= 20.0) {
  1039.             if (x >= 1.0) {
  1040.                 /*
  1041.                  * From the recurrence relation
  1042.                  * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
  1043.                  * then
  1044.                  * Gamma(t) = 1 / [1 + invGamma1pm1(t - 1)],
  1045.                  * where t = x - n. This means that t must satisfy
  1046.                  * -0.5 <= t - 1 <= 1.5.
  1047.                  */
  1048.                 double prod = 1.0;
  1049.                 double t = x;
  1050.                 while (t > 2.5) {
  1051.                     t -= 1.0;
  1052.                     prod *= t;
  1053.                 }
  1054.                 ret = prod / (1.0 + invGamma1pm1(t - 1.0));
  1055.             } else {
  1056.                 /*
  1057.                  * From the recurrence relation
  1058.                  * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
  1059.                  * then
  1060.                  * Gamma(x + n + 1) = 1 / [1 + invGamma1pm1(x + n)],
  1061.                  * which requires -0.5 <= x + n <= 1.5.
  1062.                  */
  1063.                 double prod = x;
  1064.                 double t = x;
  1065.                 while (t < -0.5) {
  1066.                     t += 1.0;
  1067.                     prod *= t;
  1068.                 }
  1069.                 ret = 1.0 / (prod * (1.0 + invGamma1pm1(t)));
  1070.             }
  1071.         } else {
  1072.             final double y = absX + LANCZOS_G + 0.5;
  1073.             final double gammaAbs = SQRT_TWO_PI / absX *
  1074.                                     FastMath.pow(y, absX + 0.5) *
  1075.                                     FastMath.exp(-y) * lanczos(absX);
  1076.             if (x > 0.0) {
  1077.                 ret = gammaAbs;
  1078.             } else {
  1079.                 /*
  1080.                  * From the reflection formula
  1081.                  * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
  1082.                  * and the recurrence relation
  1083.                  * Gamma(1 - x) = -x * Gamma(-x),
  1084.                  * it is found
  1085.                  * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
  1086.                  */
  1087.                 ret = -FastMath.PI /
  1088.                       (x * FastMath.sin(FastMath.PI * x) * gammaAbs);
  1089.             }
  1090.         }
  1091.         return ret;
  1092.     }

  1093.     /**
  1094.      * Returns the value of Γ(x). Based on the <em>NSWC Library of
  1095.      * Mathematics Subroutines</em> double precision implementation,
  1096.      * {@code DGAMMA}.
  1097.      *
  1098.      * @param x Argument.
  1099.      * @param <T> Type of the field elements.
  1100.      * @return the value of {@code Gamma(x)}.
  1101.      */
  1102.     public static <T extends CalculusFieldElement<T>> T gamma(final T x) {
  1103.         final T one = x.getField().getOne();

  1104.         if ((x.getReal() == x.rint().getReal()) && (x.getReal() <= 0.0)) {
  1105.             return one.multiply(Double.NaN);
  1106.         }

  1107.         final T ret;
  1108.         final T absX = x.abs();
  1109.         if (absX.getReal() <= 20.0) {
  1110.             if (x.getReal() >= 1.0) {
  1111.                 /*
  1112.                  * From the recurrence relation
  1113.                  * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
  1114.                  * then
  1115.                  * Gamma(t) = 1 / [1 + invGamma1pm1(t - 1)],
  1116.                  * where t = x - n. This means that t must satisfy
  1117.                  * -0.5 <= t - 1 <= 1.5.
  1118.                  */
  1119.                 T prod = one;
  1120.                 T t    = x;
  1121.                 while (t.getReal() > 2.5) {
  1122.                     t    = t.subtract(1.0);
  1123.                     prod = prod.multiply(t);
  1124.                 }
  1125.                 ret = prod.divide(invGamma1pm1(t.subtract(1.0)).add(1.0));
  1126.             }
  1127.             else {
  1128.                 /*
  1129.                  * From the recurrence relation
  1130.                  * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
  1131.                  * then
  1132.                  * Gamma(x + n + 1) = 1 / [1 + invGamma1pm1(x + n)],
  1133.                  * which requires -0.5 <= x + n <= 1.5.
  1134.                  */
  1135.                 T prod = x;
  1136.                 T t    = x;
  1137.                 while (t.getReal() < -0.5) {
  1138.                     t    = t.add(1.0);
  1139.                     prod = prod.multiply(t);
  1140.                 }
  1141.                 ret = prod.multiply(invGamma1pm1(t).add(1)).reciprocal();
  1142.             }
  1143.         }
  1144.         else {
  1145.             final T y = absX.add(LANCZOS_G + 0.5);
  1146.             final T gammaAbs = absX.reciprocal().multiply(SQRT_TWO_PI).multiply(y.pow(absX.add(0.5)))
  1147.                                    .multiply(y.negate().exp()).multiply(lanczos(absX));
  1148.             if (x.getReal() > 0.0) {
  1149.                 ret = gammaAbs;
  1150.             }
  1151.             else {
  1152.                 /*
  1153.                  * From the reflection formula
  1154.                  * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
  1155.                  * and the recurrence relation
  1156.                  * Gamma(1 - x) = -x * Gamma(-x),
  1157.                  * it is found
  1158.                  * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
  1159.                  */
  1160.                 ret = x.multiply(x.multiply(FastMath.PI).sin()).multiply(gammaAbs).reciprocal().multiply(-FastMath.PI);
  1161.             }
  1162.         }
  1163.         return ret;
  1164.     }
  1165. }