    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math.special;
    018
    019    import org.apache.commons.math.MathException;
    020    import org.apache.commons.math.MaxIterationsExceededException;
    021    import org.apache.commons.math.util.ContinuedFraction;
    022    import org.apache.commons.math.util.FastMath;
    023
    024    /**
    025     * This is a utility class that provides computation methods related to the
    026     * Gamma family of functions.
    027     *
    028     * @version $Revision: 1042510 $ $Date: 2010-12-06 02:54:18 +0100 (lun. 06 d??c. 2010) $
    029     */
    030    public class Gamma {
    031
    032        /**
    033         * <a href="http://en.wikipedia.org/wiki/Euler-Mascheroni_constant">Euler-Mascheroni constant</a>
    034         * @since 2.0
    035         */
    036        public static final double GAMMA = 0.577215664901532860606512090082;
    037
    038        /** Maximum allowed numerical error. */
    039        private static final double DEFAULT_EPSILON = 10e-15;
    040
    041        /** Lanczos coefficients */
    042        private static final double[] LANCZOS =
    043        {
    044            0.99999999999999709182,
    045            57.156235665862923517,
    046            -59.597960355475491248,
    047            14.136097974741747174,
    048            -0.49191381609762019978,
    049            .33994649984811888699e-4,
    050            .46523628927048575665e-4,
    051            -.98374475304879564677e-4,
    052            .15808870322491248884e-3,
    053            -.21026444172410488319e-3,
    054            .21743961811521264320e-3,
    055            -.16431810653676389022e-3,
    056            .84418223983852743293e-4,
    057            -.26190838401581408670e-4,
    058            .36899182659531622704e-5,
    059        };
    060
    061        /** Avoid repeated computation of log of 2 PI in logGamma */
    062        private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(2.0 * FastMath.PI);
    063
    064        // limits for switching algorithm in digamma
    065        /** C limit. */
    066        private static final double C_LIMIT = 49;
    067
    068        /** S limit. */
    069        private static final double S_LIMIT = 1e-5;
    070
    071        /**
    072         * Default constructor.  Prohibit instantiation.
    073         */
    074        private Gamma() {
    075            super();
    076        }
    077
    078        /**
    079         * Returns the natural logarithm of the gamma function &#915;(x).
    080         *
    081         * The implementation of this method is based on:
    082         * <ul>
    083         * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">
    084         * Gamma Function</a>, equation (28).</li>
    085         * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
    086         * Lanczos Approximation</a>, equations (1) through (5).</li>
    087         * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
    088         * the computation of the convergent Lanczos complex Gamma approximation
    089         * </a></li>
    090         * </ul>
    091         *
    092         * @param x the value.
    093         * @return log(&#915;(x))
    094         */
    095        public static double logGamma(double x) {
    096            double ret;
    097
    098            if (Double.isNaN(x) || (x <= 0.0)) {
    099                ret = Double.NaN;
    100            } else {
    101                double g = 607.0 / 128.0;
    102
    103                double sum = 0.0;
    104                for (int i = LANCZOS.length - 1; i > 0; --i) {
    105                    sum = sum + (LANCZOS[i] / (x + i));
    106                }
    107                sum = sum + LANCZOS[0];
    108
    109                double tmp = x + g + .5;
    110                ret = ((x + .5) * FastMath.log(tmp)) - tmp +
    111                    HALF_LOG_2_PI + FastMath.log(sum / x);
    112            }
    113
    114            return ret;
    115        }
    116
    117        /**
    118         * Returns the regularized gamma function P(a, x).
    119         *
    120         * @param a the a parameter.
    121         * @param x the value.
    122         * @return the regularized gamma function P(a, x)
    123         * @throws MathException if the algorithm fails to converge.
    124         */
    125        public static double regularizedGammaP(double a, double x)
    126            throws MathException
    127        {
    128            return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
    129        }
    130
    131
    132        /**
    133         * Returns the regularized gamma function P(a, x).
    134         *
    135         * The implementation of this method is based on:
    136         * <ul>
    137         * <li>
    138         * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
    139         * Regularized Gamma Function</a>, equation (1).</li>
    140         * <li>
    141         * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
    142         * Incomplete Gamma Function</a>, equation (4).</li>
    143         * <li>
    144         * <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html">
    145         * Confluent Hypergeometric Function of the First Kind</a>, equation (1).
    146         * </li>
    147         * </ul>
    148         *
    149         * @param a the a parameter.
    150         * @param x the value.
    151         * @param epsilon When the absolute value of the nth item in the
    152         *                series is less than epsilon the approximation ceases
    153         *                to calculate further elements in the series.
    154         * @param maxIterations Maximum number of "iterations" to complete.
    155         * @return the regularized gamma function P(a, x)
    156         * @throws MathException if the algorithm fails to converge.
    157         */
    158        public static double regularizedGammaP(double a,
    159                                               double x,
    160                                               double epsilon,
    161                                               int maxIterations)
    162            throws MathException
    163        {
    164            double ret;
    165
    166            if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
    167                ret = Double.NaN;
    168            } else if (x == 0.0) {
    169                ret = 0.0;
    170            } else if (x >= a + 1) {
    171                // use regularizedGammaQ because it should converge faster in this
    172                // case.
    173                ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
    174            } else {
    175                // calculate series
    176                double n = 0.0; // current element index
    177                double an = 1.0 / a; // n-th element in the series
    178                double sum = an; // partial sum
    179                while (FastMath.abs(an/sum) > epsilon && n < maxIterations && sum < Double.POSITIVE_INFINITY) {
    180                    // compute next element in the series
    181                    n = n + 1.0;
    182                    an = an * (x / (a + n));
    183
    184                    // update partial sum
    185                    sum = sum + an;
    186                }
    187                if (n >= maxIterations) {
    188                    throw new MaxIterationsExceededException(maxIterations);
    189                } else if (Double.isInfinite(sum)) {
    190                    ret = 1.0;
    191                } else {
    192                    ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * sum;
    193                }
    194            }
    195
    196            return ret;
    197        }
    198
    199        /**
    200         * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
    201         *
    202         * @param a the a parameter.
    203         * @param x the value.
    204         * @return the regularized gamma function Q(a, x)
    205         * @throws MathException if the algorithm fails to converge.
    206         */
    207        public static double regularizedGammaQ(double a, double x)
    208            throws MathException
    209        {
    210            return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
    211        }
    212
    213        /**
    214         * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
    215         *
    216         * The implementation of this method is based on:
    217         * <ul>
    218         * <li>
    219         * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
    220         * Regularized Gamma Function</a>, equation (1).</li>
    221         * <li>
    222         * <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
    223         * Regularized incomplete gamma function: Continued fraction representations  (formula 06.08.10.0003)</a></li>
    224         * </ul>
    225         *
    226         * @param a the a parameter.
    227         * @param x the value.
    228         * @param epsilon When the absolute value of the nth item in the
    229         *                series is less than epsilon the approximation ceases
    230         *                to calculate further elements in the series.
    231         * @param maxIterations Maximum number of "iterations" to complete.
    232         * @return the regularized gamma function P(a, x)
    233         * @throws MathException if the algorithm fails to converge.
    234         */
    235        public static double regularizedGammaQ(final double a,
    236                                               double x,
    237                                               double epsilon,
    238                                               int maxIterations)
    239            throws MathException
    240        {
    241            double ret;
    242
    243            if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
    244                ret = Double.NaN;
    245            } else if (x == 0.0) {
    246                ret = 1.0;
    247            } else if (x < a + 1.0) {
    248                // use regularizedGammaP because it should converge faster in this
    249                // case.
    250                ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations);
    251            } else {
    252                // create continued fraction
    253                ContinuedFraction cf = new ContinuedFraction() {
    254
    255                    @Override
    256                    protected double getA(int n, double x) {
    257                        return ((2.0 * n) + 1.0) - a + x;
    258                    }
    259
    260                    @Override
    261                    protected double getB(int n, double x) {
    262                        return n * (a - n);
    263                    }
    264                };
    265
    266                ret = 1.0 / cf.evaluate(x, epsilon, maxIterations);
    267                ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * ret;
    268            }
    269
    270            return ret;
    271        }
    272
    273
    274        /**
    275         * <p>Computes the digamma function of x.</p>
    276         *
    277         * <p>This is an independently written implementation of the algorithm described in
    278         * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p>
    279         *
    280         * <p>Some of the constants have been changed to increase accuracy at the moderate expense
    281         * of run-time.  The result should be accurate to within 10^-8 absolute tolerance for
    282         * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.</p>
    283         *
    284         * <p>Performance for large negative values of x will be quite expensive (proportional to
    285         * |x|).  Accuracy for negative values of x should be about 10^-8 absolute for results
    286         * less than 10^5 and 10^-8 relative for results larger than that.</p>
    287         *
    288         * @param x  the argument
    289         * @return   digamma(x) to within 10-8 relative or absolute error whichever is smaller
    290         * @see <a href="http://en.wikipedia.org/wiki/Digamma_function"> Digamma at wikipedia </a>
    291         * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf"> Bernardo&apos;s original article </a>
    292         * @since 2.0
    293         */
    294        public static double digamma(double x) {
    295            if (x > 0 && x <= S_LIMIT) {
    296                // use method 5 from Bernardo AS103
    297                // accurate to O(x)
    298                return -GAMMA - 1 / x;
    299            }
    300
    301            if (x >= C_LIMIT) {
    302                // use method 4 (accurate to O(1/x^8)
    303                double inv = 1 / (x * x);
    304                //            1       1        1         1
    305                // log(x) -  --- - ------ + ------- - -------
    306                //           2 x   12 x^2   120 x^4   252 x^6
    307                return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252));
    308            }
    309
    310            return digamma(x + 1) - 1 / x;
    311        }
    312
    313        /**
    314         * <p>Computes the trigamma function of x.  This function is derived by taking the derivative of
    315         * the implementation of digamma.</p>
    316         *
    317         * @param x  the argument
    318         * @return   trigamma(x) to within 10-8 relative or absolute error whichever is smaller
    319         * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function"> Trigamma at wikipedia </a>
    320         * @see Gamma#digamma(double)
    321         * @since 2.0
    322         */
    323        public static double trigamma(double x) {
    324            if (x > 0 && x <= S_LIMIT) {
    325                return 1 / (x * x);
    326            }
    327
    328            if (x >= C_LIMIT) {
    329                double inv = 1 / (x * x);
    330                //  1    1      1       1       1
    331                //  - + ---- + ---- - ----- + -----
    332                //  x      2      3       5       7
    333                //      2 x    6 x    30 x    42 x
    334                return 1 / x + inv / 2 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv / 42));
    335            }
    336
    337            return trigamma(x + 1) + 1 / (x * x);
    338        }
    339    }