< prev index next >

test/java/lang/StrictMath/FdlibmTranslit.java

Print this page

        

@@ -71,10 +71,71 @@
     public static double hypot(double x, double y) {
         return Hypot.compute(x, y);
     }
 
     /**
+     * cbrt(x)
+     * Return cube root of x
+     */
+    public static class Cbrt {
+        // unsigned
+        private static final int B1 = 715094163; /* B1 = (682-0.03306235651)*2**20 */
+        private static final int B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
+
+        private static final double C =  5.42857142857142815906e-01; /* 19/35     = 0x3FE15F15, 0xF15F15F1 */
+        private static final double D = -7.05306122448979611050e-01; /* -864/1225 = 0xBFE691DE, 0x2532C834 */
+        private static final double E =  1.41428571428571436819e+00; /* 99/70     = 0x3FF6A0EA, 0x0EA0EA0F */
+        private static final double F =  1.60714285714285720630e+00; /* 45/28     = 0x3FF9B6DB, 0x6DB6DB6E */
+        private static final double G =  3.57142857142857150787e-01; /* 5/14      = 0x3FD6DB6D, 0xB6DB6DB7 */
+
+        public static strictfp double compute(double x) {
+            int     hx;
+            double  r, s, t=0.0, w;
+            int sign; // unsigned
+
+            hx = __HI(x);           // high word of x
+            sign = hx & 0x80000000;             // sign= sign(x)
+            hx  ^= sign;
+            if (hx >= 0x7ff00000)
+                return (x+x); // cbrt(NaN,INF) is itself
+            if ((hx | __LO(x)) == 0)
+                return(x);          // cbrt(0) is itself
+
+            x = __HI(x, hx);   // x <- |x|
+            // rough cbrt to 5 bits
+            if (hx < 0x00100000) {               // subnormal number
+                t = __HI(t, 0x43500000);          // set t= 2**54
+                t *= x;
+                t = __HI(t, __HI(t)/3+B2);
+            } else {
+                t = __HI(t, hx/3+B1);
+            }
+
+            // new cbrt to 23 bits, may be implemented in single precision
+            r = t * t/x;
+            s = C + r*t;
+            t *= G + F/(s + E + D/s);
+
+            // chopped to 20 bits and make it larger than cbrt(x)
+            t = __LO(t, 0);
+            t = __HI(t, __HI(t)+0x00000001);
+
+
+            // one step newton iteration to 53 bits with error less than 0.667 ulps
+            s = t * t;          // t*t is exact
+            r = x / s;
+            w = t + t;
+            r= (r - t)/(w + r);  // r-s is exact
+            t= t + t*r;
+
+            // retore the sign bit
+            t = __HI(t, __HI(t) | sign);
+            return(t);
+        }
+    }
+
+    /**
      * hypot(x,y)
      *
      * Method :
      *      If (assume round-to-nearest) z = x*x + y*y
      *      has error less than sqrt(2)/2 ulp, than
< prev index next >