1 /* 2 * Copyright (c) 1998, 2001, Oracle and/or its affiliates. All rights reserved. 3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. 4 * 5 * This code is free software; you can redistribute it and/or modify it 6 * under the terms of the GNU General Public License version 2 only, as 7 * published by the Free Software Foundation. Oracle designates this 8 * particular file as subject to the "Classpath" exception as provided 9 * by Oracle in the LICENSE file that accompanied this code. 10 * 11 * This code is distributed in the hope that it will be useful, but WITHOUT 12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 14 * version 2 for more details (a copy is included in the LICENSE file that 15 * accompanied this code). 16 * 17 * You should have received a copy of the GNU General Public License version 18 * 2 along with this work; if not, write to the Free Software Foundation, 19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 20 * 21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA 22 * or visit www.oracle.com if you need additional information or have any 23 * questions. 24 */ 25 26 /* __ieee754_exp(x) 27 * Returns the exponential of x. 28 * 29 * Method 30 * 1. Argument reduction: 31 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. 32 * Given x, find r and integer k such that 33 * 34 * x = k*ln2 + r, |r| <= 0.5*ln2. 35 * 36 * Here r will be represented as r = hi-lo for better 37 * accuracy. 38 * 39 * 2. Approximation of exp(r) by a special rational function on 40 * the interval [0,0.34658]: 41 * Write 42 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... 43 * We use a special Reme algorithm on [0,0.34658] to generate 44 * a polynomial of degree 5 to approximate R. The maximum error 45 * of this polynomial approximation is bounded by 2**-59. In 46 * other words, 68 * Special cases: 69 * exp(INF) is INF, exp(NaN) is NaN; 70 * exp(-INF) is 0, and 71 * for finite argument, only exp(0)=1 is exact. 72 * 73 * Accuracy: 74 * according to an error analysis, the error is always less than 75 * 1 ulp (unit in the last place). 76 * 77 * Misc. info. 78 * For IEEE double 79 * if x > 7.09782712893383973096e+02 then exp(x) overflow 80 * if x < -7.45133219101941108420e+02 then exp(x) underflow 81 * 82 * Constants: 83 * The hexadecimal values are the intended ones for the following 84 * constants. The decimal values may be used, provided that the 85 * compiler will convert from decimal to binary accurately enough 86 * to produce the hexadecimal values shown. 87 */ 88 89 #include "fdlibm.h" 90 91 #ifdef __STDC__ 92 static const double 93 #else 94 static double 95 #endif 96 one = 1.0, 97 halF[2] = {0.5,-0.5,}, 98 huge = 1.0e+300, 99 twom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/ 100 o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ 101 u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ 102 ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ 103 -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */ 104 ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ 105 -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */ 106 invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ 107 P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ 108 P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ 109 P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ 110 P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ 111 P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ 112 113 114 #ifdef __STDC__ 115 double __ieee754_exp(double x) /* default IEEE double exp */ 116 #else 117 double __ieee754_exp(x) /* default IEEE double exp */ 118 double x; 119 #endif 120 { 121 double y,hi=0,lo=0,c,t; 122 int k=0,xsb; 123 unsigned hx; 124 125 hx = __HI(x); /* high word of x */ 126 xsb = (hx>>31)&1; /* sign bit of x */ 127 hx &= 0x7fffffff; /* high word of |x| */ 128 129 /* filter out non-finite argument */ 130 if(hx >= 0x40862E42) { /* if |x|>=709.78... */ 131 if(hx>=0x7ff00000) { 132 if(((hx&0xfffff)|__LO(x))!=0) 133 return x+x; /* NaN */ 134 else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ 135 } 136 if(x > o_threshold) return huge*huge; /* overflow */ 137 if(x < u_threshold) return twom1000*twom1000; /* underflow */ 138 } 139 140 /* argument reduction */ 141 if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ 142 if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ 143 hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; 144 } else { 145 k = invln2*x+halF[xsb]; 146 t = k; 147 hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */ 148 lo = t*ln2LO[0]; 149 } 150 x = hi - lo; 151 } 152 else if(hx < 0x3e300000) { /* when |x|<2**-28 */ 153 if(huge+x>one) return one+x;/* trigger inexact */ 154 } 155 else k = 0; 156 157 /* x is now in primary range */ 158 t = x*x; 159 c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); 160 if(k==0) return one-((x*c)/(c-2.0)-x); 161 else y = one-((lo-(x*c)/(2.0-c))-hi); 162 if(k >= -1021) { 163 __HI(y) += (k<<20); /* add k to y's exponent */ 164 return y; 165 } else { 166 __HI(y) += ((k+1000)<<20);/* add k to y's exponent */ 167 return y*twom1000; 168 } 169 } | 1 /* 2 * Copyright (c) 1998, 2016, Oracle and/or its affiliates. All rights reserved. 3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. 4 * 5 * This code is free software; you can redistribute it and/or modify it 6 * under the terms of the GNU General Public License version 2 only, as 7 * published by the Free Software Foundation. Oracle designates this 8 * particular file as subject to the "Classpath" exception as provided 9 * by Oracle in the LICENSE file that accompanied this code. 10 * 11 * This code is distributed in the hope that it will be useful, but WITHOUT 12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 14 * version 2 for more details (a copy is included in the LICENSE file that 15 * accompanied this code). 16 * 17 * You should have received a copy of the GNU General Public License version 18 * 2 along with this work; if not, write to the Free Software Foundation, 19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 20 * 21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA 22 * or visit www.oracle.com if you need additional information or have any 23 * questions. 24 */ 25 26 /** 27 * A transliteration of the "Freely Distributable Math Library" 28 * algorithms from C into Java. That is, this port of the algorithms 29 * is as close to the C originals as possible while still being 30 * readable legal Java. 31 */ 32 public class FdlibmTranslit { 33 private FdlibmTranslit() { 34 throw new UnsupportedOperationException("No FdLibmTranslit instances for you."); 35 } 36 37 /** 38 * Return the low-order 32 bits of the double argument as an int. 39 */ 40 private static int __LO(double x) { 41 long transducer = Double.doubleToRawLongBits(x); 42 return (int)transducer; 43 } 44 45 /** 46 * Return a double with its low-order bits of the second argument 47 * and the high-order bits of the first argument.. 48 */ 49 private static double __LO(double x, int low) { 50 long transX = Double.doubleToRawLongBits(x); 51 return Double.longBitsToDouble((transX & 0xFFFF_FFFF_0000_0000L)|low ); 52 } 53 54 /** 55 * Return the high-order 32 bits of the double argument as an int. 56 */ 57 private static int __HI(double x) { 58 long transducer = Double.doubleToRawLongBits(x); 59 return (int)(transducer >> 32); 60 } 61 62 /** 63 * Return a double with its high-order bits of the second argument 64 * and the low-order bits of the first argument.. 65 */ 66 private static double __HI(double x, int high) { 67 long transX = Double.doubleToRawLongBits(x); 68 return Double.longBitsToDouble((transX & 0x0000_0000_FFFF_FFFFL) | 69 ( ((long)high)) << 32 ); 70 } 71 72 public static double hypot(double x, double y) { 73 return Hypot.compute(x, y); 74 } 75 76 /** 77 * cbrt(x) 78 * Return cube root of x 79 */ 80 public static class Cbrt { 81 // unsigned 82 private static final int B1 = 715094163; /* B1 = (682-0.03306235651)*2**20 */ 83 private static final int B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ 84 85 private static final double C = 5.42857142857142815906e-01; /* 19/35 = 0x3FE15F15, 0xF15F15F1 */ 86 private static final double D = -7.05306122448979611050e-01; /* -864/1225 = 0xBFE691DE, 0x2532C834 */ 87 private static final double E = 1.41428571428571436819e+00; /* 99/70 = 0x3FF6A0EA, 0x0EA0EA0F */ 88 private static final double F = 1.60714285714285720630e+00; /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */ 89 private static final double G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */ 90 91 public static strictfp double compute(double x) { 92 int hx; 93 double r, s, t=0.0, w; 94 int sign; // unsigned 95 96 hx = __HI(x); // high word of x 97 sign = hx & 0x80000000; // sign= sign(x) 98 hx ^= sign; 99 if (hx >= 0x7ff00000) 100 return (x+x); // cbrt(NaN,INF) is itself 101 if ((hx | __LO(x)) == 0) 102 return(x); // cbrt(0) is itself 103 104 x = __HI(x, hx); // x <- |x| 105 // rough cbrt to 5 bits 106 if (hx < 0x00100000) { // subnormal number 107 t = __HI(t, 0x43500000); // set t= 2**54 108 t *= x; 109 t = __HI(t, __HI(t)/3+B2); 110 } else { 111 t = __HI(t, hx/3+B1); 112 } 113 114 // new cbrt to 23 bits, may be implemented in single precision 115 r = t * t/x; 116 s = C + r*t; 117 t *= G + F/(s + E + D/s); 118 119 // chopped to 20 bits and make it larger than cbrt(x) 120 t = __LO(t, 0); 121 t = __HI(t, __HI(t)+0x00000001); 122 123 124 // one step newton iteration to 53 bits with error less than 0.667 ulps 125 s = t * t; // t*t is exact 126 r = x / s; 127 w = t + t; 128 r= (r - t)/(w + r); // r-s is exact 129 t= t + t*r; 130 131 // retore the sign bit 132 t = __HI(t, __HI(t) | sign); 133 return(t); 134 } 135 } 136 137 /** 138 * hypot(x,y) 139 * 140 * Method : 141 * If (assume round-to-nearest) z = x*x + y*y 142 * has error less than sqrt(2)/2 ulp, than 143 * sqrt(z) has error less than 1 ulp (exercise). 144 * 145 * So, compute sqrt(x*x + y*y) with some care as 146 * follows to get the error below 1 ulp: 147 * 148 * Assume x > y > 0; 149 * (if possible, set rounding to round-to-nearest) 150 * 1. if x > 2y use 151 * x1*x1 + (y*y + (x2*(x + x1))) for x*x + y*y 152 * where x1 = x with lower 32 bits cleared, x2 = x - x1; else 153 * 2. if x <= 2y use 154 * t1*y1 + ((x-y) * (x-y) + (t1*y2 + t2*y)) 155 * where t1 = 2x with lower 32 bits cleared, t2 = 2x - t1, 156 * y1= y with lower 32 bits chopped, y2 = y - y1. 157 * 158 * NOTE: scaling may be necessary if some argument is too 159 * large or too tiny 160 * 161 * Special cases: 162 * hypot(x,y) is INF if x or y is +INF or -INF; else 163 * hypot(x,y) is NAN if x or y is NAN. 164 * 165 * Accuracy: 166 * hypot(x,y) returns sqrt(x^2 + y^2) with error less 167 * than 1 ulps (units in the last place) 168 */ 169 static class Hypot { 170 public static double compute(double x, double y) { 171 double a = x; 172 double b = y; 173 double t1, t2, y1, y2, w; 174 int j, k, ha, hb; 175 176 ha = __HI(x) & 0x7fffffff; // high word of x 177 hb = __HI(y) & 0x7fffffff; // high word of y 178 if(hb > ha) { 179 a = y; 180 b = x; 181 j = ha; 182 ha = hb; 183 hb = j; 184 } else { 185 a = x; 186 b = y; 187 } 188 a = __HI(a, ha); // a <- |a| 189 b = __HI(b, hb); // b <- |b| 190 if ((ha - hb) > 0x3c00000) { 191 return a + b; // x / y > 2**60 192 } 193 k=0; 194 if (ha > 0x5f300000) { // a>2**500 195 if (ha >= 0x7ff00000) { // Inf or NaN 196 w = a + b; // for sNaN 197 if (((ha & 0xfffff) | __LO(a)) == 0) 198 w = a; 199 if (((hb ^ 0x7ff00000) | __LO(b)) == 0) 200 w = b; 201 return w; 202 } 203 // scale a and b by 2**-600 204 ha -= 0x25800000; 205 hb -= 0x25800000; 206 k += 600; 207 a = __HI(a, ha); 208 b = __HI(b, hb); 209 } 210 if (hb < 0x20b00000) { // b < 2**-500 211 if (hb <= 0x000fffff) { // subnormal b or 0 */ 212 if ((hb | (__LO(b))) == 0) 213 return a; 214 t1 = 0; 215 t1 = __HI(t1, 0x7fd00000); // t1=2^1022 216 b *= t1; 217 a *= t1; 218 k -= 1022; 219 } else { // scale a and b by 2^600 220 ha += 0x25800000; // a *= 2^600 221 hb += 0x25800000; // b *= 2^600 222 k -= 600; 223 a = __HI(a, ha); 224 b = __HI(b, hb); 225 } 226 } 227 // medium size a and b 228 w = a - b; 229 if (w > b) { 230 t1 = 0; 231 t1 = __HI(t1, ha); 232 t2 = a - t1; 233 w = Math.sqrt(t1*t1 - (b*(-b) - t2 * (a + t1))); 234 } else { 235 a = a + a; 236 y1 = 0; 237 y1 = __HI(y1, hb); 238 y2 = b - y1; 239 t1 = 0; 240 t1 = __HI(t1, ha + 0x00100000); 241 t2 = a - t1; 242 w = Math.sqrt(t1*y1 - (w*(-w) - (t1*y2 + t2*b))); 243 } 244 if (k != 0) { 245 t1 = 1.0; 246 int t1_hi = __HI(t1); 247 t1_hi += (k << 20); 248 t1 = __HI(t1, t1_hi); 249 return t1 * w; 250 } else 251 return w; 252 } 253 } 254 255 /** 256 * Returns the exponential of x. 257 * 258 * Method 259 * 1. Argument reduction: 260 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. 261 * Given x, find r and integer k such that 262 * 263 * x = k*ln2 + r, |r| <= 0.5*ln2. 264 * 265 * Here r will be represented as r = hi-lo for better 266 * accuracy. 267 * 268 * 2. Approximation of exp(r) by a special rational function on 269 * the interval [0,0.34658]: 270 * Write 271 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... 272 * We use a special Reme algorithm on [0,0.34658] to generate 273 * a polynomial of degree 5 to approximate R. The maximum error 274 * of this polynomial approximation is bounded by 2**-59. In 275 * other words, 297 * Special cases: 298 * exp(INF) is INF, exp(NaN) is NaN; 299 * exp(-INF) is 0, and 300 * for finite argument, only exp(0)=1 is exact. 301 * 302 * Accuracy: 303 * according to an error analysis, the error is always less than 304 * 1 ulp (unit in the last place). 305 * 306 * Misc. info. 307 * For IEEE double 308 * if x > 7.09782712893383973096e+02 then exp(x) overflow 309 * if x < -7.45133219101941108420e+02 then exp(x) underflow 310 * 311 * Constants: 312 * The hexadecimal values are the intended ones for the following 313 * constants. The decimal values may be used, provided that the 314 * compiler will convert from decimal to binary accurately enough 315 * to produce the hexadecimal values shown. 316 */ 317 static class Exp { 318 private static final double one = 1.0; 319 private static final double[] halF = {0.5,-0.5,}; 320 private static final double huge = 1.0e+300; 321 private static final double twom1000= 9.33263618503218878990e-302; /* 2**-1000=0x01700000,0*/ 322 private static final double o_threshold= 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */ 323 private static final double u_threshold= -7.45133219101941108420e+02; /* 0xc0874910, 0xD52D3051 */ 324 private static final double[] ln2HI ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ 325 -6.93147180369123816490e-01}; /* 0xbfe62e42, 0xfee00000 */ 326 private static final double[] ln2LO ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ 327 -1.90821492927058770002e-10,}; /* 0xbdea39ef, 0x35793c76 */ 328 private static final double invln2 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */ 329 private static final double P1 = 1.66666666666666019037e-01; /* 0x3FC55555, 0x5555553E */ 330 private static final double P2 = -2.77777777770155933842e-03; /* 0xBF66C16C, 0x16BEBD93 */ 331 private static final double P3 = 6.61375632143793436117e-05; /* 0x3F11566A, 0xAF25DE2C */ 332 private static final double P4 = -1.65339022054652515390e-06; /* 0xBEBBBD41, 0xC5D26BF1 */ 333 private static final double P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ 334 335 public static strictfp double compute(double x) { 336 double y,hi=0,lo=0,c,t; 337 int k=0,xsb; 338 /*unsigned*/ int hx; 339 340 hx = __HI(x); /* high word of x */ 341 xsb = (hx>>31)&1; /* sign bit of x */ 342 hx &= 0x7fffffff; /* high word of |x| */ 343 344 /* filter out non-finite argument */ 345 if(hx >= 0x40862E42) { /* if |x|>=709.78... */ 346 if(hx>=0x7ff00000) { 347 if(((hx&0xfffff)|__LO(x))!=0) 348 return x+x; /* NaN */ 349 else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ 350 } 351 if(x > o_threshold) return huge*huge; /* overflow */ 352 if(x < u_threshold) return twom1000*twom1000; /* underflow */ 353 } 354 355 /* argument reduction */ 356 if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ 357 if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ 358 hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; 359 } else { 360 k = (int)(invln2*x+halF[xsb]); 361 t = k; 362 hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */ 363 lo = t*ln2LO[0]; 364 } 365 x = hi - lo; 366 } 367 else if(hx < 0x3e300000) { /* when |x|<2**-28 */ 368 if(huge+x>one) return one+x;/* trigger inexact */ 369 } 370 else k = 0; 371 372 /* x is now in primary range */ 373 t = x*x; 374 c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); 375 if(k==0) return one-((x*c)/(c-2.0)-x); 376 else y = one-((lo-(x*c)/(2.0-c))-hi); 377 if(k >= -1021) { 378 y = __HI(y, __HI(y) + (k<<20)); /* add k to y's exponent */ 379 return y; 380 } else { 381 y = __HI(y, __HI(y) + ((k+1000)<<20));/* add k to y's exponent */ 382 return y*twom1000; 383 } 384 } 385 } 386 } |