Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

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.

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1371670 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2012-08-10 11:54:46 +00:00
parent f465d0aa69
commit 746c998e1a
2 changed files with 52 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
src/main/java/org/apache/commons/math3/util/FastMath.java
View File

@ -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;
// 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;
e = e >> 1;
// 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;
}
return result;
// 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 resultHigh + resultLow;
}
/**

17
src/test/java/org/apache/commons/math3/util/FastMathTest.java
View File

@ -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);
}
}
}