commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From celes...@apache.org
Subject svn commit: r1336483 - /commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java
Date Thu, 10 May 2012 05:08:02 GMT
Author: celestin
Date: Thu May 10 05:08:01 2012
New Revision: 1336483

URL: http://svn.apache.org/viewvc?rev=1336483&view=rev
Log:
In o.a.c.m3.special.Gamma, exposed the Lanczos approximation (MATH-753).

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java?rev=1336483&r1=1336482&r2=1336483&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java Thu
May 10 05:08:01 2012
@@ -32,6 +32,13 @@ public class Gamma {
      * @since 2.0
      */
     public static final double GAMMA = 0.577215664901532860606512090082;
+
+    /**
+     * The value of the {@code g} constant in the Lanczos approximation, see
+     * {@link #lanczos(double)}.
+     */
+    public static final double LANCZOS_G = 607.0 / 128.0;
+
     /** Maximum allowed numerical error. */
     private static final double DEFAULT_EPSILON = 10e-15;
     /** Lanczos coefficients */
@@ -89,13 +96,7 @@ public class Gamma {
             ret = Double.NaN;
         } else {
             double g = 607.0 / 128.0;
-
-            double sum = 0.0;
-            for (int i = LANCZOS.length - 1; i > 0; --i) {
-                sum = sum + (LANCZOS[i] / (x + i));
-            }
-            sum = sum + LANCZOS[0];
-
+            double sum = lanczos(x);
             double tmp = x + g + .5;
             ret = ((x + .5) * FastMath.log(tmp)) - tmp +
                 HALF_LOG_2_PI + FastMath.log(sum / x);
@@ -325,4 +326,31 @@ public class Gamma {
 
         return trigamma(x + 1) + 1 / (x * x);
     }
+
+    /**
+     * <p>
+     * Returns the Lanczos approximation used to compute the gamma function.
+     * The Lanczos approximation is related to the Gamma function by the
+     * following equation
+     * <center>
+     * {@code gamma(x) = sqrt(2 * pi) / x * (x + g + 0.5) ^ (x + 0.5)
+     *                   * exp(-x - g - 0.5) * lanczos(x)},
+     * </center>
+     * where {@code g} is a constant, returned by {@link #getLanczosG()}.
+     * </p>
+     *
+     * @param x the argument
+     * @return the Lanczos approximation
+     * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos
Approximation</a>,
+     * equations (1) through (5), and Paul Godfrey's
+     * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation
+     * of the convergent Lanczos complex Gamma approximation</a>
+     */
+    public static double lanczos(final double x) {
+        double sum = 0.0;
+        for (int i = LANCZOS.length - 1; i > 0; --i) {
+            sum = sum + (LANCZOS[i] / (x + i));
+        }
+        return sum + LANCZOS[0];
+    }
 }



Mime
View raw message