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 * </h3>
24 *
25 * <p>
26 * This is a utility class that provides computation methods related to the
27 * Γ (Gamma) family of functions.
28 * </p>
29 */
30 final class InternalGamma { // Class is package-private on purpose; do not make it public.
31 /**
32 * Constant \( g = \frac{607}{128} \) in the Lanczos approximation.
33 */
34 public static final double LANCZOS_G = 607.0 / 128.0;
35
36 /** Lanczos coefficients. */
37 private static final double[] LANCZOS_COEFFICIENTS = {
38 0.99999999999999709182,
39 57.156235665862923517,
40 -59.597960355475491248,
41 14.136097974741747174,
42 -0.49191381609762019978,
43 .33994649984811888699e-4,
44 .46523628927048575665e-4,
45 -.98374475304879564677e-4,
46 .15808870322491248884e-3,
47 -.21026444172410488319e-3,
48 .21743961811521264320e-3,
49 -.16431810653676389022e-3,
50 .84418223983852743293e-4,
51 -.26190838401581408670e-4,
52 .36899182659531622704e-5,
53 };
54
55 /** Avoid repeated computation of log of 2 PI in logGamma. */
56 private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI);
57
58 /**
59 * Class contains only static methods.
60 */
61 private InternalGamma() {}
62
63 /**
64 * Computes the function \( \ln \Gamma(x) \) for \( x > 0 \).
65 *
66 * <p>
67 * For \( x \leq 8 \), the implementation is based on the double precision
68 * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
69 * {@code DGAMLN}. For \( x \geq 8 \), the implementation is based on
70 * </p>
71 *
72 * <ul>
73 * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
74 * Function</a>, equation (28).</li>
75 * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
76 * Lanczos Approximation</a>, equations (1) through (5).</li>
77 * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
78 * the computation of the convergent Lanczos complex Gamma
79 * approximation</a></li>
80 * </ul>
81 *
82 * @param x Argument.
83 * @return \( \ln \Gamma(x) \), or {@code NaN} if {@code x <= 0}.
84 */
85 public static double logGamma(double x) {
86 // Stripped-down version of the same method defined in "Commons Math":
87 // Unused "if" branches (for when x < 8) have been removed here since
88 // this method is only used (by class "InternalUtils") in order to
89 // compute log(n!) for x > 20.
90
91 final double sum = lanczos(x);
92 final double tmp = x + LANCZOS_G + 0.5;
93 return (x + 0.5) * Math.log(tmp) - tmp + HALF_LOG_2_PI + Math.log(sum / x);
94 }
95
96 /**
97 * Computes the Lanczos approximation used to compute the gamma function.
98 *
99 * <p>
100 * The Lanczos approximation is related to the Gamma function by the
101 * following equation
102 * \[
103 * \Gamma(x) = \sqrt{2\pi} \, \frac{(g + x + \frac{1}{2})^{x + \frac{1}{2}} \, e^{-(g + x + \frac{1}{2})} \, \mathrm{lanczos}(x)}
104 * {x}
105 * \]
106 * where \(g\) is the Lanczos constant.
107 * </p>
108 *
109 * @param x Argument.
110 * @return The Lanczos approximation.
111 *
112 * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos Approximation</a>
113 * equations (1) through (5), and Paul Godfrey's
114 * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation
115 * of the convergent Lanczos complex Gamma approximation</a>
116 */
117 private static double lanczos(final double x) {
118 double sum = 0.0;
119 for (int i = LANCZOS_COEFFICIENTS.length - 1; i > 0; --i) {
120 sum += LANCZOS_COEFFICIENTS[i] / (x + i);
121 }
122 return sum + LANCZOS_COEFFICIENTS[0];
123 }
124 }