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 */
017package org.apache.commons.math3.special;
018
019import org.apache.commons.math3.exception.MaxCountExceededException;
020import org.apache.commons.math3.exception.NumberIsTooLargeException;
021import org.apache.commons.math3.exception.NumberIsTooSmallException;
022import org.apache.commons.math3.util.ContinuedFraction;
023import org.apache.commons.math3.util.FastMath;
024
025/**
026 * <p>
027 * This is a utility class that provides computation methods related to the
028 * &Gamma; (Gamma) family of functions.
029 * </p>
030 * <p>
031 * Implementation of {@link #invGamma1pm1(double)} and
032 * {@link #logGamma1p(double)} is based on the algorithms described in
033 * <ul>
034 * <li><a href="http://dx.doi.org/10.1145/22721.23109">Didonato and Morris
035 * (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios and
036 *     their Inverse</em>, TOMS 12(4), 377-393,</li>
037 * <li><a href="http://dx.doi.org/10.1145/131766.131776">Didonato and Morris
038 * (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the
039 *     Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li>
040 * </ul>
041 * and implemented in the
042 * <a href="http://www.dtic.mil/docs/citations/ADA476840">NSWC Library of Mathematical Functions</a>,
043 * available
044 * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html">here</a>.
045 * This library is "approved for public release", and the
046 * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf">Copyright guidance</a>
047 * indicates that unless otherwise stated in the code, all FORTRAN functions in
048 * this library are license free. Since no such notice appears in the code these
049 * functions can safely be ported to Commons-Math.
050 * </p>
051 *
052 */
053public class Gamma {
054    /**
055     * <a href="http://en.wikipedia.org/wiki/Euler-Mascheroni_constant">Euler-Mascheroni constant</a>
056     * @since 2.0
057     */
058    public static final double GAMMA = 0.577215664901532860606512090082;
059
060    /**
061     * The value of the {@code g} constant in the Lanczos approximation, see
062     * {@link #lanczos(double)}.
063     * @since 3.1
064     */
065    public static final double LANCZOS_G = 607.0 / 128.0;
066
067    /** Maximum allowed numerical error. */
068    private static final double DEFAULT_EPSILON = 10e-15;
069
070    /** Lanczos coefficients */
071    private static final double[] LANCZOS = {
072        0.99999999999999709182,
073        57.156235665862923517,
074        -59.597960355475491248,
075        14.136097974741747174,
076        -0.49191381609762019978,
077        .33994649984811888699e-4,
078        .46523628927048575665e-4,
079        -.98374475304879564677e-4,
080        .15808870322491248884e-3,
081        -.21026444172410488319e-3,
082        .21743961811521264320e-3,
083        -.16431810653676389022e-3,
084        .84418223983852743293e-4,
085        -.26190838401581408670e-4,
086        .36899182659531622704e-5,
087    };
088
089    /** Avoid repeated computation of log of 2 PI in logGamma */
090    private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(2.0 * FastMath.PI);
091
092    /** The constant value of &radic;(2&pi;). */
093    private static final double SQRT_TWO_PI = 2.506628274631000502;
094
095    // limits for switching algorithm in digamma
096    /** C limit. */
097    private static final double C_LIMIT = 49;
098
099    /** S limit. */
100    private static final double S_LIMIT = 1e-5;
101
102    /*
103     * Constants for the computation of double invGamma1pm1(double).
104     * Copied from DGAM1 in the NSWC library.
105     */
106
107    /** The constant {@code A0} defined in {@code DGAM1}. */
108    private static final double INV_GAMMA1P_M1_A0 = .611609510448141581788E-08;
109
110    /** The constant {@code A1} defined in {@code DGAM1}. */
111    private static final double INV_GAMMA1P_M1_A1 = .624730830116465516210E-08;
112
113    /** The constant {@code B1} defined in {@code DGAM1}. */
114    private static final double INV_GAMMA1P_M1_B1 = .203610414066806987300E+00;
115
116    /** The constant {@code B2} defined in {@code DGAM1}. */
117    private static final double INV_GAMMA1P_M1_B2 = .266205348428949217746E-01;
118
119    /** The constant {@code B3} defined in {@code DGAM1}. */
120    private static final double INV_GAMMA1P_M1_B3 = .493944979382446875238E-03;
121
122    /** The constant {@code B4} defined in {@code DGAM1}. */
123    private static final double INV_GAMMA1P_M1_B4 = -.851419432440314906588E-05;
124
125    /** The constant {@code B5} defined in {@code DGAM1}. */
126    private static final double INV_GAMMA1P_M1_B5 = -.643045481779353022248E-05;
127
128    /** The constant {@code B6} defined in {@code DGAM1}. */
129    private static final double INV_GAMMA1P_M1_B6 = .992641840672773722196E-06;
130
131    /** The constant {@code B7} defined in {@code DGAM1}. */
132    private static final double INV_GAMMA1P_M1_B7 = -.607761895722825260739E-07;
133
134    /** The constant {@code B8} defined in {@code DGAM1}. */
135    private static final double INV_GAMMA1P_M1_B8 = .195755836614639731882E-09;
136
137    /** The constant {@code P0} defined in {@code DGAM1}. */
138    private static final double INV_GAMMA1P_M1_P0 = .6116095104481415817861E-08;
139
140    /** The constant {@code P1} defined in {@code DGAM1}. */
141    private static final double INV_GAMMA1P_M1_P1 = .6871674113067198736152E-08;
142
143    /** The constant {@code P2} defined in {@code DGAM1}. */
144    private static final double INV_GAMMA1P_M1_P2 = .6820161668496170657918E-09;
145
146    /** The constant {@code P3} defined in {@code DGAM1}. */
147    private static final double INV_GAMMA1P_M1_P3 = .4686843322948848031080E-10;
148
149    /** The constant {@code P4} defined in {@code DGAM1}. */
150    private static final double INV_GAMMA1P_M1_P4 = .1572833027710446286995E-11;
151
152    /** The constant {@code P5} defined in {@code DGAM1}. */
153    private static final double INV_GAMMA1P_M1_P5 = -.1249441572276366213222E-12;
154
155    /** The constant {@code P6} defined in {@code DGAM1}. */
156    private static final double INV_GAMMA1P_M1_P6 = .4343529937408594255178E-14;
157
158    /** The constant {@code Q1} defined in {@code DGAM1}. */
159    private static final double INV_GAMMA1P_M1_Q1 = .3056961078365221025009E+00;
160
161    /** The constant {@code Q2} defined in {@code DGAM1}. */
162    private static final double INV_GAMMA1P_M1_Q2 = .5464213086042296536016E-01;
163
164    /** The constant {@code Q3} defined in {@code DGAM1}. */
165    private static final double INV_GAMMA1P_M1_Q3 = .4956830093825887312020E-02;
166
167    /** The constant {@code Q4} defined in {@code DGAM1}. */
168    private static final double INV_GAMMA1P_M1_Q4 = .2692369466186361192876E-03;
169
170    /** The constant {@code C} defined in {@code DGAM1}. */
171    private static final double INV_GAMMA1P_M1_C = -.422784335098467139393487909917598E+00;
172
173    /** The constant {@code C0} defined in {@code DGAM1}. */
174    private static final double INV_GAMMA1P_M1_C0 = .577215664901532860606512090082402E+00;
175
176    /** The constant {@code C1} defined in {@code DGAM1}. */
177    private static final double INV_GAMMA1P_M1_C1 = -.655878071520253881077019515145390E+00;
178
179    /** The constant {@code C2} defined in {@code DGAM1}. */
180    private static final double INV_GAMMA1P_M1_C2 = -.420026350340952355290039348754298E-01;
181
182    /** The constant {@code C3} defined in {@code DGAM1}. */
183    private static final double INV_GAMMA1P_M1_C3 = .166538611382291489501700795102105E+00;
184
185    /** The constant {@code C4} defined in {@code DGAM1}. */
186    private static final double INV_GAMMA1P_M1_C4 = -.421977345555443367482083012891874E-01;
187
188    /** The constant {@code C5} defined in {@code DGAM1}. */
189    private static final double INV_GAMMA1P_M1_C5 = -.962197152787697356211492167234820E-02;
190
191    /** The constant {@code C6} defined in {@code DGAM1}. */
192    private static final double INV_GAMMA1P_M1_C6 = .721894324666309954239501034044657E-02;
193
194    /** The constant {@code C7} defined in {@code DGAM1}. */
195    private static final double INV_GAMMA1P_M1_C7 = -.116516759185906511211397108401839E-02;
196
197    /** The constant {@code C8} defined in {@code DGAM1}. */
198    private static final double INV_GAMMA1P_M1_C8 = -.215241674114950972815729963053648E-03;
199
200    /** The constant {@code C9} defined in {@code DGAM1}. */
201    private static final double INV_GAMMA1P_M1_C9 = .128050282388116186153198626328164E-03;
202
203    /** The constant {@code C10} defined in {@code DGAM1}. */
204    private static final double INV_GAMMA1P_M1_C10 = -.201348547807882386556893914210218E-04;
205
206    /** The constant {@code C11} defined in {@code DGAM1}. */
207    private static final double INV_GAMMA1P_M1_C11 = -.125049348214267065734535947383309E-05;
208
209    /** The constant {@code C12} defined in {@code DGAM1}. */
210    private static final double INV_GAMMA1P_M1_C12 = .113302723198169588237412962033074E-05;
211
212    /** The constant {@code C13} defined in {@code DGAM1}. */
213    private static final double INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06;
214
215    /**
216     * Default constructor.  Prohibit instantiation.
217     */
218    private Gamma() {}
219
220    /**
221     * <p>
222     * Returns the value of log&nbsp;&Gamma;(x) for x&nbsp;&gt;&nbsp;0.
223     * </p>
224     * <p>
225     * For x &le; 8, the implementation is based on the double precision
226     * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
227     * {@code DGAMLN}. For x &gt; 8, the implementation is based on
228     * </p>
229     * <ul>
230     * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
231     *     Function</a>, equation (28).</li>
232     * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
233     *     Lanczos Approximation</a>, equations (1) through (5).</li>
234     * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
235     *     the computation of the convergent Lanczos complex Gamma
236     *     approximation</a></li>
237     * </ul>
238     *
239     * @param x Argument.
240     * @return the value of {@code log(Gamma(x))}, {@code Double.NaN} if
241     * {@code x <= 0.0}.
242     */
243    public static double logGamma(double x) {
244        double ret;
245
246        if (Double.isNaN(x) || (x <= 0.0)) {
247            ret = Double.NaN;
248        } else if (x < 0.5) {
249            return logGamma1p(x) - FastMath.log(x);
250        } else if (x <= 2.5) {
251            return logGamma1p((x - 0.5) - 0.5);
252        } else if (x <= 8.0) {
253            final int n = (int) FastMath.floor(x - 1.5);
254            double prod = 1.0;
255            for (int i = 1; i <= n; i++) {
256                prod *= x - i;
257            }
258            return logGamma1p(x - (n + 1)) + FastMath.log(prod);
259        } else {
260            double sum = lanczos(x);
261            double tmp = x + LANCZOS_G + .5;
262            ret = ((x + .5) * FastMath.log(tmp)) - tmp +
263                HALF_LOG_2_PI + FastMath.log(sum / x);
264        }
265
266        return ret;
267    }
268
269    /**
270     * Returns the regularized gamma function P(a, x).
271     *
272     * @param a Parameter.
273     * @param x Value.
274     * @return the regularized gamma function P(a, x).
275     * @throws MaxCountExceededException if the algorithm fails to converge.
276     */
277    public static double regularizedGammaP(double a, double x) {
278        return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
279    }
280
281    /**
282     * Returns the regularized gamma function P(a, x).
283     *
284     * The implementation of this method is based on:
285     * <ul>
286     *  <li>
287     *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
288     *   Regularized Gamma Function</a>, equation (1)
289     *  </li>
290     *  <li>
291     *   <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
292     *   Incomplete Gamma Function</a>, equation (4).
293     *  </li>
294     *  <li>
295     *   <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html">
296     *   Confluent Hypergeometric Function of the First Kind</a>, equation (1).
297     *  </li>
298     * </ul>
299     *
300     * @param a the a parameter.
301     * @param x the value.
302     * @param epsilon When the absolute value of the nth item in the
303     * series is less than epsilon the approximation ceases to calculate
304     * further elements in the series.
305     * @param maxIterations Maximum number of "iterations" to complete.
306     * @return the regularized gamma function P(a, x)
307     * @throws MaxCountExceededException if the algorithm fails to converge.
308     */
309    public static double regularizedGammaP(double a,
310                                           double x,
311                                           double epsilon,
312                                           int maxIterations) {
313        double ret;
314
315        if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
316            ret = Double.NaN;
317        } else if (x == 0.0) {
318            ret = 0.0;
319        } else if (x >= a + 1) {
320            // use regularizedGammaQ because it should converge faster in this
321            // case.
322            ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
323        } else {
324            // calculate series
325            double n = 0.0; // current element index
326            double an = 1.0 / a; // n-th element in the series
327            double sum = an; // partial sum
328            while (FastMath.abs(an/sum) > epsilon &&
329                   n < maxIterations &&
330                   sum < Double.POSITIVE_INFINITY) {
331                // compute next element in the series
332                n += 1.0;
333                an *= x / (a + n);
334
335                // update partial sum
336                sum += an;
337            }
338            if (n >= maxIterations) {
339                throw new MaxCountExceededException(maxIterations);
340            } else if (Double.isInfinite(sum)) {
341                ret = 1.0;
342            } else {
343                ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * sum;
344            }
345        }
346
347        return ret;
348    }
349
350    /**
351     * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
352     *
353     * @param a the a parameter.
354     * @param x the value.
355     * @return the regularized gamma function Q(a, x)
356     * @throws MaxCountExceededException if the algorithm fails to converge.
357     */
358    public static double regularizedGammaQ(double a, double x) {
359        return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
360    }
361
362    /**
363     * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
364     *
365     * The implementation of this method is based on:
366     * <ul>
367     *  <li>
368     *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
369     *   Regularized Gamma Function</a>, equation (1).
370     *  </li>
371     *  <li>
372     *   <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
373     *   Regularized incomplete gamma function: Continued fraction representations
374     *   (formula 06.08.10.0003)</a>
375     *  </li>
376     * </ul>
377     *
378     * @param a the a parameter.
379     * @param x the value.
380     * @param epsilon When the absolute value of the nth item in the
381     * series is less than epsilon the approximation ceases to calculate
382     * further elements in the series.
383     * @param maxIterations Maximum number of "iterations" to complete.
384     * @return the regularized gamma function P(a, x)
385     * @throws MaxCountExceededException if the algorithm fails to converge.
386     */
387    public static double regularizedGammaQ(final double a,
388                                           double x,
389                                           double epsilon,
390                                           int maxIterations) {
391        double ret;
392
393        if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
394            ret = Double.NaN;
395        } else if (x == 0.0) {
396            ret = 1.0;
397        } else if (x < a + 1.0) {
398            // use regularizedGammaP because it should converge faster in this
399            // case.
400            ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations);
401        } else {
402            // create continued fraction
403            ContinuedFraction cf = new ContinuedFraction() {
404
405                @Override
406                protected double getA(int n, double x) {
407                    return ((2.0 * n) + 1.0) - a + x;
408                }
409
410                @Override
411                protected double getB(int n, double x) {
412                    return n * (a - n);
413                }
414            };
415
416            ret = 1.0 / cf.evaluate(x, epsilon, maxIterations);
417            ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * ret;
418        }
419
420        return ret;
421    }
422
423
424    /**
425     * <p>Computes the digamma function of x.</p>
426     *
427     * <p>This is an independently written implementation of the algorithm described in
428     * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p>
429     *
430     * <p>Some of the constants have been changed to increase accuracy at the moderate expense
431     * of run-time.  The result should be accurate to within 10^-8 absolute tolerance for
432     * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.</p>
433     *
434     * <p>Performance for large negative values of x will be quite expensive (proportional to
435     * |x|).  Accuracy for negative values of x should be about 10^-8 absolute for results
436     * less than 10^5 and 10^-8 relative for results larger than that.</p>
437     *
438     * @param x Argument.
439     * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller.
440     * @see <a href="http://en.wikipedia.org/wiki/Digamma_function">Digamma</a>
441     * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf">Bernardo&apos;s original article </a>
442     * @since 2.0
443     */
444    public static double digamma(double x) {
445        if (x > 0 && x <= S_LIMIT) {
446            // use method 5 from Bernardo AS103
447            // accurate to O(x)
448            return -GAMMA - 1 / x;
449        }
450
451        if (x >= C_LIMIT) {
452            // use method 4 (accurate to O(1/x^8)
453            double inv = 1 / (x * x);
454            //            1       1        1         1
455            // log(x) -  --- - ------ + ------- - -------
456            //           2 x   12 x^2   120 x^4   252 x^6
457            return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252));
458        }
459
460        return digamma(x + 1) - 1 / x;
461    }
462
463    /**
464     * Computes the trigamma function of x.
465     * This function is derived by taking the derivative of the implementation
466     * of digamma.
467     *
468     * @param x Argument.
469     * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller
470     * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma</a>
471     * @see Gamma#digamma(double)
472     * @since 2.0
473     */
474    public static double trigamma(double x) {
475        if (x > 0 && x <= S_LIMIT) {
476            return 1 / (x * x);
477        }
478
479        if (x >= C_LIMIT) {
480            double inv = 1 / (x * x);
481            //  1    1      1       1       1
482            //  - + ---- + ---- - ----- + -----
483            //  x      2      3       5       7
484            //      2 x    6 x    30 x    42 x
485            return 1 / x + inv / 2 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv / 42));
486        }
487
488        return trigamma(x + 1) + 1 / (x * x);
489    }
490
491    /**
492     * <p>
493     * Returns the Lanczos approximation used to compute the gamma function.
494     * The Lanczos approximation is related to the Gamma function by the
495     * following equation
496     * <center>
497     * {@code gamma(x) = sqrt(2 * pi) / x * (x + g + 0.5) ^ (x + 0.5)
498     *                   * exp(-x - g - 0.5) * lanczos(x)},
499     * </center>
500     * where {@code g} is the Lanczos constant.
501     * </p>
502     *
503     * @param x Argument.
504     * @return The Lanczos approximation.
505     * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos Approximation</a>
506     * equations (1) through (5), and Paul Godfrey's
507     * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation
508     * of the convergent Lanczos complex Gamma approximation</a>
509     * @since 3.1
510     */
511    public static double lanczos(final double x) {
512        double sum = 0.0;
513        for (int i = LANCZOS.length - 1; i > 0; --i) {
514            sum += LANCZOS[i] / (x + i);
515        }
516        return sum + LANCZOS[0];
517    }
518
519    /**
520     * Returns the value of 1 / &Gamma;(1 + x) - 1 for -0&#46;5 &le; x &le;
521     * 1&#46;5. This implementation is based on the double precision
522     * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
523     * {@code DGAM1}.
524     *
525     * @param x Argument.
526     * @return The value of {@code 1.0 / Gamma(1.0 + x) - 1.0}.
527     * @throws NumberIsTooSmallException if {@code x < -0.5}
528     * @throws NumberIsTooLargeException if {@code x > 1.5}
529     * @since 3.1
530     */
531    public static double invGamma1pm1(final double x) {
532
533        if (x < -0.5) {
534            throw new NumberIsTooSmallException(x, -0.5, true);
535        }
536        if (x > 1.5) {
537            throw new NumberIsTooLargeException(x, 1.5, true);
538        }
539
540        final double ret;
541        final double t = x <= 0.5 ? x : (x - 0.5) - 0.5;
542        if (t < 0.0) {
543            final double a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1;
544            double b = INV_GAMMA1P_M1_B8;
545            b = INV_GAMMA1P_M1_B7 + t * b;
546            b = INV_GAMMA1P_M1_B6 + t * b;
547            b = INV_GAMMA1P_M1_B5 + t * b;
548            b = INV_GAMMA1P_M1_B4 + t * b;
549            b = INV_GAMMA1P_M1_B3 + t * b;
550            b = INV_GAMMA1P_M1_B2 + t * b;
551            b = INV_GAMMA1P_M1_B1 + t * b;
552            b = 1.0 + t * b;
553
554            double c = INV_GAMMA1P_M1_C13 + t * (a / b);
555            c = INV_GAMMA1P_M1_C12 + t * c;
556            c = INV_GAMMA1P_M1_C11 + t * c;
557            c = INV_GAMMA1P_M1_C10 + t * c;
558            c = INV_GAMMA1P_M1_C9 + t * c;
559            c = INV_GAMMA1P_M1_C8 + t * c;
560            c = INV_GAMMA1P_M1_C7 + t * c;
561            c = INV_GAMMA1P_M1_C6 + t * c;
562            c = INV_GAMMA1P_M1_C5 + t * c;
563            c = INV_GAMMA1P_M1_C4 + t * c;
564            c = INV_GAMMA1P_M1_C3 + t * c;
565            c = INV_GAMMA1P_M1_C2 + t * c;
566            c = INV_GAMMA1P_M1_C1 + t * c;
567            c = INV_GAMMA1P_M1_C + t * c;
568            if (x > 0.5) {
569                ret = t * c / x;
570            } else {
571                ret = x * ((c + 0.5) + 0.5);
572            }
573        } else {
574            double p = INV_GAMMA1P_M1_P6;
575            p = INV_GAMMA1P_M1_P5 + t * p;
576            p = INV_GAMMA1P_M1_P4 + t * p;
577            p = INV_GAMMA1P_M1_P3 + t * p;
578            p = INV_GAMMA1P_M1_P2 + t * p;
579            p = INV_GAMMA1P_M1_P1 + t * p;
580            p = INV_GAMMA1P_M1_P0 + t * p;
581
582            double q = INV_GAMMA1P_M1_Q4;
583            q = INV_GAMMA1P_M1_Q3 + t * q;
584            q = INV_GAMMA1P_M1_Q2 + t * q;
585            q = INV_GAMMA1P_M1_Q1 + t * q;
586            q = 1.0 + t * q;
587
588            double c = INV_GAMMA1P_M1_C13 + (p / q) * t;
589            c = INV_GAMMA1P_M1_C12 + t * c;
590            c = INV_GAMMA1P_M1_C11 + t * c;
591            c = INV_GAMMA1P_M1_C10 + t * c;
592            c = INV_GAMMA1P_M1_C9 + t * c;
593            c = INV_GAMMA1P_M1_C8 + t * c;
594            c = INV_GAMMA1P_M1_C7 + t * c;
595            c = INV_GAMMA1P_M1_C6 + t * c;
596            c = INV_GAMMA1P_M1_C5 + t * c;
597            c = INV_GAMMA1P_M1_C4 + t * c;
598            c = INV_GAMMA1P_M1_C3 + t * c;
599            c = INV_GAMMA1P_M1_C2 + t * c;
600            c = INV_GAMMA1P_M1_C1 + t * c;
601            c = INV_GAMMA1P_M1_C0 + t * c;
602
603            if (x > 0.5) {
604                ret = (t / x) * ((c - 0.5) - 0.5);
605            } else {
606                ret = x * c;
607            }
608        }
609
610        return ret;
611    }
612
613    /**
614     * Returns the value of log &Gamma;(1 + x) for -0&#46;5 &le; x &le; 1&#46;5.
615     * This implementation is based on the double precision implementation in
616     * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}.
617     *
618     * @param x Argument.
619     * @return The value of {@code log(Gamma(1 + x))}.
620     * @throws NumberIsTooSmallException if {@code x < -0.5}.
621     * @throws NumberIsTooLargeException if {@code x > 1.5}.
622     * @since 3.1
623     */
624    public static double logGamma1p(final double x)
625        throws NumberIsTooSmallException, NumberIsTooLargeException {
626
627        if (x < -0.5) {
628            throw new NumberIsTooSmallException(x, -0.5, true);
629        }
630        if (x > 1.5) {
631            throw new NumberIsTooLargeException(x, 1.5, true);
632        }
633
634        return -FastMath.log1p(invGamma1pm1(x));
635    }
636
637
638    /**
639     * Returns the value of Γ(x). Based on the <em>NSWC Library of
640     * Mathematics Subroutines</em> double precision implementation,
641     * {@code DGAMMA}.
642     *
643     * @param x Argument.
644     * @return the value of {@code Gamma(x)}.
645     * @since 3.1
646     */
647    public static double gamma(final double x) {
648
649        if ((x == FastMath.rint(x)) && (x <= 0.0)) {
650            return Double.NaN;
651        }
652
653        final double ret;
654        final double absX = FastMath.abs(x);
655        if (absX <= 20.0) {
656            if (x >= 1.0) {
657                /*
658                 * From the recurrence relation
659                 * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
660                 * then
661                 * Gamma(t) = 1 / [1 + invGamma1pm1(t - 1)],
662                 * where t = x - n. This means that t must satisfy
663                 * -0.5 <= t - 1 <= 1.5.
664                 */
665                double prod = 1.0;
666                double t = x;
667                while (t > 2.5) {
668                    t -= 1.0;
669                    prod *= t;
670                }
671                ret = prod / (1.0 + invGamma1pm1(t - 1.0));
672            } else {
673                /*
674                 * From the recurrence relation
675                 * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
676                 * then
677                 * Gamma(x + n + 1) = 1 / [1 + invGamma1pm1(x + n)],
678                 * which requires -0.5 <= x + n <= 1.5.
679                 */
680                double prod = x;
681                double t = x;
682                while (t < -0.5) {
683                    t += 1.0;
684                    prod *= t;
685                }
686                ret = 1.0 / (prod * (1.0 + invGamma1pm1(t)));
687            }
688        } else {
689            final double y = absX + LANCZOS_G + 0.5;
690            final double gammaAbs = SQRT_TWO_PI / x *
691                                    FastMath.pow(y, absX + 0.5) *
692                                    FastMath.exp(-y) * lanczos(absX);
693            if (x > 0.0) {
694                ret = gammaAbs;
695            } else {
696                /*
697                 * From the reflection formula
698                 * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
699                 * and the recurrence relation
700                 * Gamma(1 - x) = -x * Gamma(-x),
701                 * it is found
702                 * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
703                 */
704                ret = -FastMath.PI /
705                      (x * FastMath.sin(FastMath.PI * x) * gammaAbs);
706            }
707        }
708        return ret;
709    }
710}