< prev index next >
test/java/lang/StrictMath/FdlibmTranslit.java
Print this page
*** 71,80 ****
--- 71,141 ----
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 >