/* * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ package org.apache.commons.math.special; import org.apache.commons.math.MathException; import org.apache.commons.math.MaxIterationsExceededException; import org.apache.commons.math.util.ContinuedFraction; import org.apache.commons.math.util.FastMath; /** * This is a utility class that provides computation methods related to the * Gamma family of functions. * * @version $Revision: 1042510 $ $Date: 2010-12-06 02:54:18 +0100 (lun. 06 déc. 2010) $ */ public class Gamma { /** * Euler-Mascheroni constant * @since 2.0 */ public static final double GAMMA = 0.577215664901532860606512090082; /** Maximum allowed numerical error. */ private static final double DEFAULT_EPSILON = 10e-15; /** Lanczos coefficients */ private static final double[] LANCZOS = { 0.99999999999999709182, 57.156235665862923517, -59.597960355475491248, 14.136097974741747174, -0.49191381609762019978, .33994649984811888699e-4, .46523628927048575665e-4, -.98374475304879564677e-4, .15808870322491248884e-3, -.21026444172410488319e-3, .21743961811521264320e-3, -.16431810653676389022e-3, .84418223983852743293e-4, -.26190838401581408670e-4, .36899182659531622704e-5, }; /** Avoid repeated computation of log of 2 PI in logGamma */ private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(2.0 * FastMath.PI); // limits for switching algorithm in digamma /** C limit. */ private static final double C_LIMIT = 49; /** S limit. */ private static final double S_LIMIT = 1e-5; /** * Default constructor. Prohibit instantiation. */ private Gamma() { super(); } /** * Returns the natural logarithm of the gamma function Γ(x). * * The implementation of this method is based on: *
Computes the digamma function of x.
* *This is an independently written implementation of the algorithm described in * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976.
* *Some of the constants have been changed to increase accuracy at the moderate expense * of run-time. The result should be accurate to within 10^-8 absolute tolerance for * x >= 10^-5 and within 10^-8 relative tolerance for x > 0.
* *Performance for large negative values of x will be quite expensive (proportional to * |x|). Accuracy for negative values of x should be about 10^-8 absolute for results * less than 10^5 and 10^-8 relative for results larger than that.
* * @param x the argument * @return digamma(x) to within 10-8 relative or absolute error whichever is smaller * @see Digamma at wikipedia * @see Bernardo's original article * @since 2.0 */ public static double digamma(double x) { if (x > 0 && x <= S_LIMIT) { // use method 5 from Bernardo AS103 // accurate to O(x) return -GAMMA - 1 / x; } if (x >= C_LIMIT) { // use method 4 (accurate to O(1/x^8) double inv = 1 / (x * x); // 1 1 1 1 // log(x) - --- - ------ + ------- - ------- // 2 x 12 x^2 120 x^4 252 x^6 return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252)); } return digamma(x + 1) - 1 / x; } /** *Computes the trigamma function of x. This function is derived by taking the derivative of * the implementation of digamma.
* * @param x the argument * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller * @see Trigamma at wikipedia * @see Gamma#digamma(double) * @since 2.0 */ public static double trigamma(double x) { if (x > 0 && x <= S_LIMIT) { return 1 / (x * x); } if (x >= C_LIMIT) { double inv = 1 / (x * x); // 1 1 1 1 1 // - + ---- + ---- - ----- + ----- // x 2 3 5 7 // 2 x 6 x 30 x 42 x return 1 / x + inv / 2 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv / 42)); } return trigamma(x + 1) + 1 / (x * x); } }