1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/src/share/vm/runtime/sharedRuntimeTrans.cpp Sat Dec 01 00:00:00 2007 +0000 1.3 @@ -0,0 +1,719 @@ 1.4 +/* 1.5 + * Copyright 2005 Sun Microsystems, Inc. All Rights Reserved. 1.6 + * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. 1.7 + * 1.8 + * This code is free software; you can redistribute it and/or modify it 1.9 + * under the terms of the GNU General Public License version 2 only, as 1.10 + * published by the Free Software Foundation. 1.11 + * 1.12 + * This code is distributed in the hope that it will be useful, but WITHOUT 1.13 + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 1.14 + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 1.15 + * version 2 for more details (a copy is included in the LICENSE file that 1.16 + * accompanied this code). 1.17 + * 1.18 + * You should have received a copy of the GNU General Public License version 1.19 + * 2 along with this work; if not, write to the Free Software Foundation, 1.20 + * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 1.21 + * 1.22 + * Please contact Sun Microsystems, Inc., 4150 Network Circle, Santa Clara, 1.23 + * CA 95054 USA or visit www.sun.com if you need additional information or 1.24 + * have any questions. 1.25 + * 1.26 + */ 1.27 + 1.28 +#include "incls/_precompiled.incl" 1.29 +#include "incls/_sharedRuntimeTrans.cpp.incl" 1.30 + 1.31 +// This file contains copies of the fdlibm routines used by 1.32 +// StrictMath. It turns out that it is almost always required to use 1.33 +// these runtime routines; the Intel CPU doesn't meet the Java 1.34 +// specification for sin/cos outside a certain limited argument range, 1.35 +// and the SPARC CPU doesn't appear to have sin/cos instructions. It 1.36 +// also turns out that avoiding the indirect call through function 1.37 +// pointer out to libjava.so in SharedRuntime speeds these routines up 1.38 +// by roughly 15% on both Win32/x86 and Solaris/SPARC. 1.39 + 1.40 +// Enabling optimizations in this file causes incorrect code to be 1.41 +// generated; can not figure out how to turn down optimization for one 1.42 +// file in the IDE on Windows 1.43 +#ifdef WIN32 1.44 +# pragma optimize ( "", off ) 1.45 +#endif 1.46 + 1.47 +#include <math.h> 1.48 + 1.49 +// VM_LITTLE_ENDIAN is #defined appropriately in the Makefiles 1.50 +// [jk] this is not 100% correct because the float word order may different 1.51 +// from the byte order (e.g. on ARM) 1.52 +#ifdef VM_LITTLE_ENDIAN 1.53 +# define __HI(x) *(1+(int*)&x) 1.54 +# define __LO(x) *(int*)&x 1.55 +#else 1.56 +# define __HI(x) *(int*)&x 1.57 +# define __LO(x) *(1+(int*)&x) 1.58 +#endif 1.59 + 1.60 +double copysign(double x, double y) { 1.61 + __HI(x) = (__HI(x)&0x7fffffff)|(__HI(y)&0x80000000); 1.62 + return x; 1.63 +} 1.64 + 1.65 +/* 1.66 + * ==================================================== 1.67 + * Copyright (C) 1998 by Sun Microsystems, Inc. All rights reserved. 1.68 + * 1.69 + * Developed at SunSoft, a Sun Microsystems, Inc. business. 1.70 + * Permission to use, copy, modify, and distribute this 1.71 + * software is freely granted, provided that this notice 1.72 + * is preserved. 1.73 + * ==================================================== 1.74 + */ 1.75 + 1.76 +/* 1.77 + * scalbn (double x, int n) 1.78 + * scalbn(x,n) returns x* 2**n computed by exponent 1.79 + * manipulation rather than by actually performing an 1.80 + * exponentiation or a multiplication. 1.81 + */ 1.82 + 1.83 +static const double 1.84 +two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ 1.85 + twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ 1.86 + hugeX = 1.0e+300, 1.87 + tiny = 1.0e-300; 1.88 + 1.89 +double scalbn (double x, int n) { 1.90 + int k,hx,lx; 1.91 + hx = __HI(x); 1.92 + lx = __LO(x); 1.93 + k = (hx&0x7ff00000)>>20; /* extract exponent */ 1.94 + if (k==0) { /* 0 or subnormal x */ 1.95 + if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ 1.96 + x *= two54; 1.97 + hx = __HI(x); 1.98 + k = ((hx&0x7ff00000)>>20) - 54; 1.99 + if (n< -50000) return tiny*x; /*underflow*/ 1.100 + } 1.101 + if (k==0x7ff) return x+x; /* NaN or Inf */ 1.102 + k = k+n; 1.103 + if (k > 0x7fe) return hugeX*copysign(hugeX,x); /* overflow */ 1.104 + if (k > 0) /* normal result */ 1.105 + {__HI(x) = (hx&0x800fffff)|(k<<20); return x;} 1.106 + if (k <= -54) { 1.107 + if (n > 50000) /* in case integer overflow in n+k */ 1.108 + return hugeX*copysign(hugeX,x); /*overflow*/ 1.109 + else return tiny*copysign(tiny,x); /*underflow*/ 1.110 + } 1.111 + k += 54; /* subnormal result */ 1.112 + __HI(x) = (hx&0x800fffff)|(k<<20); 1.113 + return x*twom54; 1.114 +} 1.115 + 1.116 +/* __ieee754_log(x) 1.117 + * Return the logrithm of x 1.118 + * 1.119 + * Method : 1.120 + * 1. Argument Reduction: find k and f such that 1.121 + * x = 2^k * (1+f), 1.122 + * where sqrt(2)/2 < 1+f < sqrt(2) . 1.123 + * 1.124 + * 2. Approximation of log(1+f). 1.125 + * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) 1.126 + * = 2s + 2/3 s**3 + 2/5 s**5 + ....., 1.127 + * = 2s + s*R 1.128 + * We use a special Reme algorithm on [0,0.1716] to generate 1.129 + * a polynomial of degree 14 to approximate R The maximum error 1.130 + * of this polynomial approximation is bounded by 2**-58.45. In 1.131 + * other words, 1.132 + * 2 4 6 8 10 12 14 1.133 + * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s 1.134 + * (the values of Lg1 to Lg7 are listed in the program) 1.135 + * and 1.136 + * | 2 14 | -58.45 1.137 + * | Lg1*s +...+Lg7*s - R(z) | <= 2 1.138 + * | | 1.139 + * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. 1.140 + * In order to guarantee error in log below 1ulp, we compute log 1.141 + * by 1.142 + * log(1+f) = f - s*(f - R) (if f is not too large) 1.143 + * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) 1.144 + * 1.145 + * 3. Finally, log(x) = k*ln2 + log(1+f). 1.146 + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) 1.147 + * Here ln2 is split into two floating point number: 1.148 + * ln2_hi + ln2_lo, 1.149 + * where n*ln2_hi is always exact for |n| < 2000. 1.150 + * 1.151 + * Special cases: 1.152 + * log(x) is NaN with signal if x < 0 (including -INF) ; 1.153 + * log(+INF) is +INF; log(0) is -INF with signal; 1.154 + * log(NaN) is that NaN with no signal. 1.155 + * 1.156 + * Accuracy: 1.157 + * according to an error analysis, the error is always less than 1.158 + * 1 ulp (unit in the last place). 1.159 + * 1.160 + * Constants: 1.161 + * The hexadecimal values are the intended ones for the following 1.162 + * constants. The decimal values may be used, provided that the 1.163 + * compiler will convert from decimal to binary accurately enough 1.164 + * to produce the hexadecimal values shown. 1.165 + */ 1.166 + 1.167 +static const double 1.168 +ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ 1.169 + ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ 1.170 + Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ 1.171 + Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ 1.172 + Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ 1.173 + Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ 1.174 + Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ 1.175 + Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 1.176 + Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ 1.177 + 1.178 +static double zero = 0.0; 1.179 + 1.180 +static double __ieee754_log(double x) { 1.181 + double hfsq,f,s,z,R,w,t1,t2,dk; 1.182 + int k,hx,i,j; 1.183 + unsigned lx; 1.184 + 1.185 + hx = __HI(x); /* high word of x */ 1.186 + lx = __LO(x); /* low word of x */ 1.187 + 1.188 + k=0; 1.189 + if (hx < 0x00100000) { /* x < 2**-1022 */ 1.190 + if (((hx&0x7fffffff)|lx)==0) 1.191 + return -two54/zero; /* log(+-0)=-inf */ 1.192 + if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ 1.193 + k -= 54; x *= two54; /* subnormal number, scale up x */ 1.194 + hx = __HI(x); /* high word of x */ 1.195 + } 1.196 + if (hx >= 0x7ff00000) return x+x; 1.197 + k += (hx>>20)-1023; 1.198 + hx &= 0x000fffff; 1.199 + i = (hx+0x95f64)&0x100000; 1.200 + __HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */ 1.201 + k += (i>>20); 1.202 + f = x-1.0; 1.203 + if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */ 1.204 + if(f==zero) { 1.205 + if (k==0) return zero; 1.206 + else {dk=(double)k; return dk*ln2_hi+dk*ln2_lo;} 1.207 + } 1.208 + R = f*f*(0.5-0.33333333333333333*f); 1.209 + if(k==0) return f-R; else {dk=(double)k; 1.210 + return dk*ln2_hi-((R-dk*ln2_lo)-f);} 1.211 + } 1.212 + s = f/(2.0+f); 1.213 + dk = (double)k; 1.214 + z = s*s; 1.215 + i = hx-0x6147a; 1.216 + w = z*z; 1.217 + j = 0x6b851-hx; 1.218 + t1= w*(Lg2+w*(Lg4+w*Lg6)); 1.219 + t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 1.220 + i |= j; 1.221 + R = t2+t1; 1.222 + if(i>0) { 1.223 + hfsq=0.5*f*f; 1.224 + if(k==0) return f-(hfsq-s*(hfsq+R)); else 1.225 + return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); 1.226 + } else { 1.227 + if(k==0) return f-s*(f-R); else 1.228 + return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); 1.229 + } 1.230 +} 1.231 + 1.232 +JRT_LEAF(jdouble, SharedRuntime::dlog(jdouble x)) 1.233 + return __ieee754_log(x); 1.234 +JRT_END 1.235 + 1.236 +/* __ieee754_log10(x) 1.237 + * Return the base 10 logarithm of x 1.238 + * 1.239 + * Method : 1.240 + * Let log10_2hi = leading 40 bits of log10(2) and 1.241 + * log10_2lo = log10(2) - log10_2hi, 1.242 + * ivln10 = 1/log(10) rounded. 1.243 + * Then 1.244 + * n = ilogb(x), 1.245 + * if(n<0) n = n+1; 1.246 + * x = scalbn(x,-n); 1.247 + * log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x)) 1.248 + * 1.249 + * Note 1: 1.250 + * To guarantee log10(10**n)=n, where 10**n is normal, the rounding 1.251 + * mode must set to Round-to-Nearest. 1.252 + * Note 2: 1.253 + * [1/log(10)] rounded to 53 bits has error .198 ulps; 1.254 + * log10 is monotonic at all binary break points. 1.255 + * 1.256 + * Special cases: 1.257 + * log10(x) is NaN with signal if x < 0; 1.258 + * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; 1.259 + * log10(NaN) is that NaN with no signal; 1.260 + * log10(10**N) = N for N=0,1,...,22. 1.261 + * 1.262 + * Constants: 1.263 + * The hexadecimal values are the intended ones for the following constants. 1.264 + * The decimal values may be used, provided that the compiler will convert 1.265 + * from decimal to binary accurately enough to produce the hexadecimal values 1.266 + * shown. 1.267 + */ 1.268 + 1.269 +static const double 1.270 +ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */ 1.271 + log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ 1.272 + log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ 1.273 + 1.274 +static double __ieee754_log10(double x) { 1.275 + double y,z; 1.276 + int i,k,hx; 1.277 + unsigned lx; 1.278 + 1.279 + hx = __HI(x); /* high word of x */ 1.280 + lx = __LO(x); /* low word of x */ 1.281 + 1.282 + k=0; 1.283 + if (hx < 0x00100000) { /* x < 2**-1022 */ 1.284 + if (((hx&0x7fffffff)|lx)==0) 1.285 + return -two54/zero; /* log(+-0)=-inf */ 1.286 + if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ 1.287 + k -= 54; x *= two54; /* subnormal number, scale up x */ 1.288 + hx = __HI(x); /* high word of x */ 1.289 + } 1.290 + if (hx >= 0x7ff00000) return x+x; 1.291 + k += (hx>>20)-1023; 1.292 + i = ((unsigned)k&0x80000000)>>31; 1.293 + hx = (hx&0x000fffff)|((0x3ff-i)<<20); 1.294 + y = (double)(k+i); 1.295 + __HI(x) = hx; 1.296 + z = y*log10_2lo + ivln10*__ieee754_log(x); 1.297 + return z+y*log10_2hi; 1.298 +} 1.299 + 1.300 +JRT_LEAF(jdouble, SharedRuntime::dlog10(jdouble x)) 1.301 + return __ieee754_log10(x); 1.302 +JRT_END 1.303 + 1.304 + 1.305 +/* __ieee754_exp(x) 1.306 + * Returns the exponential of x. 1.307 + * 1.308 + * Method 1.309 + * 1. Argument reduction: 1.310 + * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. 1.311 + * Given x, find r and integer k such that 1.312 + * 1.313 + * x = k*ln2 + r, |r| <= 0.5*ln2. 1.314 + * 1.315 + * Here r will be represented as r = hi-lo for better 1.316 + * accuracy. 1.317 + * 1.318 + * 2. Approximation of exp(r) by a special rational function on 1.319 + * the interval [0,0.34658]: 1.320 + * Write 1.321 + * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... 1.322 + * We use a special Reme algorithm on [0,0.34658] to generate 1.323 + * a polynomial of degree 5 to approximate R. The maximum error 1.324 + * of this polynomial approximation is bounded by 2**-59. In 1.325 + * other words, 1.326 + * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 1.327 + * (where z=r*r, and the values of P1 to P5 are listed below) 1.328 + * and 1.329 + * | 5 | -59 1.330 + * | 2.0+P1*z+...+P5*z - R(z) | <= 2 1.331 + * | | 1.332 + * The computation of exp(r) thus becomes 1.333 + * 2*r 1.334 + * exp(r) = 1 + ------- 1.335 + * R - r 1.336 + * r*R1(r) 1.337 + * = 1 + r + ----------- (for better accuracy) 1.338 + * 2 - R1(r) 1.339 + * where 1.340 + * 2 4 10 1.341 + * R1(r) = r - (P1*r + P2*r + ... + P5*r ). 1.342 + * 1.343 + * 3. Scale back to obtain exp(x): 1.344 + * From step 1, we have 1.345 + * exp(x) = 2^k * exp(r) 1.346 + * 1.347 + * Special cases: 1.348 + * exp(INF) is INF, exp(NaN) is NaN; 1.349 + * exp(-INF) is 0, and 1.350 + * for finite argument, only exp(0)=1 is exact. 1.351 + * 1.352 + * Accuracy: 1.353 + * according to an error analysis, the error is always less than 1.354 + * 1 ulp (unit in the last place). 1.355 + * 1.356 + * Misc. info. 1.357 + * For IEEE double 1.358 + * if x > 7.09782712893383973096e+02 then exp(x) overflow 1.359 + * if x < -7.45133219101941108420e+02 then exp(x) underflow 1.360 + * 1.361 + * Constants: 1.362 + * The hexadecimal values are the intended ones for the following 1.363 + * constants. The decimal values may be used, provided that the 1.364 + * compiler will convert from decimal to binary accurately enough 1.365 + * to produce the hexadecimal values shown. 1.366 + */ 1.367 + 1.368 +static const double 1.369 +one = 1.0, 1.370 + halF[2] = {0.5,-0.5,}, 1.371 + twom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/ 1.372 + o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ 1.373 + u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ 1.374 + ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ 1.375 + -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */ 1.376 + ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ 1.377 + -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */ 1.378 + invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ 1.379 + P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ 1.380 + P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ 1.381 + P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ 1.382 + P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ 1.383 + P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ 1.384 + 1.385 +static double __ieee754_exp(double x) { 1.386 + double y,hi=0,lo=0,c,t; 1.387 + int k=0,xsb; 1.388 + unsigned hx; 1.389 + 1.390 + hx = __HI(x); /* high word of x */ 1.391 + xsb = (hx>>31)&1; /* sign bit of x */ 1.392 + hx &= 0x7fffffff; /* high word of |x| */ 1.393 + 1.394 + /* filter out non-finite argument */ 1.395 + if(hx >= 0x40862E42) { /* if |x|>=709.78... */ 1.396 + if(hx>=0x7ff00000) { 1.397 + if(((hx&0xfffff)|__LO(x))!=0) 1.398 + return x+x; /* NaN */ 1.399 + else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ 1.400 + } 1.401 + if(x > o_threshold) return hugeX*hugeX; /* overflow */ 1.402 + if(x < u_threshold) return twom1000*twom1000; /* underflow */ 1.403 + } 1.404 + 1.405 + /* argument reduction */ 1.406 + if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ 1.407 + if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ 1.408 + hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; 1.409 + } else { 1.410 + k = (int)(invln2*x+halF[xsb]); 1.411 + t = k; 1.412 + hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */ 1.413 + lo = t*ln2LO[0]; 1.414 + } 1.415 + x = hi - lo; 1.416 + } 1.417 + else if(hx < 0x3e300000) { /* when |x|<2**-28 */ 1.418 + if(hugeX+x>one) return one+x;/* trigger inexact */ 1.419 + } 1.420 + else k = 0; 1.421 + 1.422 + /* x is now in primary range */ 1.423 + t = x*x; 1.424 + c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); 1.425 + if(k==0) return one-((x*c)/(c-2.0)-x); 1.426 + else y = one-((lo-(x*c)/(2.0-c))-hi); 1.427 + if(k >= -1021) { 1.428 + __HI(y) += (k<<20); /* add k to y's exponent */ 1.429 + return y; 1.430 + } else { 1.431 + __HI(y) += ((k+1000)<<20);/* add k to y's exponent */ 1.432 + return y*twom1000; 1.433 + } 1.434 +} 1.435 + 1.436 +JRT_LEAF(jdouble, SharedRuntime::dexp(jdouble x)) 1.437 + return __ieee754_exp(x); 1.438 +JRT_END 1.439 + 1.440 +/* __ieee754_pow(x,y) return x**y 1.441 + * 1.442 + * n 1.443 + * Method: Let x = 2 * (1+f) 1.444 + * 1. Compute and return log2(x) in two pieces: 1.445 + * log2(x) = w1 + w2, 1.446 + * where w1 has 53-24 = 29 bit trailing zeros. 1.447 + * 2. Perform y*log2(x) = n+y' by simulating muti-precision 1.448 + * arithmetic, where |y'|<=0.5. 1.449 + * 3. Return x**y = 2**n*exp(y'*log2) 1.450 + * 1.451 + * Special cases: 1.452 + * 1. (anything) ** 0 is 1 1.453 + * 2. (anything) ** 1 is itself 1.454 + * 3. (anything) ** NAN is NAN 1.455 + * 4. NAN ** (anything except 0) is NAN 1.456 + * 5. +-(|x| > 1) ** +INF is +INF 1.457 + * 6. +-(|x| > 1) ** -INF is +0 1.458 + * 7. +-(|x| < 1) ** +INF is +0 1.459 + * 8. +-(|x| < 1) ** -INF is +INF 1.460 + * 9. +-1 ** +-INF is NAN 1.461 + * 10. +0 ** (+anything except 0, NAN) is +0 1.462 + * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 1.463 + * 12. +0 ** (-anything except 0, NAN) is +INF 1.464 + * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF 1.465 + * 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) 1.466 + * 15. +INF ** (+anything except 0,NAN) is +INF 1.467 + * 16. +INF ** (-anything except 0,NAN) is +0 1.468 + * 17. -INF ** (anything) = -0 ** (-anything) 1.469 + * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) 1.470 + * 19. (-anything except 0 and inf) ** (non-integer) is NAN 1.471 + * 1.472 + * Accuracy: 1.473 + * pow(x,y) returns x**y nearly rounded. In particular 1.474 + * pow(integer,integer) 1.475 + * always returns the correct integer provided it is 1.476 + * representable. 1.477 + * 1.478 + * Constants : 1.479 + * The hexadecimal values are the intended ones for the following 1.480 + * constants. The decimal values may be used, provided that the 1.481 + * compiler will convert from decimal to binary accurately enough 1.482 + * to produce the hexadecimal values shown. 1.483 + */ 1.484 + 1.485 +static const double 1.486 +bp[] = {1.0, 1.5,}, 1.487 + dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ 1.488 + dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ 1.489 + zeroX = 0.0, 1.490 + two = 2.0, 1.491 + two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ 1.492 + /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ 1.493 + L1X = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ 1.494 + L2X = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ 1.495 + L3X = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */ 1.496 + L4X = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */ 1.497 + L5X = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */ 1.498 + L6X = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */ 1.499 + lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ 1.500 + lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */ 1.501 + lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */ 1.502 + ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */ 1.503 + cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */ 1.504 + cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */ 1.505 + cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/ 1.506 + ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */ 1.507 + ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/ 1.508 + ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ 1.509 + 1.510 +double __ieee754_pow(double x, double y) { 1.511 + double z,ax,z_h,z_l,p_h,p_l; 1.512 + double y1,t1,t2,r,s,t,u,v,w; 1.513 + int i0,i1,i,j,k,yisint,n; 1.514 + int hx,hy,ix,iy; 1.515 + unsigned lx,ly; 1.516 + 1.517 + i0 = ((*(int*)&one)>>29)^1; i1=1-i0; 1.518 + hx = __HI(x); lx = __LO(x); 1.519 + hy = __HI(y); ly = __LO(y); 1.520 + ix = hx&0x7fffffff; iy = hy&0x7fffffff; 1.521 + 1.522 + /* y==zero: x**0 = 1 */ 1.523 + if((iy|ly)==0) return one; 1.524 + 1.525 + /* +-NaN return x+y */ 1.526 + if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) || 1.527 + iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) 1.528 + return x+y; 1.529 + 1.530 + /* determine if y is an odd int when x < 0 1.531 + * yisint = 0 ... y is not an integer 1.532 + * yisint = 1 ... y is an odd int 1.533 + * yisint = 2 ... y is an even int 1.534 + */ 1.535 + yisint = 0; 1.536 + if(hx<0) { 1.537 + if(iy>=0x43400000) yisint = 2; /* even integer y */ 1.538 + else if(iy>=0x3ff00000) { 1.539 + k = (iy>>20)-0x3ff; /* exponent */ 1.540 + if(k>20) { 1.541 + j = ly>>(52-k); 1.542 + if((unsigned)(j<<(52-k))==ly) yisint = 2-(j&1); 1.543 + } else if(ly==0) { 1.544 + j = iy>>(20-k); 1.545 + if((j<<(20-k))==iy) yisint = 2-(j&1); 1.546 + } 1.547 + } 1.548 + } 1.549 + 1.550 + /* special value of y */ 1.551 + if(ly==0) { 1.552 + if (iy==0x7ff00000) { /* y is +-inf */ 1.553 + if(((ix-0x3ff00000)|lx)==0) 1.554 + return y - y; /* inf**+-1 is NaN */ 1.555 + else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */ 1.556 + return (hy>=0)? y: zeroX; 1.557 + else /* (|x|<1)**-,+inf = inf,0 */ 1.558 + return (hy<0)?-y: zeroX; 1.559 + } 1.560 + if(iy==0x3ff00000) { /* y is +-1 */ 1.561 + if(hy<0) return one/x; else return x; 1.562 + } 1.563 + if(hy==0x40000000) return x*x; /* y is 2 */ 1.564 + if(hy==0x3fe00000) { /* y is 0.5 */ 1.565 + if(hx>=0) /* x >= +0 */ 1.566 + return sqrt(x); 1.567 + } 1.568 + } 1.569 + 1.570 + ax = fabsd(x); 1.571 + /* special value of x */ 1.572 + if(lx==0) { 1.573 + if(ix==0x7ff00000||ix==0||ix==0x3ff00000){ 1.574 + z = ax; /*x is +-0,+-inf,+-1*/ 1.575 + if(hy<0) z = one/z; /* z = (1/|x|) */ 1.576 + if(hx<0) { 1.577 + if(((ix-0x3ff00000)|yisint)==0) { 1.578 + z = (z-z)/(z-z); /* (-1)**non-int is NaN */ 1.579 + } else if(yisint==1) 1.580 + z = -1.0*z; /* (x<0)**odd = -(|x|**odd) */ 1.581 + } 1.582 + return z; 1.583 + } 1.584 + } 1.585 + 1.586 + n = (hx>>31)+1; 1.587 + 1.588 + /* (x<0)**(non-int) is NaN */ 1.589 + if((n|yisint)==0) return (x-x)/(x-x); 1.590 + 1.591 + s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ 1.592 + if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */ 1.593 + 1.594 + /* |y| is huge */ 1.595 + if(iy>0x41e00000) { /* if |y| > 2**31 */ 1.596 + if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */ 1.597 + if(ix<=0x3fefffff) return (hy<0)? hugeX*hugeX:tiny*tiny; 1.598 + if(ix>=0x3ff00000) return (hy>0)? hugeX*hugeX:tiny*tiny; 1.599 + } 1.600 + /* over/underflow if x is not close to one */ 1.601 + if(ix<0x3fefffff) return (hy<0)? s*hugeX*hugeX:s*tiny*tiny; 1.602 + if(ix>0x3ff00000) return (hy>0)? s*hugeX*hugeX:s*tiny*tiny; 1.603 + /* now |1-x| is tiny <= 2**-20, suffice to compute 1.604 + log(x) by x-x^2/2+x^3/3-x^4/4 */ 1.605 + t = ax-one; /* t has 20 trailing zeros */ 1.606 + w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25)); 1.607 + u = ivln2_h*t; /* ivln2_h has 21 sig. bits */ 1.608 + v = t*ivln2_l-w*ivln2; 1.609 + t1 = u+v; 1.610 + __LO(t1) = 0; 1.611 + t2 = v-(t1-u); 1.612 + } else { 1.613 + double ss,s2,s_h,s_l,t_h,t_l; 1.614 + n = 0; 1.615 + /* take care subnormal number */ 1.616 + if(ix<0x00100000) 1.617 + {ax *= two53; n -= 53; ix = __HI(ax); } 1.618 + n += ((ix)>>20)-0x3ff; 1.619 + j = ix&0x000fffff; 1.620 + /* determine interval */ 1.621 + ix = j|0x3ff00000; /* normalize ix */ 1.622 + if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */ 1.623 + else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */ 1.624 + else {k=0;n+=1;ix -= 0x00100000;} 1.625 + __HI(ax) = ix; 1.626 + 1.627 + /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ 1.628 + u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ 1.629 + v = one/(ax+bp[k]); 1.630 + ss = u*v; 1.631 + s_h = ss; 1.632 + __LO(s_h) = 0; 1.633 + /* t_h=ax+bp[k] High */ 1.634 + t_h = zeroX; 1.635 + __HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18); 1.636 + t_l = ax - (t_h-bp[k]); 1.637 + s_l = v*((u-s_h*t_h)-s_h*t_l); 1.638 + /* compute log(ax) */ 1.639 + s2 = ss*ss; 1.640 + r = s2*s2*(L1X+s2*(L2X+s2*(L3X+s2*(L4X+s2*(L5X+s2*L6X))))); 1.641 + r += s_l*(s_h+ss); 1.642 + s2 = s_h*s_h; 1.643 + t_h = 3.0+s2+r; 1.644 + __LO(t_h) = 0; 1.645 + t_l = r-((t_h-3.0)-s2); 1.646 + /* u+v = ss*(1+...) */ 1.647 + u = s_h*t_h; 1.648 + v = s_l*t_h+t_l*ss; 1.649 + /* 2/(3log2)*(ss+...) */ 1.650 + p_h = u+v; 1.651 + __LO(p_h) = 0; 1.652 + p_l = v-(p_h-u); 1.653 + z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ 1.654 + z_l = cp_l*p_h+p_l*cp+dp_l[k]; 1.655 + /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */ 1.656 + t = (double)n; 1.657 + t1 = (((z_h+z_l)+dp_h[k])+t); 1.658 + __LO(t1) = 0; 1.659 + t2 = z_l-(((t1-t)-dp_h[k])-z_h); 1.660 + } 1.661 + 1.662 + /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ 1.663 + y1 = y; 1.664 + __LO(y1) = 0; 1.665 + p_l = (y-y1)*t1+y*t2; 1.666 + p_h = y1*t1; 1.667 + z = p_l+p_h; 1.668 + j = __HI(z); 1.669 + i = __LO(z); 1.670 + if (j>=0x40900000) { /* z >= 1024 */ 1.671 + if(((j-0x40900000)|i)!=0) /* if z > 1024 */ 1.672 + return s*hugeX*hugeX; /* overflow */ 1.673 + else { 1.674 + if(p_l+ovt>z-p_h) return s*hugeX*hugeX; /* overflow */ 1.675 + } 1.676 + } else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */ 1.677 + if(((j-0xc090cc00)|i)!=0) /* z < -1075 */ 1.678 + return s*tiny*tiny; /* underflow */ 1.679 + else { 1.680 + if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ 1.681 + } 1.682 + } 1.683 + /* 1.684 + * compute 2**(p_h+p_l) 1.685 + */ 1.686 + i = j&0x7fffffff; 1.687 + k = (i>>20)-0x3ff; 1.688 + n = 0; 1.689 + if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */ 1.690 + n = j+(0x00100000>>(k+1)); 1.691 + k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */ 1.692 + t = zeroX; 1.693 + __HI(t) = (n&~(0x000fffff>>k)); 1.694 + n = ((n&0x000fffff)|0x00100000)>>(20-k); 1.695 + if(j<0) n = -n; 1.696 + p_h -= t; 1.697 + } 1.698 + t = p_l+p_h; 1.699 + __LO(t) = 0; 1.700 + u = t*lg2_h; 1.701 + v = (p_l-(t-p_h))*lg2+t*lg2_l; 1.702 + z = u+v; 1.703 + w = v-(z-u); 1.704 + t = z*z; 1.705 + t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); 1.706 + r = (z*t1)/(t1-two)-(w+z*w); 1.707 + z = one-(r-z); 1.708 + j = __HI(z); 1.709 + j += (n<<20); 1.710 + if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */ 1.711 + else __HI(z) += (n<<20); 1.712 + return s*z; 1.713 +} 1.714 + 1.715 + 1.716 +JRT_LEAF(jdouble, SharedRuntime::dpow(jdouble x, jdouble y)) 1.717 + return __ieee754_pow(x, y); 1.718 +JRT_END 1.719 + 1.720 +#ifdef WIN32 1.721 +# pragma optimize ( "", on ) 1.722 +#endif