commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1371670 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math3/util/FastMath.java test/java/org/apache/commons/math3/util/FastMathTest.java
Date Fri, 10 Aug 2012 11:54:47 GMT
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);
         }
     }
+
 }



Mime
View raw message