MATH-738: in class o.a.c.m3.special.Gamma, implemented function
logGammaMinusLogGammaSum(double, double), which computes accurately log(Gamma(b)) - log(Gamma(a + b)) for a >= 0 and b >= 8. Based on the NSWC library. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1413366 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Sebastien Brisard
parent
1d25ac89e3
commit
27178b80c4
|
|
@ -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;
|
||||
|
||||
/**
|
||||
* <p>
|
||||
* The d<sub>0</sub> coefficient of the minimax approximation of the Δ
|
||||
* function. This function is defined as follows
|
||||
* </p>
|
||||
* <center>Δ(x) = log Γ(x) - (x - 0.5) log a + a - 0.5 log 2π,</center>
|
||||
* <p>
|
||||
* The minimax approximation is defined by the following sum
|
||||
* </p>
|
||||
* <pre>
|
||||
* 5
|
||||
* ====
|
||||
* \ - 2 n - 1
|
||||
* Δ(x) = > d x
|
||||
* / n
|
||||
* ====
|
||||
* n = 0
|
||||
* <pre>
|
||||
* <p>
|
||||
* see equationS (23) and (25) in Didonato and Morris (1992).
|
||||
* </p>
|
||||
*/
|
||||
private static final double D0 = .833333333333333E-01;
|
||||
|
||||
/**
|
||||
* The d<sub>1</sub> coefficent of the minimax approximation of the Δ
|
||||
* function (see {@link #D0}).
|
||||
*/
|
||||
private static final double D1 = -.277777777760991E-02;
|
||||
|
||||
/**
|
||||
* The d<sub>2</sub> coefficent of the minimax approximation of the Δ
|
||||
* function (see {@link #D0}).
|
||||
*/
|
||||
private static final double D2 = .793650666825390E-03;
|
||||
|
||||
/**
|
||||
* The d<sub>3</sub> coefficent of the minimax approximation of the Δ
|
||||
* function (see {@link #D0}).
|
||||
*/
|
||||
private static final double D3 = -.595202931351870E-03;
|
||||
|
||||
/**
|
||||
* The d<sub>4</sub> coefficent of the minimax approximation of the Δ
|
||||
* function (see {@link #D0}).
|
||||
*/
|
||||
private static final double D4 = .837308034031215E-03;
|
||||
|
||||
/**
|
||||
* The d<sub>5</sub> 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 <em>NSWC Library of Mathematics Subroutines</em>, {@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;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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));
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue