in commons-numbers-gamma/src/main/java/org/apache/commons/numbers/gamma/SpecialMath.java [157:218]
private static double log1pmxSmall(double x, double a) {
// Use a direct Taylor series:
// ln(1 + x) = x - x^2/2 + x^3/3 - x^4/4 + ...
// Reverse the summation (small to large) for a marginal increase in precision.
// To stop the Taylor series the next term must be less than 1 ulp from the
// answer.
// x^n/n < |log(1+x)-x| * eps
// eps = machine epsilon = 2^-53
// x^n < |log(1+x)-x| * eps
// n < (log(|log(1+x)-x|) + log(eps)) / log(x)
// In practice this is a conservative limit.
final double x2 = x * x;
if (a < TWO_POW_M53) {
// Below machine epsilon. Addition of x^3/3 is not possible.
// Subtract from zero to prevent creating -0.0 for x=0.
return 0 - x2 / 2;
}
final double x4 = x2 * x2;
// +/-9.5367431640625e-07: log1pmx = -4.547470617660916e-13 :
// -4.5474764000725028e-13
// n = 4.69
if (a < TWO_POW_M20) {
// n=5
return x * x4 / 5 -
x4 / 4 +
x * x2 / 3 -
x2 / 2;
}
// +/-2.44140625E-4: log1pmx = -2.9797472637290841e-08 : -2.9807173914456693e-08
// n = 6.49
if (a < TWO_POW_M12) {
// n=7
return x * x2 * x4 / 7 -
x2 * x4 / 6 +
x * x4 / 5 -
x4 / 4 +
x * x2 / 3 -
x2 / 2;
}
// Assume |x| < 2^-6
// +/-0.015625: log1pmx = -0.00012081346403474586 : -0.00012335696813916864
// n = 10.9974
// n=11
final double x8 = x4 * x4;
return x * x2 * x8 / 11 -
x2 * x8 / 10 +
x * x8 / 9 -
x8 / 8 +
x * x2 * x4 / 7 -
x2 * x4 / 6 +
x * x4 / 5 -
x4 / 4 +
x * x2 / 3 -
x2 / 2;
}