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    *      http://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  package org.apache.commons.rng.sampling.distribution;
18  
19  /**
20   * <h3>
21   *  Adapted and stripped down copy of class
22   *  {@code "org.apache.commons.math4.special.Gamma"}.
23   *  TODO: Include it in a "core" component upon which high-level functionality
24   *  such as sampling can depend.
25   * </h3>
26   *
27   * <p>
28   * This is a utility class that provides computation methods related to the
29   * &Gamma; (Gamma) family of functions.
30   * </p>
31   * <p>
32   * Implementation of {@link #invGamma1pm1(double)} and
33   * {@link #logGamma1p(double)} is based on the algorithms described in
34   * <ul>
35   * <li><a href="http://dx.doi.org/10.1145/22721.23109">Didonato and Morris
36   * (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios and
37   *     their Inverse</em>, TOMS 12(4), 377-393,</li>
38   * <li><a href="http://dx.doi.org/10.1145/131766.131776">Didonato and Morris
39   * (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the
40   *     Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li>
41   * </ul>
42   * and implemented in the
43   * <a href="http://www.dtic.mil/docs/citations/ADA476840">NSWC Library of Mathematical Functions</a>,
44   * available
45   * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html">here</a>.
46   * This library is "approved for public release", and the
47   * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf">Copyright guidance</a>
48   * indicates that unless otherwise stated in the code, all FORTRAN functions in
49   * this library are license free. Since no such notice appears in the code these
50   * functions can safely be ported to Commons-Math.
51   * </p>
52   */
53  class InternalGamma { // Class is package-private on purpose; do not make it public.
54      /**
55       * Constant \( g = \frac{607}{128} \) in the {@link #lanczos(double) Lanczos approximation}.
56       */
57      public static final double LANCZOS_G = 607.0 / 128.0;
58  
59      /** Lanczos coefficients */
60      private static final double[] LANCZOS = {
61          0.99999999999999709182,
62          57.156235665862923517,
63          -59.597960355475491248,
64          14.136097974741747174,
65          -0.49191381609762019978,
66          .33994649984811888699e-4,
67          .46523628927048575665e-4,
68          -.98374475304879564677e-4,
69          .15808870322491248884e-3,
70          -.21026444172410488319e-3,
71          .21743961811521264320e-3,
72          -.16431810653676389022e-3,
73          .84418223983852743293e-4,
74          -.26190838401581408670e-4,
75          .36899182659531622704e-5,
76      };
77  
78      /** Avoid repeated computation of log of 2 PI in logGamma */
79      private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI);
80  
81      /**
82       * Class contains only static methods.
83       */
84      private InternalGamma() {}
85  
86      /**
87       * Computes the function \( \ln \Gamma(x) \) for \( x > 0 \).
88       *
89       * <p>
90       * For \( x \leq 8 \), the implementation is based on the double precision
91       * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
92       * {@code DGAMLN}. For \( x \geq 8 \), the implementation is based on
93       * </p>
94       *
95       * <ul>
96       * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
97       *     Function</a>, equation (28).</li>
98       * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
99       *     Lanczos Approximation</a>, equations (1) through (5).</li>
100      * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
101      *     the computation of the convergent Lanczos complex Gamma
102      *     approximation</a></li>
103      * </ul>
104      *
105      * @param x Argument.
106      * @return \( \ln \Gamma(x) \), or {@code NaN} if {@code x <= 0}.
107      */
108     public static double logGamma(double x) {
109         // Stripped-down version of the same method defined in "Commons Math":
110         // Unused "if" branches (for when x < 8) have been removed here since
111         // this method is only used (by class "InternalUtils") in order to
112         // compute log(n!) for x > 20.
113 
114         final double sum = lanczos(x);
115         final double tmp = x + LANCZOS_G + 0.5;
116         return (x + 0.5) * Math.log(tmp) - tmp +  HALF_LOG_2_PI + Math.log(sum / x);
117     }
118 
119     /**
120      * Computes the Lanczos approximation used to compute the gamma function.
121      *
122      * <p>
123      * The Lanczos approximation is related to the Gamma function by the
124      * following equation
125      * \[
126      * \Gamma(x) = \sqrt{2\pi} \, \frac{(g + x + \frac{1}{2})^{x + \frac{1}{2}} \, e^{-(g + x + \frac{1}{2})} \, \mathrm{lanczos}(x)}
127      *                                 {x}
128      * \]
129      * where \(g\) is the Lanczos constant.
130      * </p>
131      *
132      * @param x Argument.
133      * @return The Lanczos approximation.
134      *
135      * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos Approximation</a>
136      * equations (1) through (5), and Paul Godfrey's
137      * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation
138      * of the convergent Lanczos complex Gamma approximation</a>
139      */
140     private static double lanczos(final double x) {
141         double sum = 0.0;
142         for (int i = LANCZOS.length - 1; i > 0; --i) {
143             sum += LANCZOS[i] / (x + i);
144         }
145         return sum + LANCZOS[0];
146     }
147 }