commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From celes...@apache.org
Subject svn commit: r1411387 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math3/special/Gamma.java test/java/org/apache/commons/math3/special/GammaTest.java
Date Mon, 19 Nov 2012 20:07:26 GMT
Author: celestin
Date: Mon Nov 19 20:07:25 2012
New Revision: 1411387

URL: http://svn.apache.org/viewvc?rev=1411387&view=rev
Log:
In class Gamma, implementation of (a, b) -> log(Gamma(a + b)).

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math3/special/Gamma.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/GammaTest.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=1411387&r1=1411386&r2=1411387&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 Mon
Nov 19 20:07:25 2012
@@ -19,6 +19,7 @@ package org.apache.commons.math3.special
 import org.apache.commons.math3.exception.MaxCountExceededException;
 import org.apache.commons.math3.exception.NumberIsTooLargeException;
 import org.apache.commons.math3.exception.NumberIsTooSmallException;
+import org.apache.commons.math3.exception.OutOfRangeException;
 import org.apache.commons.math3.util.ContinuedFraction;
 import org.apache.commons.math3.util.FastMath;
 
@@ -701,4 +702,33 @@ public class Gamma {
         }
         return ret;
     }
+
+    /**
+     * Returns the value of log Γ(a + b) for 1 ≤ a, b ≤ 2.
+     *
+     * @param a First argument.
+     * @param b Second argument.
+     * @return the value of {@code log(Gamma(a + b))}.
+     * @throws OutOfRangeException if {@code a} or {@code b} is lower than
+     * {@code 1.0} or greater than {@code 2.0}.
+     */
+    public static double logGammaSum(final double a, final double b)
+        throws OutOfRangeException {
+
+        if ((a < 1.0) || (a > 2.0)) {
+            throw new OutOfRangeException(a, 1.0, 2.0);
+        }
+        if ((b < 1.0) || (b > 2.0)) {
+            throw new OutOfRangeException(b, 1.0, 2.0);
+        }
+
+        final double x = a + b - 2.0;
+        if (x <= 0.25) {
+            return Gamma.logGamma1p(1.0 + x);
+        } else if (x <= 1.25) {
+            return Gamma.logGamma1p(x) + FastMath.log1p(x);
+        } else {
+            return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
+        }
+    }
 }

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/GammaTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/GammaTest.java?rev=1411387&r1=1411386&r2=1411387&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/GammaTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/special/GammaTest.java
Mon Nov 19 20:07:25 2012
@@ -465,7 +465,7 @@ public class GammaTest {
      *     print("{", float(x[i]), ",", float(gamma(x[i])), "},");
      * </pre>
      */
-    private final static double[][] GAMMA_REF = {
+    private static final double[][] GAMMA_REF = {
         { - 19.875 , 4.920331854832504e-18 },
         { - 19.75 , 3.879938752480031e-18 },
         { - 19.625 , 4.323498423815027e-18 },
@@ -954,6 +954,127 @@ public class GammaTest {
         }
     }
 
+    /**
+     * Reference data for the {@link Gamma#logGammaSum(double, double)}
+     * function. This data was generated with the following
+     * <a href="http://maxima.sourceforge.net/">Maxima</a> script.
+     *
+     * <pre>
+     * kill(all);
+     *
+     * fpprec : 64;
+     * gsumln(a, b) := log(gamma(a + b));
+     *
+     * x : [1.0b0, 1.125b0, 1.25b0, 1.375b0, 1.5b0, 1.625b0, 1.75b0, 1.875b0, 2.0b0];
+     *
+     * for i : 1 while i <= length(x) do
+     *   for j : 1 while j <= length(x) do block(
+     *     a : x[i],
+     *     b : x[j],
+     *     print("{", float(a), ",", float(b), ",", float(gsumln(a, b)), "},")
+     *   );
+     * </pre>
+     */
+    private static final double[][] LOG_GAMMA_SUM_REF = {
+        { 1.0 , 1.0 , 0.0 },
+        { 1.0 , 1.125 , .05775985153034387 },
+        { 1.0 , 1.25 , .1248717148923966 },
+        { 1.0 , 1.375 , .2006984603774558 },
+        { 1.0 , 1.5 , .2846828704729192 },
+        { 1.0 , 1.625 , .3763336820249054 },
+        { 1.0 , 1.75 , .4752146669149371 },
+        { 1.0 , 1.875 , .5809359740231859 },
+        { 1.0 , 2.0 , .6931471805599453 },
+        { 1.125 , 1.0 , .05775985153034387 },
+        { 1.125 , 1.125 , .1248717148923966 },
+        { 1.125 , 1.25 , .2006984603774558 },
+        { 1.125 , 1.375 , .2846828704729192 },
+        { 1.125 , 1.5 , .3763336820249054 },
+        { 1.125 , 1.625 , .4752146669149371 },
+        { 1.125 , 1.75 , .5809359740231859 },
+        { 1.125 , 1.875 , .6931471805599453 },
+        { 1.125 , 2.0 , 0.811531653906724 },
+        { 1.25 , 1.0 , .1248717148923966 },
+        { 1.25 , 1.125 , .2006984603774558 },
+        { 1.25 , 1.25 , .2846828704729192 },
+        { 1.25 , 1.375 , .3763336820249054 },
+        { 1.25 , 1.5 , .4752146669149371 },
+        { 1.25 , 1.625 , .5809359740231859 },
+        { 1.25 , 1.75 , .6931471805599453 },
+        { 1.25 , 1.875 , 0.811531653906724 },
+        { 1.25 , 2.0 , .9358019311087253 },
+        { 1.375 , 1.0 , .2006984603774558 },
+        { 1.375 , 1.125 , .2846828704729192 },
+        { 1.375 , 1.25 , .3763336820249054 },
+        { 1.375 , 1.375 , .4752146669149371 },
+        { 1.375 , 1.5 , .5809359740231859 },
+        { 1.375 , 1.625 , .6931471805599453 },
+        { 1.375 , 1.75 , 0.811531653906724 },
+        { 1.375 , 1.875 , .9358019311087253 },
+        { 1.375 , 2.0 , 1.06569589786406 },
+        { 1.5 , 1.0 , .2846828704729192 },
+        { 1.5 , 1.125 , .3763336820249054 },
+        { 1.5 , 1.25 , .4752146669149371 },
+        { 1.5 , 1.375 , .5809359740231859 },
+        { 1.5 , 1.5 , .6931471805599453 },
+        { 1.5 , 1.625 , 0.811531653906724 },
+        { 1.5 , 1.75 , .9358019311087253 },
+        { 1.5 , 1.875 , 1.06569589786406 },
+        { 1.5 , 2.0 , 1.200973602347074 },
+        { 1.625 , 1.0 , .3763336820249054 },
+        { 1.625 , 1.125 , .4752146669149371 },
+        { 1.625 , 1.25 , .5809359740231859 },
+        { 1.625 , 1.375 , .6931471805599453 },
+        { 1.625 , 1.5 , 0.811531653906724 },
+        { 1.625 , 1.625 , .9358019311087253 },
+        { 1.625 , 1.75 , 1.06569589786406 },
+        { 1.625 , 1.875 , 1.200973602347074 },
+        { 1.625 , 2.0 , 1.341414578068493 },
+        { 1.75 , 1.0 , .4752146669149371 },
+        { 1.75 , 1.125 , .5809359740231859 },
+        { 1.75 , 1.25 , .6931471805599453 },
+        { 1.75 , 1.375 , 0.811531653906724 },
+        { 1.75 , 1.5 , .9358019311087253 },
+        { 1.75 , 1.625 , 1.06569589786406 },
+        { 1.75 , 1.75 , 1.200973602347074 },
+        { 1.75 , 1.875 , 1.341414578068493 },
+        { 1.75 , 2.0 , 1.486815578593417 },
+        { 1.875 , 1.0 , .5809359740231859 },
+        { 1.875 , 1.125 , .6931471805599453 },
+        { 1.875 , 1.25 , 0.811531653906724 },
+        { 1.875 , 1.375 , .9358019311087253 },
+        { 1.875 , 1.5 , 1.06569589786406 },
+        { 1.875 , 1.625 , 1.200973602347074 },
+        { 1.875 , 1.75 , 1.341414578068493 },
+        { 1.875 , 1.875 , 1.486815578593417 },
+        { 1.875 , 2.0 , 1.6369886482725 },
+        { 2.0 , 1.0 , .6931471805599453 },
+        { 2.0 , 1.125 , 0.811531653906724 },
+        { 2.0 , 1.25 , .9358019311087253 },
+        { 2.0 , 1.375 , 1.06569589786406 },
+        { 2.0 , 1.5 , 1.200973602347074 },
+        { 2.0 , 1.625 , 1.341414578068493 },
+        { 2.0 , 1.75 , 1.486815578593417 },
+        { 2.0 , 1.875 , 1.6369886482725 },
+        { 2.0 , 2.0 , 1.791759469228055 },
+    };
+
+    @Test
+    public void testLogGammaSum() {
+        final int ulps = 5;
+        for (int i = 0; i < LOG_GAMMA_SUM_REF.length; i++) {
+            final double[] ref = LOG_GAMMA_SUM_REF[i];
+            final double a = ref[0];
+            final double b = ref[1];
+            final double expected = ref[2];
+            final double actual = Gamma.logGammaSum(a, b);
+            final double tol = ulps * FastMath.ulp(expected);
+            final StringBuilder builder = new StringBuilder();
+            builder.append(a).append(", ").append(b);
+            Assert.assertEquals(builder.toString(), expected, actual, tol);
+        }
+    }
+
     private void checkRelativeError(String msg, double expected, double actual, double tolerance)
{
         Assert.assertEquals(msg, expected, actual, FastMath.abs(tolerance * actual));
     }



Mime
View raw message