in commons-math-core/src/main/java/org/apache/commons/math4/core/jdkmath/AccurateMath.java [1133:1353]
private static double log(final double x, final double[] hiPrec) {
if (x == 0) { // Handle special case of +0/-0
return Double.NEGATIVE_INFINITY;
}
long bits = Double.doubleToRawLongBits(x);
/* Handle special cases of negative input, and NaN */
if (((bits & 0x8000000000000000L) != 0 || Double.isNaN(x)) && x != 0.0) {
if (hiPrec != null) {
hiPrec[0] = Double.NaN;
}
return Double.NaN;
}
/* Handle special cases of Positive infinity. */
if (x == Double.POSITIVE_INFINITY) {
if (hiPrec != null) {
hiPrec[0] = Double.POSITIVE_INFINITY;
}
return Double.POSITIVE_INFINITY;
}
/* Extract the exponent */
int exp = (int) (bits >> 52) - 1023;
if ((bits & 0x7ff0000000000000L) == 0) {
// Subnormal!
if (x == 0) {
// Zero
if (hiPrec != null) {
hiPrec[0] = Double.NEGATIVE_INFINITY;
}
return Double.NEGATIVE_INFINITY;
}
/* Normalize the subnormal number. */
bits <<= 1;
while ((bits & 0x0010000000000000L) == 0) {
--exp;
bits <<= 1;
}
}
if ((exp == -1 || exp == 0) && x < 1.01 && x > 0.99 && hiPrec == null) {
/* The normal method doesn't work well in the range [0.99, 1.01], so call do a straight
polynomial expansion in higer precision. */
/* Compute x - 1.0 and split it */
double xa = x - 1.0;
double xb = xa - x + 1.0;
double tmp = xa * HEX_40000000;
double aa = xa + tmp - tmp;
double ab = xa - aa;
xa = aa;
xb = ab;
final double[] lnCoef_last = LN_QUICK_COEF[LN_QUICK_COEF.length - 1];
double ya = lnCoef_last[0];
double yb = lnCoef_last[1];
for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
/* Multiply a = y * x */
aa = ya * xa;
ab = ya * xb + yb * xa + yb * xb;
/* split, so now y = a */
tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
yb = aa - ya + ab;
/* Add a = y + lnQuickCoef */
final double[] lnCoef_i = LN_QUICK_COEF[i];
aa = ya + lnCoef_i[0];
ab = yb + lnCoef_i[1];
/* Split y = a */
tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
yb = aa - ya + ab;
}
/* Multiply a = y * x */
aa = ya * xa;
ab = ya * xb + yb * xa + yb * xb;
/* split, so now y = a */
tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
yb = aa - ya + ab;
return ya + yb;
}
// lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2)
final double[] lnm = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
/*
double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L);
epsilon -= 1.0;
*/
// y is the most significant 10 bits of the mantissa
//double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L);
//double epsilon = (x - y) / y;
final double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
double lnza = 0.0;
double lnzb = 0.0;
if (hiPrec != null) {
/* split epsilon -> x */
double tmp = epsilon * HEX_40000000;
double aa = epsilon + tmp - tmp;
double ab = epsilon - aa;
double xa = aa;
double xb = ab;
/* Need a more accurate epsilon, so adjust the division. */
final double numer = bits & 0x3ffffffffffL;
final double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
aa = numer - xa * denom - xb * denom;
xb += aa / denom;
/* Remez polynomial evaluation */
final double[] lnCoef_last = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length - 1];
double ya = lnCoef_last[0];
double yb = lnCoef_last[1];
for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
/* Multiply a = y * x */
aa = ya * xa;
ab = ya * xb + yb * xa + yb * xb;
/* split, so now y = a */
tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
yb = aa - ya + ab;
/* Add a = y + lnHiPrecCoef */
final double[] lnCoef_i = LN_HI_PREC_COEF[i];
aa = ya + lnCoef_i[0];
ab = yb + lnCoef_i[1];
/* Split y = a */
tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
yb = aa - ya + ab;
}
/* Multiply a = y * x */
aa = ya * xa;
ab = ya * xb + yb * xa + yb * xb;
/* split, so now lnz = a */
/*
tmp = aa * 1073741824.0;
lnza = aa + tmp - tmp;
lnzb = aa - lnza + ab;
*/
lnza = aa + ab;
lnzb = -(lnza - aa - ab);
} else {
/* High precision not required. Eval Remez polynomial
using standard double precision */
lnza = -0.16624882440418567;
lnza = lnza * epsilon + 0.19999954120254515;
lnza = lnza * epsilon + -0.2499999997677497;
lnza = lnza * epsilon + 0.3333333333332802;
lnza = lnza * epsilon + -0.5;
lnza = lnza * epsilon + 1.0;
lnza *= epsilon;
}
/* Relative sizes:
* lnzb [0, 2.33E-10]
* lnm[1] [0, 1.17E-7]
* ln2B*exp [0, 1.12E-4]
* lnza [0, 9.7E-4]
* lnm[0] [0, 0.692]
* ln2A*exp [0, 709]
*/
/* Compute the following sum:
* lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
*/
//return lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
double a = LN_2_A * exp;
double b = 0.0;
double c = a + lnm[0];
double d = -(c - a - lnm[0]);
a = c;
b += d;
c = a + lnza;
d = -(c - a - lnza);
a = c;
b += d;
c = a + LN_2_B * exp;
d = -(c - a - LN_2_B * exp);
a = c;
b += d;
c = a + lnm[1];
d = -(c - a - lnm[1]);
a = c;
b += d;
c = a + lnzb;
d = -(c - a - lnzb);
a = c;
b += d;
if (hiPrec != null) {
hiPrec[0] = a;
hiPrec[1] = b;
}
return a + b;
}