Author: luc
Date: Fri Aug 10 11:54:46 2012
New Revision: 1371670
URL: http://svn.apache.org/viewvc?rev=1371670&view=rev
Log:
Fixed accuracy issues in FastMath.pow(double, int).
The fixed version is slightly slower, but still much faster than
FastMath.pow(double, double). Some random testing showed that the
accuracy is now always better than 0.5ulp, even for large exponent.
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/FastMath.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/util/FastMathTest.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/FastMath.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/FastMath.java?rev=1371670&r1=1371669&r2=1371670&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/FastMath.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/FastMath.java
Fri Aug 10 11:54:46 2012
@@ -1589,6 +1589,7 @@ public class FastMath {
* @return d<sup>e</sup>
*/
public static double pow(double d, int e) {
+
if (e == 0) {
return 1.0;
} else if (e < 0) {
@@ -1596,17 +1597,54 @@ public class FastMath {
d = 1.0 / d;
}
- double result = 1;
- double d2p = d;
+ // split d as two 26 bits numbers
+ // beware the following expressions must NOT be simplified, they rely
on floating point arithmetic properties
+ final int splitFactor = 0x8000001;
+ final double cd = splitFactor * d;
+ final double d1High = cd - (cd - d);
+ final double d1Low = d - d1High;
+
+ // prepare result
+ double resultHigh = 1;
+ double resultLow = 0;
+
+ // d^(2p)
+ double d2p = d;
+ double d2pHigh = d1High;
+ double d2pLow = d1Low;
+
while (e != 0) {
+
if ((e & 0x1) != 0) {
- result *= d2p;
- }
- d2p *= d2p;
+ // accurate multiplication result = result * d^(2p) using
Veltkamp TwoProduct algorithm
+ // beware the following expressions must NOT be simplified,
they rely on floating point arithmetic properties
+ final double tmpHigh = resultHigh * d2p;
+ final double cRH = splitFactor * resultHigh;
+ final double rHH = cRH - (cRH - resultHigh);
+ final double rHL = resultHigh - rHH;
+ final double tmpLow = rHL * d2pLow - (((tmpHigh - rHH *
d2pHigh) - rHL * d2pHigh) - rHH * d2pLow);
+ resultHigh = tmpHigh;
+ resultLow = resultLow * d2p + tmpLow;
+ }
+
+ // accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp
TwoProduct algorithm
+ // beware the following expressions must NOT be simplified, they
rely on floating point arithmetic properties
+ final double tmpHigh = d2pHigh * d2p;
+ final double cD2pH = splitFactor * d2pHigh;
+ final double d2pHH = cD2pH - (cD2pH - d2pHigh);
+ final double d2pHL = d2pHigh - d2pHH;
+ final double tmpLow = d2pHL * d2pLow - (((tmpHigh - d2pHH *
d2pHigh) - d2pHL * d2pHigh) - d2pHH * d2pLow);
+ final double cTmpH = splitFactor * tmpHigh;
+ d2pHigh = cTmpH - (cTmpH - tmpHigh);
+ d2pLow = d2pLow * d2p + tmpLow + (tmpHigh - d2pHigh);
+ d2p = d2pHigh + d2pLow;
+
e = e >> 1;
+
}
- return result;
+ return resultHigh + resultLow;
+
}
/**
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/util/FastMathTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/util/FastMathTest.java?rev=1371670&r1=1371669&r2=1371670&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/util/FastMathTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/util/FastMathTest.java
Fri Aug 10 11:54:46 2012
@@ -20,12 +20,12 @@ import java.lang.reflect.Method;
import java.lang.reflect.Modifier;
import java.lang.reflect.Type;
+import org.apache.commons.math3.TestUtils;
import org.apache.commons.math3.dfp.Dfp;
import org.apache.commons.math3.dfp.DfpField;
import org.apache.commons.math3.dfp.DfpMath;
import org.apache.commons.math3.random.MersenneTwister;
import org.apache.commons.math3.random.RandomGenerator;
-import org.apache.commons.math3.TestUtils;
import org.junit.Assert;
import org.junit.Before;
import org.junit.Ignore;
@@ -1108,15 +1108,16 @@ public class FastMathTest {
@Test
public void testIntPow() {
- final double base = 1.23456789;
final int maxExp = 300;
-
+ DfpField field = new DfpField(40);
+ final double base = 1.23456789;
+ Dfp baseDfp = field.newDfp(base);
+ Dfp dfpPower = field.getOne();
for (int i = 0; i < maxExp; i++) {
- final double expected = FastMath.pow(base, (double) i);
- Assert.assertEquals("exp=" + i,
- expected,
- FastMath.pow(base, i),
- 60 * Math.ulp(expected));
+ Assert.assertEquals("exp=" + i, dfpPower.toDouble(),
FastMath.pow(base, i),
+ 0.6 * FastMath.ulp(dfpPower.toDouble()));
+ dfpPower = dfpPower.multiply(baseDfp);
}
}
+
}
