View Javadoc
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  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.special;
23  
24  import org.hipparchus.CalculusFieldElement;
25  import org.hipparchus.Field;
26  import org.hipparchus.exception.LocalizedCoreFormats;
27  import org.hipparchus.exception.MathIllegalArgumentException;
28  import org.hipparchus.exception.MathIllegalStateException;
29  import org.hipparchus.util.ContinuedFraction;
30  import org.hipparchus.util.FastMath;
31  import org.hipparchus.util.FieldContinuedFraction;
32  
33  /**
34   * <p>
35   * This is a utility class that provides computation methods related to the
36   * &Gamma; (Gamma) family of functions.
37   * </p>
38   * <p>
39   * Implementation of {@link #invGamma1pm1(double)} and
40   * {@link #logGamma1p(double)} is based on the algorithms described in
41   * </p>
42   * <ul>
43   * <li><a href="http://dx.doi.org/10.1145/22721.23109">Didonato and Morris
44   * (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios and
45   *     their Inverse</em>, TOMS 12(4), 377-393,</li>
46   * <li><a href="http://dx.doi.org/10.1145/131766.131776">Didonato and Morris
47   * (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the
48   *     Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li>
49   * </ul>
50   * <p>
51   * and implemented in the
52   * <a href="http://www.dtic.mil/docs/citations/ADA476840">NSWC Library of Mathematical Functions</a>,
53   * available
54   * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html">here</a>.
55   * This library is "approved for public release", and the
56   * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf">Copyright guidance</a>
57   * indicates that unless otherwise stated in the code, all FORTRAN functions in
58   * this library are license free. Since no such notice appears in the code these
59   * functions can safely be ported to Hipparchus.
60   * </p>
61   *
62   */
63  public class Gamma {
64      /**
65       * <a href="http://en.wikipedia.org/wiki/Euler-Mascheroni_constant">Euler-Mascheroni constant</a>
66       */
67      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
68  
69      /**
70       * The value of the {@code g} constant in the Lanczos approximation, see
71       * {@link #lanczos(double)}.
72       */
73      public static final double LANCZOS_G = 607.0 / 128.0;
74  
75      /** Maximum allowed numerical error. */
76      private static final double DEFAULT_EPSILON = 10e-15;
77  
78      /** Lanczos coefficients */
79      private static final double[] LANCZOS = {
80          0.99999999999999709182,
81          57.156235665862923517,
82          -59.597960355475491248,
83          14.136097974741747174,
84          -0.49191381609762019978,
85          .33994649984811888699e-4,
86          .46523628927048575665e-4,
87          -.98374475304879564677e-4,
88          .15808870322491248884e-3,
89          -.21026444172410488319e-3,
90          .21743961811521264320e-3,
91          -.16431810653676389022e-3,
92          .84418223983852743293e-4,
93          -.26190838401581408670e-4,
94          .36899182659531622704e-5,
95      };
96  
97      /** Avoid repeated computation of log of 2 PI in logGamma */
98      private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(2.0 * FastMath.PI);
99  
100     /** The constant value of &radic;(2&pi;). */
101     private static final double SQRT_TWO_PI = 2.506628274631000502;
102 
103     // limits for switching algorithm in digamma
104     /** C limit. */
105     private static final double C_LIMIT = 49;
106 
107     /** S limit. */
108     private static final double S_LIMIT = 1e-8;
109 
110     /*
111      * Constants for the computation of double invGamma1pm1(double).
112      * Copied from DGAM1 in the NSWC library.
113      */
114 
115     /** The constant {@code A0} defined in {@code DGAM1}. */
116     private static final double INV_GAMMA1P_M1_A0 = .611609510448141581788E-08;
117 
118     /** The constant {@code A1} defined in {@code DGAM1}. */
119     private static final double INV_GAMMA1P_M1_A1 = .624730830116465516210E-08;
120 
121     /** The constant {@code B1} defined in {@code DGAM1}. */
122     private static final double INV_GAMMA1P_M1_B1 = .203610414066806987300E+00;
123 
124     /** The constant {@code B2} defined in {@code DGAM1}. */
125     private static final double INV_GAMMA1P_M1_B2 = .266205348428949217746E-01;
126 
127     /** The constant {@code B3} defined in {@code DGAM1}. */
128     private static final double INV_GAMMA1P_M1_B3 = .493944979382446875238E-03;
129 
130     /** The constant {@code B4} defined in {@code DGAM1}. */
131     private static final double INV_GAMMA1P_M1_B4 = -.851419432440314906588E-05;
132 
133     /** The constant {@code B5} defined in {@code DGAM1}. */
134     private static final double INV_GAMMA1P_M1_B5 = -.643045481779353022248E-05;
135 
136     /** The constant {@code B6} defined in {@code DGAM1}. */
137     private static final double INV_GAMMA1P_M1_B6 = .992641840672773722196E-06;
138 
139     /** The constant {@code B7} defined in {@code DGAM1}. */
140     private static final double INV_GAMMA1P_M1_B7 = -.607761895722825260739E-07;
141 
142     /** The constant {@code B8} defined in {@code DGAM1}. */
143     private static final double INV_GAMMA1P_M1_B8 = .195755836614639731882E-09;
144 
145     /** The constant {@code P0} defined in {@code DGAM1}. */
146     private static final double INV_GAMMA1P_M1_P0 = .6116095104481415817861E-08;
147 
148     /** The constant {@code P1} defined in {@code DGAM1}. */
149     private static final double INV_GAMMA1P_M1_P1 = .6871674113067198736152E-08;
150 
151     /** The constant {@code P2} defined in {@code DGAM1}. */
152     private static final double INV_GAMMA1P_M1_P2 = .6820161668496170657918E-09;
153 
154     /** The constant {@code P3} defined in {@code DGAM1}. */
155     private static final double INV_GAMMA1P_M1_P3 = .4686843322948848031080E-10;
156 
157     /** The constant {@code P4} defined in {@code DGAM1}. */
158     private static final double INV_GAMMA1P_M1_P4 = .1572833027710446286995E-11;
159 
160     /** The constant {@code P5} defined in {@code DGAM1}. */
161     private static final double INV_GAMMA1P_M1_P5 = -.1249441572276366213222E-12;
162 
163     /** The constant {@code P6} defined in {@code DGAM1}. */
164     private static final double INV_GAMMA1P_M1_P6 = .4343529937408594255178E-14;
165 
166     /** The constant {@code Q1} defined in {@code DGAM1}. */
167     private static final double INV_GAMMA1P_M1_Q1 = .3056961078365221025009E+00;
168 
169     /** The constant {@code Q2} defined in {@code DGAM1}. */
170     private static final double INV_GAMMA1P_M1_Q2 = .5464213086042296536016E-01;
171 
172     /** The constant {@code Q3} defined in {@code DGAM1}. */
173     private static final double INV_GAMMA1P_M1_Q3 = .4956830093825887312020E-02;
174 
175     /** The constant {@code Q4} defined in {@code DGAM1}. */
176     private static final double INV_GAMMA1P_M1_Q4 = .2692369466186361192876E-03;
177 
178     /** The constant {@code C} defined in {@code DGAM1}. */
179     private static final double INV_GAMMA1P_M1_C = -.422784335098467139393487909917598E+00;
180 
181     /** The constant {@code C0} defined in {@code DGAM1}. */
182     private static final double INV_GAMMA1P_M1_C0 = .577215664901532860606512090082402E+00;
183 
184     /** The constant {@code C1} defined in {@code DGAM1}. */
185     private static final double INV_GAMMA1P_M1_C1 = -.655878071520253881077019515145390E+00;
186 
187     /** The constant {@code C2} defined in {@code DGAM1}. */
188     private static final double INV_GAMMA1P_M1_C2 = -.420026350340952355290039348754298E-01;
189 
190     /** The constant {@code C3} defined in {@code DGAM1}. */
191     private static final double INV_GAMMA1P_M1_C3 = .166538611382291489501700795102105E+00;
192 
193     /** The constant {@code C4} defined in {@code DGAM1}. */
194     private static final double INV_GAMMA1P_M1_C4 = -.421977345555443367482083012891874E-01;
195 
196     /** The constant {@code C5} defined in {@code DGAM1}. */
197     private static final double INV_GAMMA1P_M1_C5 = -.962197152787697356211492167234820E-02;
198 
199     /** The constant {@code C6} defined in {@code DGAM1}. */
200     private static final double INV_GAMMA1P_M1_C6 = .721894324666309954239501034044657E-02;
201 
202     /** The constant {@code C7} defined in {@code DGAM1}. */
203     private static final double INV_GAMMA1P_M1_C7 = -.116516759185906511211397108401839E-02;
204 
205     /** The constant {@code C8} defined in {@code DGAM1}. */
206     private static final double INV_GAMMA1P_M1_C8 = -.215241674114950972815729963053648E-03;
207 
208     /** The constant {@code C9} defined in {@code DGAM1}. */
209     private static final double INV_GAMMA1P_M1_C9 = .128050282388116186153198626328164E-03;
210 
211     /** The constant {@code C10} defined in {@code DGAM1}. */
212     private static final double INV_GAMMA1P_M1_C10 = -.201348547807882386556893914210218E-04;
213 
214     /** The constant {@code C11} defined in {@code DGAM1}. */
215     private static final double INV_GAMMA1P_M1_C11 = -.125049348214267065734535947383309E-05;
216 
217     /** The constant {@code C12} defined in {@code DGAM1}. */
218     private static final double INV_GAMMA1P_M1_C12 = .113302723198169588237412962033074E-05;
219 
220     /** The constant {@code C13} defined in {@code DGAM1}. */
221     private static final double INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06;
222 
223     /**
224      * Default constructor.  Prohibit instantiation.
225      */
226     private Gamma() {}
227 
228     /**
229      * <p>
230      * Returns the value of log&nbsp;&Gamma;(x) for x&nbsp;&gt;&nbsp;0.
231      * </p>
232      * <p>
233      * For x &le; 8, the implementation is based on the double precision
234      * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
235      * {@code DGAMLN}. For x &gt; 8, the implementation is based on
236      * </p>
237      * <ul>
238      * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
239      *     Function</a>, equation (28).</li>
240      * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
241      *     Lanczos Approximation</a>, equations (1) through (5).</li>
242      * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
243      *     the computation of the convergent Lanczos complex Gamma
244      *     approximation</a></li>
245      * </ul>
246      *
247      * @param x Argument.
248      * @return the value of {@code log(Gamma(x))}, {@code Double.NaN} if
249      * {@code x <= 0.0}.
250      */
251     public static double logGamma(double x) {
252         double ret;
253 
254         if (Double.isNaN(x) || (x <= 0.0)) {
255             ret = Double.NaN;
256         } else if (x < 0.5) {
257             return logGamma1p(x) - FastMath.log(x);
258         } else if (x <= 2.5) {
259             return logGamma1p((x - 0.5) - 0.5);
260         } else if (x <= 8.0) {
261             final int n = (int) FastMath.floor(x - 1.5);
262             double prod = 1.0;
263             for (int i = 1; i <= n; i++) {
264                 prod *= x - i;
265             }
266             return logGamma1p(x - (n + 1)) + FastMath.log(prod);
267         } else {
268             double sum = lanczos(x);
269             double tmp = x + LANCZOS_G + .5;
270             ret = ((x + .5) * FastMath.log(tmp)) - tmp +
271                 HALF_LOG_2_PI + FastMath.log(sum / x);
272         }
273 
274         return ret;
275     }
276 
277     /**
278      * <p>
279      * Returns the value of log&nbsp;&Gamma;(x) for x&nbsp;&gt;&nbsp;0.
280      * </p>
281      * <p>
282      * For x &le; 8, the implementation is based on the double precision
283      * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
284      * {@code DGAMLN}. For x &gt; 8, the implementation is based on
285      * </p>
286      * <ul>
287      * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
288      *     Function</a>, equation (28).</li>
289      * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
290      *     Lanczos Approximation</a>, equations (1) through (5).</li>
291      * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
292      *     the computation of the convergent Lanczos complex Gamma
293      *     approximation</a></li>
294      * </ul>
295      *
296      * @param x Argument.
297      * @param <T> Type of the field elements.
298      * @return the value of {@code log(Gamma(x))}, {@code Double.NaN} if
299      * {@code x <= 0.0}.
300      */
301     public static <T extends CalculusFieldElement<T>> T logGamma(T x) {
302         final Field<T> field = x.getField();
303         T              ret;
304 
305         if (x.isNaN() || (x.getReal() <= 0.0)) {
306             ret = field.getOne().multiply(Double.NaN);
307         }
308         else if (x.getReal() < 0.5) {
309             return logGamma1p(x).subtract(x.log());
310         }
311         else if (x.getReal() <= 2.5) {
312             return logGamma1p(x.subtract(1));
313         }
314         else if (x.getReal() <= 8.0) {
315             final int n    = (int) x.subtract(1.5).floor().getReal();
316             T         prod = field.getOne();
317             for (int i = 1; i <= n; i++) {
318                 prod = prod.multiply(x.subtract(i));
319             }
320             return logGamma1p(x.subtract(n + 1)).add(prod.log());
321         }
322         else {
323             T sum = lanczos(x);
324             T tmp = x.add(LANCZOS_G + .5);
325             ret = x.add(.5).multiply(tmp.log()).subtract(tmp).add(HALF_LOG_2_PI).add(sum.divide(x).log());
326         }
327 
328         return ret;
329     }
330 
331     /**
332      * Returns the regularized gamma function P(a, x).
333      *
334      * @param a Parameter.
335      * @param x Value.
336      * @return the regularized gamma function P(a, x).
337      * @throws MathIllegalStateException if the algorithm fails to converge.
338      */
339     public static double regularizedGammaP(double a, double x) {
340         return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
341     }
342 
343     /**
344      * Returns the regularized gamma function P(a, x).
345      *
346      * @param a Parameter.
347      * @param x Value.
348      * @param <T> Type of the field elements.
349      * @return the regularized gamma function P(a, x).
350      * @throws MathIllegalStateException if the algorithm fails to converge.
351      */
352     public static <T extends CalculusFieldElement<T>> T regularizedGammaP(T a, T x) {
353         return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
354     }
355 
356     /**
357      * Returns the regularized gamma function P(a, x).
358      * <p>
359      * The implementation of this method is based on:
360      * <ul>
361      *  <li>
362      *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
363      *   Regularized Gamma Function</a>, equation (1)
364      *  </li>
365      *  <li>
366      *   <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
367      *   Incomplete Gamma Function</a>, equation (4).
368      *  </li>
369      *  <li>
370      *   <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html">
371      *   Confluent Hypergeometric Function of the First Kind</a>, equation (1).
372      *  </li>
373      * </ul>
374      *
375      * @param a the a parameter.
376      * @param x the value.
377      * @param epsilon When the absolute value of the nth item in the
378      * series is less than epsilon the approximation ceases to calculate
379      * further elements in the series.
380      * @param maxIterations Maximum number of "iterations" to complete.
381      * @return the regularized gamma function P(a, x)
382      * @throws MathIllegalStateException if the algorithm fails to converge.
383      */
384     public static double regularizedGammaP(double a,
385                                            double x,
386                                            double epsilon,
387                                            int maxIterations) {
388         double ret;
389 
390         if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
391             ret = Double.NaN;
392         } else if (x == 0.0) {
393             ret = 0.0;
394         } else if (x >= a + 1) {
395             // use regularizedGammaQ because it should converge faster in this
396             // case.
397             ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
398         } else {
399             // calculate series
400             double n = 0.0; // current element index
401             double an = 1.0 / a; // n-th element in the series
402             double sum = an; // partial sum
403             while (FastMath.abs(an/sum) > epsilon &&
404                    n < maxIterations &&
405                    sum < Double.POSITIVE_INFINITY) {
406                 // compute next element in the series
407                 n += 1.0;
408                 an *= x / (a + n);
409 
410                 // update partial sum
411                 sum += an;
412             }
413             if (n >= maxIterations) {
414                 throw new MathIllegalStateException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, maxIterations);
415             } else if (Double.isInfinite(sum)) {
416                 ret = 1.0;
417             } else {
418                 ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * sum;
419             }
420         }
421 
422         return ret;
423     }
424 
425     /**
426      * Returns the regularized gamma function P(a, x).
427      * <p>
428      * The implementation of this method is based on:
429      * <ul>
430      *  <li>
431      *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
432      *   Regularized Gamma Function</a>, equation (1)
433      *  </li>
434      *  <li>
435      *   <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
436      *   Incomplete Gamma Function</a>, equation (4).
437      *  </li>
438      *  <li>
439      *   <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html">
440      *   Confluent Hypergeometric Function of the First Kind</a>, equation (1).
441      *  </li>
442      * </ul>
443      *
444      * @param a the a parameter.
445      * @param x the value.
446      * @param epsilon When the absolute value of the nth item in the
447      * series is less than epsilon the approximation ceases to calculate
448      * further elements in the series.
449      * @param maxIterations Maximum number of "iterations" to complete.
450      * @param <T> Type of the field elements.
451      * @return the regularized gamma function P(a, x)
452      * @throws MathIllegalStateException if the algorithm fails to converge.
453      */
454     public static <T extends CalculusFieldElement<T>> T regularizedGammaP(T a,
455                                            T x,
456                                            double epsilon,
457                                            int maxIterations) {
458         final Field<T> field = x.getField();
459         final T        zero  = field.getZero();
460         final T        one   = field.getOne();
461 
462         T ret;
463 
464         if (a.isNaN() || x.isNaN() || (a.getReal() <= 0.0) || (x.getReal() < 0.0)) {
465             ret = one.multiply(Double.NaN);
466         }
467         else if (x.getReal() == 0.0) {
468             ret = zero;
469         }
470         else if (x.getReal() >= a.add(1).getReal()) {
471             // use regularizedGammaQ because it should converge faster in this
472             // case.
473             ret = one.subtract(regularizedGammaQ(a, x, epsilon, maxIterations));
474         }
475         else {
476             // calculate series
477             double n   = 0.0; // current element index
478             T      an  = one.divide(a); // n-th element in the series
479             T      sum = an; // partial sum
480             while (an.divide(sum).abs().getReal() > epsilon &&
481                     n < maxIterations &&
482                     sum.getReal() < Double.POSITIVE_INFINITY) {
483                 // compute next element in the series
484                 n += 1.0;
485                 an = an.multiply(x.divide(a.add(n)));
486 
487                 // update partial sum
488                 sum = sum.add(an);
489             }
490             if (n >= maxIterations) {
491                 throw new MathIllegalStateException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, maxIterations);
492             }
493             else if (sum.isInfinite()) {
494                 ret = one;
495             }
496             else {
497                 ret = a.multiply(x.log()).subtract(logGamma(a)).subtract(x).exp().multiply(sum);
498             }
499         }
500 
501         return ret;
502     }
503 
504     /**
505      * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
506      *
507      * @param a the a parameter.
508      * @param x the value.
509      * @return the regularized gamma function Q(a, x)
510      * @throws MathIllegalStateException if the algorithm fails to converge.
511      */
512     public static double regularizedGammaQ(double a, double x) {
513         return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
514     }
515 
516     /**
517      * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
518      *
519      * @param a the a parameter.
520      * @param x the value.
521      * @param <T> Type of the field elements.
522      * @return the regularized gamma function Q(a, x)
523      * @throws MathIllegalStateException if the algorithm fails to converge.
524      */
525     public static <T extends CalculusFieldElement<T>> T regularizedGammaQ(T a, T x) {
526         return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
527     }
528 
529     /**
530      * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
531      * <p>
532      * The implementation of this method is based on:
533      * <ul>
534      *  <li>
535      *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
536      *   Regularized Gamma Function</a>, equation (1).
537      *  </li>
538      *  <li>
539      *   <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
540      *   Regularized incomplete gamma function: Continued fraction representations
541      *   (formula 06.08.10.0003)</a>
542      *  </li>
543      * </ul>
544      *
545      * @param a the a parameter.
546      * @param x the value.
547      * @param epsilon When the absolute value of the nth item in the
548      * series is less than epsilon the approximation ceases to calculate
549      * further elements in the series.
550      * @param maxIterations Maximum number of "iterations" to complete.
551      * @return the regularized gamma function P(a, x)
552      * @throws MathIllegalStateException if the algorithm fails to converge.
553      */
554     public static double regularizedGammaQ(final double a,
555                                            double x,
556                                            double epsilon,
557                                            int maxIterations) {
558         double ret;
559 
560         if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
561             ret = Double.NaN;
562         } else if (x == 0.0) {
563             ret = 1.0;
564         } else if (x < a + 1.0) {
565             // use regularizedGammaP because it should converge faster in this
566             // case.
567             ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations);
568         } else {
569             // create continued fraction
570             ContinuedFraction cf = new ContinuedFraction() {
571 
572                 /** {@inheritDoc} */
573                 @Override
574                 protected double getA(int n, double x) {
575                     return ((2.0 * n) + 1.0) - a + x;
576                 }
577 
578                 /** {@inheritDoc} */
579                 @Override
580                 protected double getB(int n, double x) {
581                     return n * (a - n);
582                 }
583             };
584 
585             ret = 1.0 / cf.evaluate(x, epsilon, maxIterations);
586             ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * ret;
587         }
588 
589         return ret;
590     }
591 
592     /**
593      * Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
594      * <p>
595      * The implementation of this method is based on:
596      * <ul>
597      *  <li>
598      *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
599      *   Regularized Gamma Function</a>, equation (1).
600      *  </li>
601      *  <li>
602      *   <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
603      *   Regularized incomplete gamma function: Continued fraction representations
604      *   (formula 06.08.10.0003)</a>
605      *  </li>
606      * </ul>
607      *
608      * @param a the a parameter.
609      * @param x the value.
610      * @param epsilon When the absolute value of the nth item in the
611      * series is less than epsilon the approximation ceases to calculate
612      * further elements in the series.
613      * @param maxIterations Maximum number of "iterations" to complete.
614      * @param <T> Type fo the field elements.
615      * @return the regularized gamma function P(a, x)
616      * @throws MathIllegalStateException if the algorithm fails to converge.
617      */
618     public static <T extends CalculusFieldElement<T>> T regularizedGammaQ(final T a,
619                                                                           T x,
620                                                                           double epsilon,
621                                                                           int maxIterations) {
622         final Field<T> field = x.getField();
623         final T        one   = field.getOne();
624 
625         T ret;
626 
627         if (a.isNaN() || x.isNaN() || a.getReal() <= 0.0 || x.getReal() < 0.0) {
628             ret = field.getOne().multiply(Double.NaN);
629         }
630         else if (x.getReal() == 0.0) {
631             ret = one;
632         }
633         else if (x.getReal() < a.add(1.0).getReal()) {
634             // use regularizedGammaP because it should converge faster in this
635             // case.
636             ret = one.subtract(regularizedGammaP(a, x, epsilon, maxIterations));
637         }
638         else {
639             // create continued fraction
640             FieldContinuedFraction cf = new FieldContinuedFraction() {
641 
642                 /** {@inheritDoc} */
643                 @Override
644                 @SuppressWarnings("unchecked")
645                 public <C extends CalculusFieldElement<C>> C getA(final int n, final C x) {
646                     return x.subtract((C) a).add((2.0 * n) + 1.0);
647                 }
648 
649                 /** {@inheritDoc} */
650                 @Override
651                 @SuppressWarnings("unchecked")
652                 public <C extends CalculusFieldElement<C>> C getB(final int n, final C x) {
653                     return (C) a.subtract(n).multiply(n);
654                 }
655             };
656 
657             ret = one.divide(cf.evaluate(x, epsilon, maxIterations));
658             ret = a.multiply(x.log()).subtract(logGamma(a)).subtract(x).exp().multiply(ret);
659         }
660 
661         return ret;
662     }
663 
664     /**
665      * <p>Computes the digamma function of x.</p>
666      *
667      * <p>This is an independently written implementation of the algorithm described in
668      * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p>
669      *
670      * <p>Some of the constants have been changed to increase accuracy at the moderate expense
671      * of run-time.  The result should be accurate to within 10^-8 absolute tolerance for
672      * x &gt;= 10^-5 and within 10^-8 relative tolerance for x &gt; 0.</p>
673      *
674      * <p>Performance for large negative values of x will be quite expensive (proportional to
675      * |x|).  Accuracy for negative values of x should be about 10^-8 absolute for results
676      * less than 10^5 and 10^-8 relative for results larger than that.</p>
677      *
678      * @param x Argument.
679      * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller.
680      * @see <a href="http://en.wikipedia.org/wiki/Digamma_function">Digamma</a>
681      * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf">Bernardo&#39;s original article </a>
682      */
683     public static double digamma(double x) {
684         if (Double.isNaN(x) || Double.isInfinite(x)) {
685             return x;
686         }
687 
688         if (x > 0 && x <= S_LIMIT) {
689             // use method 5 from Bernardo AS103
690             // accurate to O(x)
691             return -GAMMA - 1 / x;
692         }
693         if (x >= C_LIMIT) {
694             // use method 8 (accurate to O(1/x^8))
695             double inv = 1 / (x * x);
696             //            1       1        1         1         1         5           691         1
697             // log(x) -  --- - ------ + ------- - ------- + ------- - ------- +  ---------- - -------
698             //           2 x   12 x^2   120 x^4   252 x^6   240 x^8   660 x^10   32760 x^12   12 x^14
699             return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv * (1.0 / 252 + inv *
700                     (1.0 / 240 - inv * (5.0 / 660 + inv * (691.0 / 32760 - inv / 12))))));
701         }
702 
703         return digamma(x + 1) - 1 / x;
704     }
705 
706     /**
707      * <p>Computes the digamma function of x.</p>
708      *
709      * <p>This is an independently written implementation of the algorithm described in
710      * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.</p>
711      *
712      * <p>Some of the constants have been changed to increase accuracy at the moderate expense
713      * of run-time.  The result should be accurate to within 10^-8 absolute tolerance for
714      * x &gt;= 10^-5 and within 10^-8 relative tolerance for x &gt; 0.</p>
715      *
716      * <p>Performance for large negative values of x will be quite expensive (proportional to
717      * |x|).  Accuracy for negative values of x should be about 10^-8 absolute for results
718      * less than 10^5 and 10^-8 relative for results larger than that.</p>
719      *
720      * @param x Argument.
721      * @param <T> Type of the field elements.
722      * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller.
723      * @see <a href="http://en.wikipedia.org/wiki/Digamma_function">Digamma</a>
724      * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf">Bernardo&#39;s original article </a>
725      */
726     public static <T extends CalculusFieldElement<T>> T digamma(T x) {
727         if (x.isNaN() || x.isInfinite()) {
728             return x;
729         }
730 
731         if (x.getReal() > 0 && x.getReal() <= S_LIMIT) {
732             // use method 5 from Bernardo AS103
733             // accurate to O(x)
734             return x.pow(-1).negate().subtract(GAMMA);
735         }
736 
737         if (x.getReal() >= C_LIMIT) {
738             // use method 8 (accurate to O(1/x^8))
739             T inv = x.multiply(x).reciprocal();
740             //            1       1        1         1         1         5           691         1
741             // log(x) -  --- - ------ + ------- - ------- + ------- - ------- +  ---------- - -------
742             //           2 x   12 x^2   120 x^4   252 x^6   240 x^8   660 x^10   32760 x^12   12 x^14
743             return x.log().subtract(x.pow(-1).multiply(0.5)).add(
744                     inv.multiply(
745                             inv.multiply(
746                                     inv.multiply(
747                                             inv.multiply(
748                                                     inv.multiply(
749                                                             inv.multiply(inv.divide(-12.)
750                                                                             .add(691. / 32760))
751                                                                .subtract(5. / 660))
752                                                        .add(1.0 / 240))
753                                                .subtract(1.0 / 252))
754                                        .add(1.0 / 120))
755                                .subtract(1.0 / 12)));
756         }
757 
758         return digamma(x.add(1.)).subtract(x.pow(-1));
759     }
760 
761     /**
762      * Computes the trigamma function of x.
763      * This function is derived by taking the derivative of the implementation
764      * of digamma.
765      *
766      * @param x Argument.
767      * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller
768      * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma</a>
769      * @see Gamma#digamma(double)
770      */
771     public static double trigamma(double x) {
772         if (Double.isNaN(x) || Double.isInfinite(x)) {
773             return x;
774         }
775 
776         if (x > 0 && x <= S_LIMIT) {
777             return 1 / (x * x);
778         }
779 
780         if (x >= C_LIMIT) {
781             double inv = 1 / (x * x);
782             //  1    1      1       1       1      1         5        691        7
783             //  - + ---- + ---- - ----- + ----- - ----- + ------- - -------- + ------
784             //  x      2      3       5       7       9        11         13       15
785             //      2 x    6 x    30 x    42 x    30 x    66 x      2730 x      6 x
786             return 1 / x + inv * 0.5 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv * (1.0 / 42 - inv * (1.0 / 30 + inv *
787                     (5.0 / 66 - inv * (691. / 2730 + inv * 7. / 15))))));
788         }
789 
790         return trigamma(x + 1) + 1 / (x * x);
791     }
792 
793     /**
794      * Computes the trigamma function of x.
795      * This function is derived by taking the derivative of the implementation
796      * of digamma.
797      *
798      * @param x Argument.
799      * @param <T> Type of the field elements.
800      * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller
801      * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma</a>
802      * @see Gamma#digamma(double)
803      */
804     public static <T extends CalculusFieldElement<T>> T trigamma(T x) {
805         if (x.isNaN() || x.isInfinite()) {
806             return x;
807         }
808 
809         if (x.getReal() > 0 && x.getReal() <= S_LIMIT) {
810             // use method 5 from Bernardo AS103
811             // accurate to O(x)
812             return x.multiply(x).reciprocal();
813         }
814 
815         if (x.getReal() >= C_LIMIT) {
816             // use method 4 (accurate to O(1/x^8)
817             T inv    = x.multiply(x).reciprocal();
818             T invCub = inv.multiply(x.reciprocal());
819             //  1    1      1       1       1      1         5        691        7
820             //  - + ---- + ---- - ----- + ----- + ----- + ------- - -------- + ------
821             //  x      2      3       5       7       9        11         13       15
822             //      2 x    6 x    30 x    42 x    30 x    66 x      2730 x      6 x
823             return x.pow(-1).add(
824                     inv.multiply(0.5)).add(
825                             invCub.multiply(
826                                     inv.multiply(
827                                             inv.multiply(
828                                                     inv.multiply(
829                                                             inv.multiply(
830                                                                     inv.multiply(inv.multiply(7. / 6)
831                                                                                     .subtract(691. / 2730))
832                                                                        .add(5. / 66))
833                                                                .subtract(1.0 / 30))
834                                                        .add(1.0 / 42))
835                                                .subtract(1.0 / 30))
836                                        .add(1.0 / 6)));
837         }
838 
839         return trigamma(x.add(1.)).add(x.multiply(x).reciprocal());
840     }
841 
842     /**
843      * <p>
844      * Returns the Lanczos approximation used to compute the gamma function.
845      * The Lanczos approximation is related to the Gamma function by the
846      * following equation
847      * \[
848      * \Gamma(x) = \frac{\sqrt{2\pi}}{x} \times (x + g + \frac{1}{2}) ^ (x + \frac{1}{2})
849      *                   \times e^{-x - g - 0.5} \times \mathrm{lanczos}(x)
850      * \]
851      * where {@code g} is the Lanczos constant.
852      * </p>
853      *
854      * @param x Argument.
855      * @return The Lanczos approximation.
856      * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos Approximation</a>
857      * equations (1) through (5), and Paul Godfrey's
858      * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation
859      * of the convergent Lanczos complex Gamma approximation</a>
860      */
861     public static double lanczos(final double x) {
862         double sum = 0.0;
863         for (int i = LANCZOS.length - 1; i > 0; --i) {
864             sum += LANCZOS[i] / (x + i);
865         }
866         return sum + LANCZOS[0];
867     }
868 
869     /**
870      * <p>
871      * Returns the Lanczos approximation used to compute the gamma function.
872      * The Lanczos approximation is related to the Gamma function by the
873      * following equation
874      * \[
875      * \Gamma(x) = \frac{\sqrt{2\pi}}{x} \times (x + g + \frac{1}{2}) ^ (x + \frac{1}{2})
876      *                   \times e^{-x - g - 0.5} \times \mathrm{lanczos}(x)
877      * \]
878      * where {@code g} is the Lanczos constant.
879      * </p>
880      *
881      * @param x Argument.
882      * @param <T> Type of the field elements.
883      * @return The Lanczos approximation.
884      * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos Approximation</a>
885      * equations (1) through (5), and Paul Godfrey's
886      * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation
887      * of the convergent Lanczos complex Gamma approximation</a>
888      */
889     public static <T extends CalculusFieldElement<T>> T lanczos(final T x) {
890         final Field<T> field = x.getField();
891         T              sum   = field.getZero();
892         for (int i = LANCZOS.length - 1; i > 0; --i) {
893             sum = sum.add(x.add(i).pow(-1.).multiply(LANCZOS[i]));
894         }
895         return sum.add(LANCZOS[0]);
896     }
897 
898     /**
899      * Returns the value of 1 / &Gamma;(1 + x) - 1 for -0&#46;5 &le; x &le;
900      * 1&#46;5. This implementation is based on the double precision
901      * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
902      * {@code DGAM1}.
903      *
904      * @param x Argument.
905      * @return The value of {@code 1.0 / Gamma(1.0 + x) - 1.0}.
906      * @throws MathIllegalArgumentException if {@code x < -0.5}
907      * @throws MathIllegalArgumentException if {@code x > 1.5}
908      */
909     public static double invGamma1pm1(final double x) {
910 
911         if (x < -0.5) {
912             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
913                                                    x, -0.5);
914         }
915         if (x > 1.5) {
916             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
917                                                    x, 1.5);
918         }
919 
920         final double ret;
921         final double t = x <= 0.5 ? x : (x - 0.5) - 0.5;
922         if (t < 0.0) {
923             final double a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1;
924             double b = INV_GAMMA1P_M1_B8;
925             b = INV_GAMMA1P_M1_B7 + t * b;
926             b = INV_GAMMA1P_M1_B6 + t * b;
927             b = INV_GAMMA1P_M1_B5 + t * b;
928             b = INV_GAMMA1P_M1_B4 + t * b;
929             b = INV_GAMMA1P_M1_B3 + t * b;
930             b = INV_GAMMA1P_M1_B2 + t * b;
931             b = INV_GAMMA1P_M1_B1 + t * b;
932             b = 1.0 + t * b;
933 
934             double c = INV_GAMMA1P_M1_C13 + t * (a / b);
935             c = INV_GAMMA1P_M1_C12 + t * c;
936             c = INV_GAMMA1P_M1_C11 + t * c;
937             c = INV_GAMMA1P_M1_C10 + t * c;
938             c = INV_GAMMA1P_M1_C9 + t * c;
939             c = INV_GAMMA1P_M1_C8 + t * c;
940             c = INV_GAMMA1P_M1_C7 + t * c;
941             c = INV_GAMMA1P_M1_C6 + t * c;
942             c = INV_GAMMA1P_M1_C5 + t * c;
943             c = INV_GAMMA1P_M1_C4 + t * c;
944             c = INV_GAMMA1P_M1_C3 + t * c;
945             c = INV_GAMMA1P_M1_C2 + t * c;
946             c = INV_GAMMA1P_M1_C1 + t * c;
947             c = INV_GAMMA1P_M1_C + t * c;
948             if (x > 0.5) {
949                 ret = t * c / x;
950             } else {
951                 ret = x * ((c + 0.5) + 0.5);
952             }
953         } else {
954             double p = INV_GAMMA1P_M1_P6;
955             p = INV_GAMMA1P_M1_P5 + t * p;
956             p = INV_GAMMA1P_M1_P4 + t * p;
957             p = INV_GAMMA1P_M1_P3 + t * p;
958             p = INV_GAMMA1P_M1_P2 + t * p;
959             p = INV_GAMMA1P_M1_P1 + t * p;
960             p = INV_GAMMA1P_M1_P0 + t * p;
961 
962             double q = INV_GAMMA1P_M1_Q4;
963             q = INV_GAMMA1P_M1_Q3 + t * q;
964             q = INV_GAMMA1P_M1_Q2 + t * q;
965             q = INV_GAMMA1P_M1_Q1 + t * q;
966             q = 1.0 + t * q;
967 
968             double c = INV_GAMMA1P_M1_C13 + (p / q) * t;
969             c = INV_GAMMA1P_M1_C12 + t * c;
970             c = INV_GAMMA1P_M1_C11 + t * c;
971             c = INV_GAMMA1P_M1_C10 + t * c;
972             c = INV_GAMMA1P_M1_C9 + t * c;
973             c = INV_GAMMA1P_M1_C8 + t * c;
974             c = INV_GAMMA1P_M1_C7 + t * c;
975             c = INV_GAMMA1P_M1_C6 + t * c;
976             c = INV_GAMMA1P_M1_C5 + t * c;
977             c = INV_GAMMA1P_M1_C4 + t * c;
978             c = INV_GAMMA1P_M1_C3 + t * c;
979             c = INV_GAMMA1P_M1_C2 + t * c;
980             c = INV_GAMMA1P_M1_C1 + t * c;
981             c = INV_GAMMA1P_M1_C0 + t * c;
982 
983             if (x > 0.5) {
984                 ret = (t / x) * ((c - 0.5) - 0.5);
985             } else {
986                 ret = x * c;
987             }
988         }
989 
990         return ret;
991     }
992 
993     /**
994      * Returns the value of 1 / &Gamma;(1 + x) - 1 for -0&#46;5 &le; x &le;
995      * 1&#46;5. This implementation is based on the double precision
996      * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
997      * {@code DGAM1}.
998      *
999      * @param x Argument.
1000      * @param <T> Type of the field elements.
1001      * @return The value of {@code 1.0 / Gamma(1.0 + x) - 1.0}.
1002      * @throws MathIllegalArgumentException if {@code x < -0.5}
1003      * @throws MathIllegalArgumentException if {@code x > 1.5}
1004      */
1005     public static <T extends CalculusFieldElement<T>> T invGamma1pm1(final T x) {
1006         final T one = x.getField().getOne();
1007 
1008         if (x.getReal() < -0.5) {
1009             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
1010                                                    x, -0.5);
1011         }
1012         if (x.getReal() > 1.5) {
1013             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
1014                                                    x, 1.5);
1015         }
1016 
1017         final T ret;
1018         final T t = x.getReal() <= 0.5 ? x : x.subtract(1);
1019         if (t.getReal() < 0.0) {
1020             final T a = one.multiply(INV_GAMMA1P_M1_A0).add(t.multiply(INV_GAMMA1P_M1_A1));
1021             T       b = one.multiply(INV_GAMMA1P_M1_B8);
1022             b = one.multiply(INV_GAMMA1P_M1_B7).add(t.multiply(b));
1023             b = one.multiply(INV_GAMMA1P_M1_B6).add(t.multiply(b));
1024             b = one.multiply(INV_GAMMA1P_M1_B5).add(t.multiply(b));
1025             b = one.multiply(INV_GAMMA1P_M1_B4).add(t.multiply(b));
1026             b = one.multiply(INV_GAMMA1P_M1_B3).add(t.multiply(b));
1027             b = one.multiply(INV_GAMMA1P_M1_B2).add(t.multiply(b));
1028             b = one.multiply(INV_GAMMA1P_M1_B1).add(t.multiply(b));
1029             b = one.add(t.multiply(b));
1030 
1031             T c = one.multiply(INV_GAMMA1P_M1_C13).add(t.multiply(a.divide(b)));
1032             c = one.multiply(INV_GAMMA1P_M1_C12).add(t.multiply(c));
1033             c = one.multiply(INV_GAMMA1P_M1_C11).add(t.multiply(c));
1034             c = one.multiply(INV_GAMMA1P_M1_C10).add(t.multiply(c));
1035             c = one.multiply(INV_GAMMA1P_M1_C9).add(t.multiply(c));
1036             c = one.multiply(INV_GAMMA1P_M1_C8).add(t.multiply(c));
1037             c = one.multiply(INV_GAMMA1P_M1_C7).add(t.multiply(c));
1038             c = one.multiply(INV_GAMMA1P_M1_C6).add(t.multiply(c));
1039             c = one.multiply(INV_GAMMA1P_M1_C5).add(t.multiply(c));
1040             c = one.multiply(INV_GAMMA1P_M1_C4).add(t.multiply(c));
1041             c = one.multiply(INV_GAMMA1P_M1_C3).add(t.multiply(c));
1042             c = one.multiply(INV_GAMMA1P_M1_C2).add(t.multiply(c));
1043             c = one.multiply(INV_GAMMA1P_M1_C1).add(t.multiply(c));
1044             c = one.multiply(INV_GAMMA1P_M1_C).add(t.multiply(c));
1045             if (x.getReal() > 0.5) {
1046                 ret = t.multiply(c).divide(x);
1047             }
1048             else {
1049                 ret = x.multiply(c.add(1));
1050             }
1051         }
1052         else {
1053             T p = one.multiply(INV_GAMMA1P_M1_P6);
1054             p = one.multiply(INV_GAMMA1P_M1_P5).add(t.multiply(p));
1055             p = one.multiply(INV_GAMMA1P_M1_P4).add(t.multiply(p));
1056             p = one.multiply(INV_GAMMA1P_M1_P3).add(t.multiply(p));
1057             p = one.multiply(INV_GAMMA1P_M1_P2).add(t.multiply(p));
1058             p = one.multiply(INV_GAMMA1P_M1_P1).add(t.multiply(p));
1059             p = one.multiply(INV_GAMMA1P_M1_P0).add(t.multiply(p));
1060 
1061             T q = one.multiply(INV_GAMMA1P_M1_Q4);
1062             q = one.multiply(INV_GAMMA1P_M1_Q3).add(t.multiply(q));
1063             q = one.multiply(INV_GAMMA1P_M1_Q2).add(t.multiply(q));
1064             q = one.multiply(INV_GAMMA1P_M1_Q1).add(t.multiply(q));
1065             q = one.add(t.multiply(q));
1066 
1067             T c = one.multiply(INV_GAMMA1P_M1_C13).add(t.multiply(p.divide(q)));
1068             c = one.multiply(INV_GAMMA1P_M1_C12).add(t.multiply(c));
1069             c = one.multiply(INV_GAMMA1P_M1_C11).add(t.multiply(c));
1070             c = one.multiply(INV_GAMMA1P_M1_C10).add(t.multiply(c));
1071             c = one.multiply(INV_GAMMA1P_M1_C9).add(t.multiply(c));
1072             c = one.multiply(INV_GAMMA1P_M1_C8).add(t.multiply(c));
1073             c = one.multiply(INV_GAMMA1P_M1_C7).add(t.multiply(c));
1074             c = one.multiply(INV_GAMMA1P_M1_C6).add(t.multiply(c));
1075             c = one.multiply(INV_GAMMA1P_M1_C5).add(t.multiply(c));
1076             c = one.multiply(INV_GAMMA1P_M1_C4).add(t.multiply(c));
1077             c = one.multiply(INV_GAMMA1P_M1_C3).add(t.multiply(c));
1078             c = one.multiply(INV_GAMMA1P_M1_C2).add(t.multiply(c));
1079             c = one.multiply(INV_GAMMA1P_M1_C1).add(t.multiply(c));
1080             c = one.multiply(INV_GAMMA1P_M1_C0).add(t.multiply(c));
1081 
1082             if (x.getReal() > 0.5) {
1083                 ret = t.divide(x).multiply(c.subtract(1));
1084             }
1085             else {
1086                 ret = x.multiply(c);
1087             }
1088         }
1089 
1090         return ret;
1091     }
1092 
1093     /**
1094      * Returns the value of log &Gamma;(1 + x) for -0&#46;5 &le; x &le; 1&#46;5.
1095      * This implementation is based on the double precision implementation in
1096      * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}.
1097      *
1098      * @param x Argument.
1099      * @return The value of {@code log(Gamma(1 + x))}.
1100      * @throws MathIllegalArgumentException if {@code x < -0.5}.
1101      * @throws MathIllegalArgumentException if {@code x > 1.5}.
1102      */
1103     public static double logGamma1p(final double x)
1104         throws MathIllegalArgumentException {
1105 
1106         if (x < -0.5) {
1107             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
1108                                                    x, -0.5);
1109         }
1110         if (x > 1.5) {
1111             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
1112                                                    x, 1.5);
1113         }
1114 
1115         return -FastMath.log1p(invGamma1pm1(x));
1116     }
1117 
1118     /**
1119      * Returns the value of log &Gamma;(1 + x) for -0&#46;5 &le; x &le; 1&#46;5.
1120      * This implementation is based on the double precision implementation in
1121      * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}.
1122      *
1123      * @param x Argument.
1124      * @param <T> Type of the field elements.
1125      * @return The value of {@code log(Gamma(1 + x))}.
1126      * @throws MathIllegalArgumentException if {@code x < -0.5}.
1127      * @throws MathIllegalArgumentException if {@code x > 1.5}.
1128      */
1129     public static <T extends CalculusFieldElement<T>> T logGamma1p(final T x)
1130             throws MathIllegalArgumentException {
1131 
1132         if (x.getReal() < -0.5) {
1133             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
1134                                                    x, -0.5);
1135         }
1136         if (x.getReal() > 1.5) {
1137             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
1138                                                    x, 1.5);
1139         }
1140 
1141         return invGamma1pm1(x).log1p().negate();
1142     }
1143 
1144 
1145     /**
1146      * Returns the value of Γ(x). Based on the <em>NSWC Library of
1147      * Mathematics Subroutines</em> double precision implementation,
1148      * {@code DGAMMA}.
1149      *
1150      * @param x Argument.
1151      * @return the value of {@code Gamma(x)}.
1152      */
1153     public static double gamma(final double x) {
1154 
1155         if ((x == FastMath.rint(x)) && (x <= 0.0)) {
1156             return Double.NaN;
1157         }
1158 
1159         final double ret;
1160         final double absX = FastMath.abs(x);
1161         if (absX <= 20.0) {
1162             if (x >= 1.0) {
1163                 /*
1164                  * From the recurrence relation
1165                  * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
1166                  * then
1167                  * Gamma(t) = 1 / [1 + invGamma1pm1(t - 1)],
1168                  * where t = x - n. This means that t must satisfy
1169                  * -0.5 <= t - 1 <= 1.5.
1170                  */
1171                 double prod = 1.0;
1172                 double t = x;
1173                 while (t > 2.5) {
1174                     t -= 1.0;
1175                     prod *= t;
1176                 }
1177                 ret = prod / (1.0 + invGamma1pm1(t - 1.0));
1178             } else {
1179                 /*
1180                  * From the recurrence relation
1181                  * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
1182                  * then
1183                  * Gamma(x + n + 1) = 1 / [1 + invGamma1pm1(x + n)],
1184                  * which requires -0.5 <= x + n <= 1.5.
1185                  */
1186                 double prod = x;
1187                 double t = x;
1188                 while (t < -0.5) {
1189                     t += 1.0;
1190                     prod *= t;
1191                 }
1192                 ret = 1.0 / (prod * (1.0 + invGamma1pm1(t)));
1193             }
1194         } else {
1195             final double y = absX + LANCZOS_G + 0.5;
1196             final double gammaAbs = SQRT_TWO_PI / absX *
1197                                     FastMath.pow(y, absX + 0.5) *
1198                                     FastMath.exp(-y) * lanczos(absX);
1199             if (x > 0.0) {
1200                 ret = gammaAbs;
1201             } else {
1202                 /*
1203                  * From the reflection formula
1204                  * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
1205                  * and the recurrence relation
1206                  * Gamma(1 - x) = -x * Gamma(-x),
1207                  * it is found
1208                  * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
1209                  */
1210                 ret = -FastMath.PI /
1211                       (x * FastMath.sin(FastMath.PI * x) * gammaAbs);
1212             }
1213         }
1214         return ret;
1215     }
1216 
1217     /**
1218      * Returns the value of Γ(x). Based on the <em>NSWC Library of
1219      * Mathematics Subroutines</em> double precision implementation,
1220      * {@code DGAMMA}.
1221      *
1222      * @param x Argument.
1223      * @param <T> Type of the field elements.
1224      * @return the value of {@code Gamma(x)}.
1225      */
1226     public static <T extends CalculusFieldElement<T>> T gamma(final T x) {
1227         final T one = x.getField().getOne();
1228 
1229         if ((x.getReal() == x.rint().getReal()) && (x.getReal() <= 0.0)) {
1230             return one.multiply(Double.NaN);
1231         }
1232 
1233         final T ret;
1234         final T absX = x.abs();
1235         if (absX.getReal() <= 20.0) {
1236             if (x.getReal() >= 1.0) {
1237                 /*
1238                  * From the recurrence relation
1239                  * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
1240                  * then
1241                  * Gamma(t) = 1 / [1 + invGamma1pm1(t - 1)],
1242                  * where t = x - n. This means that t must satisfy
1243                  * -0.5 <= t - 1 <= 1.5.
1244                  */
1245                 T prod = one;
1246                 T t    = x;
1247                 while (t.getReal() > 2.5) {
1248                     t    = t.subtract(1.0);
1249                     prod = prod.multiply(t);
1250                 }
1251                 ret = prod.divide(invGamma1pm1(t.subtract(1.0)).add(1.0));
1252             }
1253             else {
1254                 /*
1255                  * From the recurrence relation
1256                  * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
1257                  * then
1258                  * Gamma(x + n + 1) = 1 / [1 + invGamma1pm1(x + n)],
1259                  * which requires -0.5 <= x + n <= 1.5.
1260                  */
1261                 T prod = x;
1262                 T t    = x;
1263                 while (t.getReal() < -0.5) {
1264                     t    = t.add(1.0);
1265                     prod = prod.multiply(t);
1266                 }
1267                 ret = prod.multiply(invGamma1pm1(t).add(1)).reciprocal();
1268             }
1269         }
1270         else {
1271             final T y = absX.add(LANCZOS_G + 0.5);
1272             final T gammaAbs = absX.reciprocal().multiply(SQRT_TWO_PI).multiply(y.pow(absX.add(0.5)))
1273                                    .multiply(y.negate().exp()).multiply(lanczos(absX));
1274             if (x.getReal() > 0.0) {
1275                 ret = gammaAbs;
1276             }
1277             else {
1278                 /*
1279                  * From the reflection formula
1280                  * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
1281                  * and the recurrence relation
1282                  * Gamma(1 - x) = -x * Gamma(-x),
1283                  * it is found
1284                  * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
1285                  */
1286                 ret = x.multiply(x.multiply(FastMath.PI).sin()).multiply(gammaAbs).reciprocal().multiply(-FastMath.PI);
1287             }
1288         }
1289         return ret;
1290     }
1291 }