diff --git a/src/main/java/org/apache/commons/math3/special/Gamma.java b/src/main/java/org/apache/commons/math3/special/Gamma.java index 889791623..69fa26f2d 100644 --- a/src/main/java/org/apache/commons/math3/special/Gamma.java +++ b/src/main/java/org/apache/commons/math3/special/Gamma.java @@ -214,6 +214,68 @@ public class Gamma { /** The constant {@code C13} defined in {@code DGAM1}. */ private static final double INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06; + /** + *

+ * The d0 coefficient of the minimax approximation of the Δ + * function. This function is defined as follows + *

+ *
Δ(x) = log Γ(x) - (x - 0.5) log a + a - 0.5 log 2π,
+ *

+ * The minimax approximation is defined by the following sum + *

+ *
+     *             5
+     *            ====
+     *            \         - 2 n - 1
+     *     Δ(x) =  >    d  x
+     *            /      n
+     *            ====
+     *            n = 0
+     * 
+     * 

+ * see equationS (23) and (25) in Didonato and Morris (1992). + *

+ */ + private static final double D0 = .833333333333333E-01; + + /** + * The d1 coefficent of the minimax approximation of the Δ + * function (see {@link #D0}). + */ + private static final double D1 = -.277777777760991E-02; + + /** + * The d2 coefficent of the minimax approximation of the Δ + * function (see {@link #D0}). + */ + private static final double D2 = .793650666825390E-03; + + /** + * The d3 coefficent of the minimax approximation of the Δ + * function (see {@link #D0}). + */ + private static final double D3 = -.595202931351870E-03; + + /** + * The d4 coefficent of the minimax approximation of the Δ + * function (see {@link #D0}). + */ + private static final double D4 = .837308034031215E-03; + + /** + * The d5 coefficent of the minimax approximation of the Δ + * function (see {@link #D0}). + */ + private static final double D5 = -.165322962780713E-02; + /* + * NOTA: the value of d0 published in Didonato and Morris (1992), eq. (25) + * and the value implemented in the NSWC library are NOT EQUAL + * - in Didonato and Morris (1992) D5 = -.125322962780713E-02, + * - while in the NSWC library D5 = -.165322962780713E-02. + * Checking the value of algdiv(1.0, 8.0)}, it seems that the second value + * leads to the smallest error. This is the one which is implemented here. + */ + /** * Default constructor. Prohibit instantiation. */ @@ -734,4 +796,70 @@ public class Gamma { return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x)); } } + + /** + * Returns the value of log[Γ(b) / Γ(a + b)] for a ≥ 0 and b ≥ 8. The + * present implementation is based on the double precision implementation in + * the NSWC Library of Mathematics Subroutines, {@code ALGDIV}. + * + * @param a First argument. + * @param b Second argument. + * @return the value of {@code log(Gamma(b) / Gamma(a + b))}. + * @throws NumberIsTooSmallException if {@code a < 0.0} or {@code b < 8.0}. + */ + public static final double logGammaMinusLogGammaSum(final double a, + final double b) + throws NumberIsTooSmallException { + + if (a < 0.0) { + throw new NumberIsTooSmallException(a, 0.0, true); + } + if (b < 8.0) { + throw new NumberIsTooSmallException(b, 8.0, true); + } + + /* + * p = a / (a + b), q = b / (a + b), d = a + b - 0.5 + */ + final double p; + final double q; + final double d; + if (a <= b) { + final double h = a / b; + p = h / (1.0 + h); + q = 1.0 / (1.0 + h); + d = b + (a - 0.5); + } else { + final double h = b / a; + p = 1.0 / (1.0 + h); + q = h / (1.0 + h); + d = a + (b - 0.5); + } + /* + * s_n = 1 + q + ... + q^(n - 1) + */ + final double q2 = q * q; + final double s3 = 1.0 + (q + q2); + final double s5 = 1.0 + (q + q2 * s3); + final double s7 = 1.0 + (q + q2 * s5); + final double s9 = 1.0 + (q + q2 * s7); + final double s11 = 1.0 + (q + q2 * s9); + /* + * w = Δ(b) - Δ(a + b) + */ + final double tmp = 1.0 / b; + final double t = tmp * tmp; + double w = D5 * s11; + w = D4 * s9 + t * w; + w = D3 * s7 + t * w; + w = D2 * s5 + t * w; + w = D1 * s3 + t * w; + w = D0 + t * w; + w *= p / b; + + final double u = d * FastMath.log1p(a / b); + final double v = a * (FastMath.log(b) - 1.0); + + return u <= v ? (w - u) - v : (w - v) - u; + } } diff --git a/src/test/java/org/apache/commons/math3/special/GammaTest.java b/src/test/java/org/apache/commons/math3/special/GammaTest.java index 2db68194b..c121fd78b 100644 --- a/src/test/java/org/apache/commons/math3/special/GammaTest.java +++ b/src/test/java/org/apache/commons/math3/special/GammaTest.java @@ -1075,6 +1075,156 @@ public class GammaTest { } } + private static final double[][] LOG_GAMMA_MINUS_LOG_GAMMA_SUM_REF = { + { 0.0 , 8.0 , 0.0 }, + { 0.0 , 9.0 , 0.0 }, + { 0.0 , 10.0 , 0.0 }, + { 0.0 , 11.0 , 0.0 }, + { 0.0 , 12.0 , 0.0 }, + { 0.0 , 13.0 , 0.0 }, + { 0.0 , 14.0 , 0.0 }, + { 0.0 , 15.0 , 0.0 }, + { 0.0 , 16.0 , 0.0 }, + { 0.0 , 17.0 , 0.0 }, + { 0.0 , 18.0 , 0.0 }, + { 1.0 , 8.0 , - 2.079441541679836 }, + { 1.0 , 9.0 , - 2.19722457733622 }, + { 1.0 , 10.0 , - 2.302585092994046 }, + { 1.0 , 11.0 , - 2.397895272798371 }, + { 1.0 , 12.0 , - 2.484906649788 }, + { 1.0 , 13.0 , - 2.564949357461537 }, + { 1.0 , 14.0 , - 2.639057329615258 }, + { 1.0 , 15.0 , - 2.70805020110221 }, + { 1.0 , 16.0 , - 2.772588722239781 }, + { 1.0 , 17.0 , - 2.833213344056216 }, + { 1.0 , 18.0 , - 2.890371757896165 }, + { 2.0 , 8.0 , - 4.276666119016055 }, + { 2.0 , 9.0 , - 4.499809670330265 }, + { 2.0 , 10.0 , - 4.700480365792417 }, + { 2.0 , 11.0 , - 4.882801922586371 }, + { 2.0 , 12.0 , - 5.049856007249537 }, + { 2.0 , 13.0 , - 5.204006687076795 }, + { 2.0 , 14.0 , - 5.347107530717468 }, + { 2.0 , 15.0 , - 5.480638923341991 }, + { 2.0 , 16.0 , - 5.605802066295998 }, + { 2.0 , 17.0 , - 5.723585101952381 }, + { 2.0 , 18.0 , - 5.834810737062605 }, + { 3.0 , 8.0 , - 6.579251212010101 }, + { 3.0 , 9.0 , - 6.897704943128636 }, + { 3.0 , 10.0 , - 7.185387015580416 }, + { 3.0 , 11.0 , - 7.447751280047908 }, + { 3.0 , 12.0 , - 7.688913336864796 }, + { 3.0 , 13.0 , - 7.912056888179006 }, + { 3.0 , 14.0 , - 8.11969625295725 }, + { 3.0 , 15.0 , - 8.313852267398207 }, + { 3.0 , 16.0 , - 8.496173824192162 }, + { 3.0 , 17.0 , - 8.668024081118821 }, + { 3.0 , 18.0 , - 8.830543010616596 }, + { 4.0 , 8.0 , - 8.977146484808472 }, + { 4.0 , 9.0 , - 9.382611592916636 }, + { 4.0 , 10.0 , - 9.750336373041954 }, + { 4.0 , 11.0 , - 10.08680860966317 }, + { 4.0 , 12.0 , - 10.39696353796701 }, + { 4.0 , 13.0 , - 10.68464561041879 }, + { 4.0 , 14.0 , - 10.95290959701347 }, + { 4.0 , 15.0 , - 11.20422402529437 }, + { 4.0 , 16.0 , - 11.4406128033586 }, + { 4.0 , 17.0 , - 11.66375635467281 }, + { 4.0 , 18.0 , - 11.87506544834002 }, + { 5.0 , 8.0 , - 11.46205313459647 }, + { 5.0 , 9.0 , - 11.94756095037817 }, + { 5.0 , 10.0 , - 12.38939370265721 }, + { 5.0 , 11.0 , - 12.79485881076538 }, + { 5.0 , 12.0 , - 13.16955226020679 }, + { 5.0 , 13.0 , - 13.517858954475 }, + { 5.0 , 14.0 , - 13.84328135490963 }, + { 5.0 , 15.0 , - 14.14866300446081 }, + { 5.0 , 16.0 , - 14.43634507691259 }, + { 5.0 , 17.0 , - 14.70827879239624 }, + { 5.0 , 18.0 , - 14.96610790169833 }, + { 6.0 , 8.0 , - 14.02700249205801 }, + { 6.0 , 9.0 , - 14.58661827999343 }, + { 6.0 , 10.0 , - 15.09744390375942 }, + { 6.0 , 11.0 , - 15.56744753300516 }, + { 6.0 , 12.0 , - 16.002765604263 }, + { 6.0 , 13.0 , - 16.40823071237117 }, + { 6.0 , 14.0 , - 16.78772033407607 }, + { 6.0 , 15.0 , - 17.14439527801481 }, + { 6.0 , 16.0 , - 17.48086751463602 }, + { 6.0 , 17.0 , - 17.79932124575455 }, + { 6.0 , 18.0 , - 18.10160211762749 }, + { 7.0 , 8.0 , - 16.66605982167327 }, + { 7.0 , 9.0 , - 17.29466848109564 }, + { 7.0 , 10.0 , - 17.8700326259992 }, + { 7.0 , 11.0 , - 18.40066087706137 }, + { 7.0 , 12.0 , - 18.89313736215917 }, + { 7.0 , 13.0 , - 19.35266969153761 }, + { 7.0 , 14.0 , - 19.78345260763006 }, + { 7.0 , 15.0 , - 20.18891771573823 }, + { 7.0 , 16.0 , - 20.57190996799433 }, + { 7.0 , 17.0 , - 20.9348154616837 }, + { 7.0 , 18.0 , - 21.27965594797543 }, + { 8.0 , 8.0 , - 19.37411002277548 }, + { 8.0 , 9.0 , - 20.06725720333542 }, + { 8.0 , 10.0 , - 20.70324597005542 }, + { 8.0 , 11.0 , - 21.29103263495754 }, + { 8.0 , 12.0 , - 21.83757634132561 }, + { 8.0 , 13.0 , - 22.3484019650916 }, + { 8.0 , 14.0 , - 22.82797504535349 }, + { 8.0 , 15.0 , - 23.27996016909654 }, + { 8.0 , 16.0 , - 23.70740418392348 }, + { 8.0 , 17.0 , - 24.11286929203165 }, + { 8.0 , 18.0 , - 24.49853177284363 }, + { 9.0 , 8.0 , - 22.14669874501526 }, + { 9.0 , 9.0 , - 22.90047054739164 }, + { 9.0 , 10.0 , - 23.59361772795159 }, + { 9.0 , 11.0 , - 24.23547161412398 }, + { 9.0 , 12.0 , - 24.8333086148796 }, + { 9.0 , 13.0 , - 25.39292440281502 }, + { 9.0 , 14.0 , - 25.9190174987118 }, + { 9.0 , 15.0 , - 26.41545438502569 }, + { 9.0 , 16.0 , - 26.88545801427143 }, + { 9.0 , 17.0 , - 27.33174511689985 }, + { 9.0 , 18.0 , - 27.75662831086511 }, + { 10.0 , 8.0 , - 24.97991208907148 }, + { 10.0 , 9.0 , - 25.7908423052878 }, + { 10.0 , 10.0 , - 26.53805670711802 }, + { 10.0 , 11.0 , - 27.23120388767797 }, + { 10.0 , 12.0 , - 27.87783105260302 }, + { 10.0 , 13.0 , - 28.48396685617334 }, + { 10.0 , 14.0 , - 29.05451171464095 }, + { 10.0 , 15.0 , - 29.59350821537364 }, + { 10.0 , 16.0 , - 30.10433383913963 }, + { 10.0 , 17.0 , - 30.58984165492133 }, + { 10.0 , 18.0 , - 31.05246517686944 }, + }; + + @Test + public void testLogGammaMinusLogGammaSum() { + final int ulps = 4; + for (int i = 0; i < LOG_GAMMA_MINUS_LOG_GAMMA_SUM_REF.length; i++) { + final double[] ref = LOG_GAMMA_MINUS_LOG_GAMMA_SUM_REF[i]; + final double a = ref[0]; + final double b = ref[1]; + final double expected = ref[2]; + final double actual = Gamma.logGammaMinusLogGammaSum(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); + } + } + + @Test(expected = NumberIsTooSmallException.class) + public void testLogGammaMinusLogGammaSumPrecondition1() { + Gamma.logGammaMinusLogGammaSum(-1.0, 8.0); + } + + @Test(expected = NumberIsTooSmallException.class) + public void testLogGammaMinusLogGammaSumPrecondition2() { + Gamma.logGammaMinusLogGammaSum(1.0, 7.0); + } + private void checkRelativeError(String msg, double expected, double actual, double tolerance) { Assert.assertEquals(msg, expected, actual, FastMath.abs(tolerance * actual)); }