added fast cubic root computation

JIRA: MATH-375

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/branches/MATH_2_X@996166 13f79535-47bb-0310-9956-ffa450edef68

src/main/java/org/apache/commons/math/util/FastMath.java

@@ -154,6 +154,13 @@ public class FastMath {
      */
     private static final double EIGHTHES[] = {0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625};
 
+    /* Table of 2^((n+2)/3) */
+    private static final double CBRTTWO[] = { 0.6299605249474366,
+                                            0.7937005259840998, 
+                                            1.0, 
+                                            1.2599210498948732, 
+                                            1.5874010519681994 };
+
     // Initialize tables
     static {
         int i;
@@ -212,14 +219,6 @@ public class FastMath {
         return Math.sqrt(a);
     }
 
-    /** Compute the cubic root of a number.
-     * @param a number on which evaluation is done
-     * @return cubic root of a
-     */
-    public static double cbrt(final double a) {
-        return Math.cbrt(a);
-    }
-
     /** Compute the hyperbolic cosine of a number.
      * @param a number on which evaluation is done
      * @return hyperbolic cosine of a
@@ -1020,11 +1019,6 @@ public class FastMath {
                 ya = aa + tmp - tmp;
                 yb = aa - ya + ab;
 
-                if (hiPrec != null) {
-                    hiPrec[0] = ya;
-                    hiPrec[1] = yb;
-                }
-
                 return ya + yb;
             }
         }
@@ -2886,6 +2880,90 @@ public class FastMath {
       return atan(ra, rb, x<0);
     }
 
+    /** Compute the cubic root of a number.
+     * @param a number on which evaluation is done
+     * @return cubic root of a
+     */
+    public static double cbrt(double x) {
+      /* Convert input double to bits */
+      long inbits = Double.doubleToLongBits(x);
+      int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
+      boolean subnormal = false;
+
+      if (exponent == -1023) {
+          if (x == 0) {
+              return x;
+          }
+
+          /* Subnormal, so normalize */
+          subnormal = true;
+          x *= 1.8014398509481984E16;  // 2^54
+          inbits = Double.doubleToLongBits(x);
+          exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
+      }
+
+      if (exponent == 1024) {
+          // Nan or infinity.  Don't care which.
+          return x;
+      }
+
+      /* Divide the exponent by 3 */
+      int exp3 = exponent / 3;
+
+      /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
+      double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) | 
+                                          (long)(((exp3 + 1023) & 0x7ff)) << 52);
+
+      /* This will be a number between 1 and 2 */
+      final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
+      
+      /* Estimate the cube root of mant by polynomial */
+      double est = -0.010714690733195933;
+      est = est * mant + 0.0875862700108075;
+      est = est * mant + -0.3058015757857271;
+      est = est * mant + 0.7249995199969751;
+      est = est * mant + 0.5039018405998233;
+
+      est *= CBRTTWO[exponent % 3 + 2];
+
+      // est should now be good to about 15 bits of precision.   Do 2 rounds of 
+      // Newton's method to get closer,  this should get us full double precision
+      // Scale down x for the purpose of doing newtons method.  This avoids over/under flows.
+      final double xs = x / (p2*p2*p2); 
+      est += (xs - est*est*est) / (3*est*est);
+      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 ya = est + temp - temp;
+      double yb = est - ya;
+
+      double za = ya * ya;
+      double zb = ya * yb * 2.0 + yb * yb;
+      temp = za * 1073741824.0;
+      double temp2 = za + temp - temp;
+      zb += (za - temp2);
+      za = temp2;
+
+      zb = za * yb + ya * zb + zb * yb;
+      za = za * ya;
+
+      double na = xs - za;
+      double nb = -(na - xs + za);
+      nb -= zb;
+
+      est += (na+nb)/(3*est*est);
+
+      /* Scale by a power of two, so this is exact. */
+      est *= p2;
+
+      if (subnormal) {
+          est *= 3.814697265625E-6;  // 2^-18
+      }
+
+      return est;
+    }
+
     /**
      *  Convert degrees to radians, with error of less than 0.5 ULP
      *  @param x angle in degrees

src/test/java/org/apache/commons/math/util/FastMathTest.java

@@ -787,6 +787,28 @@ public class FastMathTest {
         Assert.assertTrue("acos() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp < MAX_ERROR_ULP);
     }
 
+    @Test
+    public void testCbrtAccuracy() {
+        double maxerrulp = 0.0;
+
+        for (int i=0; i<10000; i++) {
+            double x = ((generator.nextDouble() * 200.0) - 100.0) * generator.nextDouble(); 
+
+            double tst = FastMath.cbrt(x);
+            double ref = cbrt(field.newDfp(x)).toDouble();
+            double err = (tst - ref) / ref;
+
+            if (err != 0) {
+                double ulp = Math.abs(ref - Double.longBitsToDouble((Double.doubleToLongBits(ref) ^ 1)));
+                double errulp = field.newDfp(tst).subtract(cbrt(field.newDfp(x))).divide(field.newDfp(ulp)).toDouble(); 
+                //System.out.println(x+"\t"+tst+"\t"+ref+"\t"+err+"\t"+errulp); 
+                maxerrulp = Math.max(maxerrulp, Math.abs(errulp));
+            }
+        }
+
+        Assert.assertTrue("cbrt() had errors in excess of " + MAX_ERROR_ULP + " ULP", maxerrulp < MAX_ERROR_ULP);
+    }
+
     private Dfp cbrt(Dfp x) {
       boolean negative=false;
 
@@ -938,6 +960,7 @@ public class FastMathTest {
             time = System.currentTimeMillis() - time;
             System.out.println("FastMath.cos " + time + "\t" + x);
 
+            x = 0;
             time = System.currentTimeMillis();
             for (int i = 0; i < numberOfRuns; i++)
                 x += StrictMath.acos(i / 10000000.0);
@@ -980,6 +1003,20 @@ public class FastMathTest {
             System.out.println("FastMath.atan " + time + "\t" + x);
 
             x = 0;
+            time = System.currentTimeMillis();
+            for (int i = 0; i < numberOfRuns; i++)
+                x += StrictMath.cbrt(i / 1000000.0);
+            time = System.currentTimeMillis() - time;
+            System.out.print("StrictMath.cbrt " + time + "\t" + x + "\t");
+
+            x = 0;
+            time = System.currentTimeMillis();
+            for (int i = 0; i < numberOfRuns; i++)
+                x += FastMath.cbrt(i / 1000000.0);
+            time = System.currentTimeMillis() - time;
+            System.out.println("FastMath.cbrt " + time + "\t" + x);
+
+           x = 0;
             time = System.currentTimeMillis();
             for (int i = 0; i < numberOfRuns; i++)
                 x += StrictMath.expm1(-i / 100000.0);