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) 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 }