You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by se...@apache.org on 2011/01/21 05:02:30 UTC
svn commit: r1061620 -
/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java
Author: sebb
Date: Fri Jan 21 04:02:29 2011
New Revision: 1061620
URL: http://svn.apache.org/viewvc?rev=1061620&view=rev
Log:
MATH-476 FastMath code contains 'magic' numbers
Extracted "splitter" value as a constant
Modified:
commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java
Modified: commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java?rev=1061620&r1=1061619&r2=1061620&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java (original)
+++ commons/proper/math/branches/MATH_2_X/src/main/java/org/apache/commons/math/util/FastMath.java Fri Jan 21 04:02:29 2011
@@ -161,6 +161,11 @@ public class FastMath {
1.2599210498948732,
1.5874010519681994 };
+ /**
+ * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
+ */
+ private static final double HEX_40000000 = 1073741824.0;
+
// Initialize tables
static {
int i;
@@ -245,13 +250,13 @@ public class FastMath {
double ya = hiPrec[0] + hiPrec[1];
double yb = -(ya - hiPrec[0] - hiPrec[1]);
- double temp = ya * 1073741824.0;
+ double temp = ya * HEX_40000000;
double yaa = ya + temp - temp;
double yab = ya - yaa;
// recip = 1/y
double recip = 1.0/ya;
- temp = recip * 1073741824.0;
+ temp = recip * HEX_40000000;
double recipa = recip + temp - temp;
double recipb = recip - recipa;
@@ -309,13 +314,13 @@ public class FastMath {
double ya = hiPrec[0] + hiPrec[1];
double yb = -(ya - hiPrec[0] - hiPrec[1]);
- double temp = ya * 1073741824.0;
+ double temp = ya * HEX_40000000;
double yaa = ya + temp - temp;
double yab = ya - yaa;
// recip = 1/y
double recip = 1.0/ya;
- temp = recip * 1073741824.0;
+ temp = recip * HEX_40000000;
double recipa = recip + temp - temp;
double recipb = recip - recipa;
@@ -350,11 +355,11 @@ public class FastMath {
double denomr = 1.0 / denom;
double denomb = -(denom - 1.0 - ya) + yb;
double ratio = ya * denomr;
- double temp = ratio * 1073741824.0;
+ double temp = ratio * HEX_40000000;
double ra = ratio + temp - temp;
double rb = ratio - ra;
- temp = denom * 1073741824.0;
+ temp = denom * HEX_40000000;
double za = denom + temp - temp;
double zb = denom - za;
@@ -434,13 +439,13 @@ public class FastMath {
db += -(temp - da - yb);
da = temp;
- temp = da * 1073741824.0;
+ temp = da * HEX_40000000;
double daa = da + temp - temp;
double dab = da - daa;
// ratio = na/da
double ratio = na/da;
- temp = ratio * 1073741824.0;
+ temp = ratio * HEX_40000000;
double ratioa = ratio + temp - temp;
double ratiob = ratio - ratioa;
@@ -473,13 +478,13 @@ public class FastMath {
db += -(temp - da - yb);
da = temp;
- temp = da * 1073741824.0;
+ temp = da * HEX_40000000;
double daa = da + temp - temp;
double dab = da - daa;
// ratio = na/da
double ratio = na/da;
- temp = ratio * 1073741824.0;
+ temp = ratio * HEX_40000000;
double ratioa = ratio + temp - temp;
double ratiob = ratio - ratioa;
@@ -814,7 +819,7 @@ public class FastMath {
tempB = -(temp - tempA - tempB);
tempA = temp;
- temp = tempA * 1073741824.0;
+ temp = tempA * HEX_40000000;
baseA = tempA + temp - temp;
baseB = tempB + (tempA - baseA);
@@ -835,7 +840,7 @@ public class FastMath {
zb = -(temp - za - zb);
za = temp;
- temp = za * 1073741824.0;
+ temp = za * HEX_40000000;
temp = za + temp - temp;
zb += za - temp;
za = temp;
@@ -882,11 +887,11 @@ public class FastMath {
double denomr = 1.0 / denom;
double denomb = -(denom - 1.0 - ya) + yb;
double ratio = ya * denomr;
- temp = ratio * 1073741824.0;
+ temp = ratio * HEX_40000000;
final double ra = ratio + temp - temp;
double rb = ratio - ra;
- temp = denom * 1073741824.0;
+ temp = denom * HEX_40000000;
za = denom + temp - temp;
zb = denom - za;
@@ -960,12 +965,12 @@ public class FastMath {
*/
private static void split(final double d, final double split[]) {
if (d < 8e298 && d > -8e298) {
- final double a = d * 1073741824.0;
+ final double a = d * HEX_40000000;
split[0] = (d + a) - a;
split[1] = d - split[0];
} else {
final double a = d * 9.31322574615478515625E-10;
- split[0] = (d + a - d) * 1073741824.0;
+ split[0] = (d + a - d) * HEX_40000000;
split[1] = d - split[0];
}
}
@@ -979,12 +984,12 @@ public class FastMath {
final double d = -(c - a[0] - a[1]);
if (c < 8e298 && c > -8e298) {
- double z = c * 1073741824.0;
+ double z = c * HEX_40000000;
a[0] = (c + z) - z;
a[1] = c - a[0] + d;
} else {
double z = c * 9.31322574615478515625E-10;
- a[0] = (c + z - c) * 1073741824.0;
+ a[0] = (c + z - c) * HEX_40000000;
a[1] = c - a[0] + d;
}
}
@@ -1235,7 +1240,7 @@ public class FastMath {
/* Compute x - 1.0 and split it */
double xa = x - 1.0;
double xb = xa - x + 1.0;
- double tmp = xa * 1073741824.0;
+ double tmp = xa * HEX_40000000;
double aa = xa + tmp - tmp;
double ab = xa - aa;
xa = aa;
@@ -1249,7 +1254,7 @@ public class FastMath {
aa = ya * xa;
ab = ya * xb + yb * xa + yb * xb;
/* split, so now y = a */
- tmp = aa * 1073741824.0;
+ tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
yb = aa - ya + ab;
@@ -1257,7 +1262,7 @@ public class FastMath {
aa = ya + LN_QUICK_COEF[i][0];
ab = yb + LN_QUICK_COEF[i][1];
/* Split y = a */
- tmp = aa * 1073741824.0;
+ tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
yb = aa - ya + ab;
}
@@ -1266,7 +1271,7 @@ public class FastMath {
aa = ya * xa;
ab = ya * xb + yb * xa + yb * xb;
/* split, so now y = a */
- tmp = aa * 1073741824.0;
+ tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
yb = aa - ya + ab;
@@ -1293,7 +1298,7 @@ public class FastMath {
if (hiPrec != null) {
/* split epsilon -> x */
- double tmp = epsilon * 1073741824.0;
+ double tmp = epsilon * HEX_40000000;
double aa = epsilon + tmp - tmp;
double ab = epsilon - aa;
double xa = aa;
@@ -1314,7 +1319,7 @@ public class FastMath {
aa = ya * xa;
ab = ya * xb + yb * xa + yb * xb;
/* split, so now y = a */
- tmp = aa * 1073741824.0;
+ tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
yb = aa - ya + ab;
@@ -1322,7 +1327,7 @@ public class FastMath {
aa = ya + LN_HI_PREC_COEF[i][0];
ab = yb + LN_HI_PREC_COEF[i][1];
/* Split y = a */
- tmp = aa * 1073741824.0;
+ tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
yb = aa - ya + ab;
}
@@ -1454,7 +1459,7 @@ public class FastMath {
return lores;
}
- final double tmp = hiPrec[0] * 1073741824.0;
+ final double tmp = hiPrec[0] * HEX_40000000;
final double lna = hiPrec[0] + tmp - tmp;
final double lnb = hiPrec[0] - lna + hiPrec[1];
@@ -1590,13 +1595,13 @@ public class FastMath {
double ya;
double yb;
if (y < 8e298 && y > -8e298) {
- double tmp1 = y * 1073741824.0;
+ double tmp1 = y * HEX_40000000;
ya = y + tmp1 - tmp1;
yb = y - ya;
} else {
double tmp1 = y * 9.31322574615478515625E-10;
double tmp2 = tmp1 * 9.31322574615478515625E-10;
- ya = (tmp1 + tmp2 - tmp1) * 1073741824.0 * 1073741824.0;
+ ya = (tmp1 + tmp2 - tmp1) * HEX_40000000 * HEX_40000000;
yb = y - ya;
}
@@ -1610,7 +1615,7 @@ public class FastMath {
double lnb = lns[1];
/* resplit lns */
- double tmp1 = lna * 1073741824.0;
+ double tmp1 = lna * HEX_40000000;
double tmp2 = lna + tmp1 - tmp1;
lnb += lna - tmp2;
lna = tmp2;
@@ -1941,7 +1946,7 @@ public class FastMath {
final double cosEpsB = polyCosine(epsilon);
// Split epsilon xa + xb = x
- final double temp = sinEpsA * 1073741824.0;
+ final double temp = sinEpsA * HEX_40000000;
double temp2 = (sinEpsA + temp) - temp;
sinEpsB += sinEpsA - temp2;
sinEpsA = temp2;
@@ -2085,7 +2090,7 @@ public class FastMath {
final double cosEpsB = polyCosine(epsilon);
// Split epsilon xa + xb = x
- double temp = sinEpsA * 1073741824.0;
+ double temp = sinEpsA * HEX_40000000;
double temp2 = (sinEpsA + temp) - temp;
sinEpsB += sinEpsA - temp2;
sinEpsA = temp2;
@@ -2177,11 +2182,11 @@ public class FastMath {
double est = sina/cosa;
/* Split the estimate to get more accurate read on division rounding */
- temp = est * 1073741824.0;
+ temp = est * HEX_40000000;
double esta = (est + temp) - temp;
double estb = est - esta;
- temp = cosa * 1073741824.0;
+ temp = cosa * HEX_40000000;
double cosaa = (cosa + temp) - temp;
double cosab = cosa - cosaa;
@@ -2765,7 +2770,7 @@ public class FastMath {
epsA = temp;
/* Compute eps = eps / (1.0 + xa*tangent) */
- temp = xa * 1073741824.0;
+ temp = xa * HEX_40000000;
double ya = xa + temp - temp;
double yb = xb + xa - ya;
xa = ya;
@@ -2791,11 +2796,11 @@ public class FastMath {
zb += xb * TANGENT_TABLE_B[idx];
ya = epsA / za;
- temp = ya * 1073741824.0;
+ temp = ya * HEX_40000000;
final double yaa = (ya + temp) - temp;
final double yab = ya - yaa;
- temp = za * 1073741824.0;
+ temp = za * HEX_40000000;
final double zaa = (za + temp) - temp;
final double zab = za - zaa;
@@ -2974,13 +2979,13 @@ public class FastMath {
}
// Split y
- double temp = x * 1073741824.0;
+ double temp = x * HEX_40000000;
final double xa = x + temp - temp;
final double xb = x - xa;
// Compute ratio r = x/y
final double r = y/x;
- temp = r * 1073741824.0;
+ temp = r * HEX_40000000;
double ra = r + temp - temp;
double rb = r - ra;
@@ -3024,7 +3029,7 @@ public class FastMath {
/* Compute asin(x) = atan(x/sqrt(1-x*x)) */
/* Split x */
- double temp = x * 1073741824.0;
+ double temp = x * HEX_40000000;
final double xa = x + temp - temp;
final double xb = x - xa;
@@ -3046,7 +3051,7 @@ public class FastMath {
/* Square root */
double y;
y = sqrt(za);
- temp = y * 1073741824.0;
+ temp = y * HEX_40000000;
ya = y + temp - temp;
yb = y - ya;
@@ -3058,7 +3063,7 @@ public class FastMath {
// Compute ratio r = x/y
double r = x/y;
- temp = r * 1073741824.0;
+ temp = r * HEX_40000000;
double ra = r + temp - temp;
double rb = r - ra;
@@ -3100,7 +3105,7 @@ public class FastMath {
/* Compute acos(x) = atan(sqrt(1-x*x)/x) */
/* Split x */
- double temp = x * 1073741824.0;
+ double temp = x * HEX_40000000;
final double xa = x + temp - temp;
final double xb = x - xa;
@@ -3121,7 +3126,7 @@ public class FastMath {
/* Square root */
double y = sqrt(za);
- temp = y * 1073741824.0;
+ temp = y * HEX_40000000;
ya = y + temp - temp;
yb = y - ya;
@@ -3141,8 +3146,8 @@ public class FastMath {
return Math.PI/2; // so return the appropriate value
}
- if (abs(r) < Double.MAX_VALUE/1073741824.0){ // is it safe to split r ?
- temp = r * 1073741824.0;
+ if (abs(r) < Double.MAX_VALUE/HEX_40000000){ // is it safe to split r ?
+ temp = r * HEX_40000000;
} else {
temp = 0.0;
}
@@ -3214,13 +3219,13 @@ public class FastMath {
est += (xs - est*est*est) / (3*est*est);
// Do one round of Newton's method in extended precision to get the last bit right.
- double temp = est * 1073741824.0;
+ double temp = est * HEX_40000000;
double ya = est + temp - temp;
double yb = est - ya;
double za = ya * ya;
double zb = ya * yb * 2.0 + yb * yb;
- temp = za * 1073741824.0;
+ temp = za * HEX_40000000;
double temp2 = za + temp - temp;
zb += za - temp2;
za = temp2;
@@ -3255,12 +3260,13 @@ public class FastMath {
return x;
}
+ // These are PI/180 split into high and low order bits
final double facta = 0.01745329052209854;
final double factb = 1.997844754509471E-9;
double temp = 0;
- if (abs(x) < Double.MAX_VALUE/1073741824.0) { // prevent overflow to infinity
- temp = x * 1073741824.0;
+ if (abs(x) < Double.MAX_VALUE/HEX_40000000) { // prevent overflow to infinity
+ temp = x * HEX_40000000;
}
double xa = x + temp - temp;
double xb = x - xa;
@@ -3283,12 +3289,13 @@ public class FastMath {
return x;
}
+ // These are 180/PI split into high and low order bits
final double facta = 57.2957763671875;
final double factb = 3.145894820876798E-6;
double temp = 0;
- if (abs(x) < Double.MAX_VALUE/1073741824.0) { // prevent overflow to infinity
- temp = x * 1073741824.0;
+ if (abs(x) < Double.MAX_VALUE/HEX_40000000) { // prevent overflow to infinity
+ temp = x * HEX_40000000;
}
double xa = x + temp - temp;
double xb = x - xa;