68 return Double.longBitsToDouble((transX & 0x0000_0000_FFFF_FFFFL)|( ((long)high)) << 32 ); |
68 return Double.longBitsToDouble((transX & 0x0000_0000_FFFF_FFFFL)|( ((long)high)) << 32 ); |
69 } |
69 } |
70 |
70 |
71 public static double hypot(double x, double y) { |
71 public static double hypot(double x, double y) { |
72 return Hypot.compute(x, y); |
72 return Hypot.compute(x, y); |
|
73 } |
|
74 |
|
75 /** |
|
76 * cbrt(x) |
|
77 * Return cube root of x |
|
78 */ |
|
79 public static class Cbrt { |
|
80 // unsigned |
|
81 private static final int B1 = 715094163; /* B1 = (682-0.03306235651)*2**20 */ |
|
82 private static final int B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ |
|
83 |
|
84 private static final double C = 5.42857142857142815906e-01; /* 19/35 = 0x3FE15F15, 0xF15F15F1 */ |
|
85 private static final double D = -7.05306122448979611050e-01; /* -864/1225 = 0xBFE691DE, 0x2532C834 */ |
|
86 private static final double E = 1.41428571428571436819e+00; /* 99/70 = 0x3FF6A0EA, 0x0EA0EA0F */ |
|
87 private static final double F = 1.60714285714285720630e+00; /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */ |
|
88 private static final double G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */ |
|
89 |
|
90 public static strictfp double compute(double x) { |
|
91 int hx; |
|
92 double r, s, t=0.0, w; |
|
93 int sign; // unsigned |
|
94 |
|
95 hx = __HI(x); // high word of x |
|
96 sign = hx & 0x80000000; // sign= sign(x) |
|
97 hx ^= sign; |
|
98 if (hx >= 0x7ff00000) |
|
99 return (x+x); // cbrt(NaN,INF) is itself |
|
100 if ((hx | __LO(x)) == 0) |
|
101 return(x); // cbrt(0) is itself |
|
102 |
|
103 x = __HI(x, hx); // x <- |x| |
|
104 // rough cbrt to 5 bits |
|
105 if (hx < 0x00100000) { // subnormal number |
|
106 t = __HI(t, 0x43500000); // set t= 2**54 |
|
107 t *= x; |
|
108 t = __HI(t, __HI(t)/3+B2); |
|
109 } else { |
|
110 t = __HI(t, hx/3+B1); |
|
111 } |
|
112 |
|
113 // new cbrt to 23 bits, may be implemented in single precision |
|
114 r = t * t/x; |
|
115 s = C + r*t; |
|
116 t *= G + F/(s + E + D/s); |
|
117 |
|
118 // chopped to 20 bits and make it larger than cbrt(x) |
|
119 t = __LO(t, 0); |
|
120 t = __HI(t, __HI(t)+0x00000001); |
|
121 |
|
122 |
|
123 // one step newton iteration to 53 bits with error less than 0.667 ulps |
|
124 s = t * t; // t*t is exact |
|
125 r = x / s; |
|
126 w = t + t; |
|
127 r= (r - t)/(w + r); // r-s is exact |
|
128 t= t + t*r; |
|
129 |
|
130 // retore the sign bit |
|
131 t = __HI(t, __HI(t) | sign); |
|
132 return(t); |
|
133 } |
73 } |
134 } |
74 |
135 |
75 /** |
136 /** |
76 * hypot(x,y) |
137 * hypot(x,y) |
77 * |
138 * |