Mercurial > hg > truffle
annotate src/share/vm/runtime/sharedRuntimeTrans.cpp @ 1721:413ad0331a0c
6977924: Changes for 6975078 produce build error with certain gcc versions
Summary: The changes introduced for 6975078 assign badHeapOopVal to the _allocation field in the ResourceObj class. In 32 bit linux builds with certain versions of gcc this assignment will be flagged as an error while compiling allocation.cpp. In 32 bit builds the constant value badHeapOopVal (which is cast to an intptr_t) is negative. The _allocation field is typed as an unsigned intptr_t and gcc catches this as an error.
Reviewed-by: jcoomes, ysr, phh
author | johnc |
---|---|
date | Wed, 18 Aug 2010 10:59:06 -0700 |
parents | 126ea7725993 |
children | f95d63e2154a |
rev | line source |
---|---|
0 | 1 /* |
1552
c18cbe5936b8
6941466: Oracle rebranding changes for Hotspot repositories
trims
parents:
0
diff
changeset
|
2 * Copyright (c) 2005, Oracle and/or its affiliates. All rights reserved. |
0 | 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. | |
8 * | |
9 * This code is distributed in the hope that it will be useful, but WITHOUT | |
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
11 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
12 * version 2 for more details (a copy is included in the LICENSE file that | |
13 * accompanied this code). | |
14 * | |
15 * You should have received a copy of the GNU General Public License version | |
16 * 2 along with this work; if not, write to the Free Software Foundation, | |
17 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. | |
18 * | |
1552
c18cbe5936b8
6941466: Oracle rebranding changes for Hotspot repositories
trims
parents:
0
diff
changeset
|
19 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA |
c18cbe5936b8
6941466: Oracle rebranding changes for Hotspot repositories
trims
parents:
0
diff
changeset
|
20 * or visit www.oracle.com if you need additional information or have any |
c18cbe5936b8
6941466: Oracle rebranding changes for Hotspot repositories
trims
parents:
0
diff
changeset
|
21 * questions. |
0 | 22 * |
23 */ | |
24 | |
25 #include "incls/_precompiled.incl" | |
26 #include "incls/_sharedRuntimeTrans.cpp.incl" | |
27 | |
28 // This file contains copies of the fdlibm routines used by | |
29 // StrictMath. It turns out that it is almost always required to use | |
30 // these runtime routines; the Intel CPU doesn't meet the Java | |
31 // specification for sin/cos outside a certain limited argument range, | |
32 // and the SPARC CPU doesn't appear to have sin/cos instructions. It | |
33 // also turns out that avoiding the indirect call through function | |
34 // pointer out to libjava.so in SharedRuntime speeds these routines up | |
35 // by roughly 15% on both Win32/x86 and Solaris/SPARC. | |
36 | |
37 // Enabling optimizations in this file causes incorrect code to be | |
38 // generated; can not figure out how to turn down optimization for one | |
39 // file in the IDE on Windows | |
40 #ifdef WIN32 | |
41 # pragma optimize ( "", off ) | |
42 #endif | |
43 | |
44 #include <math.h> | |
45 | |
46 // VM_LITTLE_ENDIAN is #defined appropriately in the Makefiles | |
47 // [jk] this is not 100% correct because the float word order may different | |
48 // from the byte order (e.g. on ARM) | |
49 #ifdef VM_LITTLE_ENDIAN | |
50 # define __HI(x) *(1+(int*)&x) | |
51 # define __LO(x) *(int*)&x | |
52 #else | |
53 # define __HI(x) *(int*)&x | |
54 # define __LO(x) *(1+(int*)&x) | |
55 #endif | |
56 | |
57 double copysign(double x, double y) { | |
58 __HI(x) = (__HI(x)&0x7fffffff)|(__HI(y)&0x80000000); | |
59 return x; | |
60 } | |
61 | |
62 /* | |
63 * ==================================================== | |
1552
c18cbe5936b8
6941466: Oracle rebranding changes for Hotspot repositories
trims
parents:
0
diff
changeset
|
64 * Copyright (c) 1998 Oracle and/or its affiliates. All rights reserved. |
0 | 65 * |
66 * Developed at SunSoft, a Sun Microsystems, Inc. business. | |
67 * Permission to use, copy, modify, and distribute this | |
68 * software is freely granted, provided that this notice | |
69 * is preserved. | |
70 * ==================================================== | |
71 */ | |
72 | |
73 /* | |
74 * scalbn (double x, int n) | |
75 * scalbn(x,n) returns x* 2**n computed by exponent | |
76 * manipulation rather than by actually performing an | |
77 * exponentiation or a multiplication. | |
78 */ | |
79 | |
80 static const double | |
81 two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ | |
82 twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ | |
83 hugeX = 1.0e+300, | |
84 tiny = 1.0e-300; | |
85 | |
86 double scalbn (double x, int n) { | |
87 int k,hx,lx; | |
88 hx = __HI(x); | |
89 lx = __LO(x); | |
90 k = (hx&0x7ff00000)>>20; /* extract exponent */ | |
91 if (k==0) { /* 0 or subnormal x */ | |
92 if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ | |
93 x *= two54; | |
94 hx = __HI(x); | |
95 k = ((hx&0x7ff00000)>>20) - 54; | |
96 if (n< -50000) return tiny*x; /*underflow*/ | |
97 } | |
98 if (k==0x7ff) return x+x; /* NaN or Inf */ | |
99 k = k+n; | |
100 if (k > 0x7fe) return hugeX*copysign(hugeX,x); /* overflow */ | |
101 if (k > 0) /* normal result */ | |
102 {__HI(x) = (hx&0x800fffff)|(k<<20); return x;} | |
103 if (k <= -54) { | |
104 if (n > 50000) /* in case integer overflow in n+k */ | |
105 return hugeX*copysign(hugeX,x); /*overflow*/ | |
106 else return tiny*copysign(tiny,x); /*underflow*/ | |
107 } | |
108 k += 54; /* subnormal result */ | |
109 __HI(x) = (hx&0x800fffff)|(k<<20); | |
110 return x*twom54; | |
111 } | |
112 | |
113 /* __ieee754_log(x) | |
114 * Return the logrithm of x | |
115 * | |
116 * Method : | |
117 * 1. Argument Reduction: find k and f such that | |
118 * x = 2^k * (1+f), | |
119 * where sqrt(2)/2 < 1+f < sqrt(2) . | |
120 * | |
121 * 2. Approximation of log(1+f). | |
122 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) | |
123 * = 2s + 2/3 s**3 + 2/5 s**5 + ....., | |
124 * = 2s + s*R | |
125 * We use a special Reme algorithm on [0,0.1716] to generate | |
126 * a polynomial of degree 14 to approximate R The maximum error | |
127 * of this polynomial approximation is bounded by 2**-58.45. In | |
128 * other words, | |
129 * 2 4 6 8 10 12 14 | |
130 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s | |
131 * (the values of Lg1 to Lg7 are listed in the program) | |
132 * and | |
133 * | 2 14 | -58.45 | |
134 * | Lg1*s +...+Lg7*s - R(z) | <= 2 | |
135 * | | | |
136 * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. | |
137 * In order to guarantee error in log below 1ulp, we compute log | |
138 * by | |
139 * log(1+f) = f - s*(f - R) (if f is not too large) | |
140 * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) | |
141 * | |
142 * 3. Finally, log(x) = k*ln2 + log(1+f). | |
143 * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) | |
144 * Here ln2 is split into two floating point number: | |
145 * ln2_hi + ln2_lo, | |
146 * where n*ln2_hi is always exact for |n| < 2000. | |
147 * | |
148 * Special cases: | |
149 * log(x) is NaN with signal if x < 0 (including -INF) ; | |
150 * log(+INF) is +INF; log(0) is -INF with signal; | |
151 * log(NaN) is that NaN with no signal. | |
152 * | |
153 * Accuracy: | |
154 * according to an error analysis, the error is always less than | |
155 * 1 ulp (unit in the last place). | |
156 * | |
157 * Constants: | |
158 * The hexadecimal values are the intended ones for the following | |
159 * constants. The decimal values may be used, provided that the | |
160 * compiler will convert from decimal to binary accurately enough | |
161 * to produce the hexadecimal values shown. | |
162 */ | |
163 | |
164 static const double | |
165 ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ | |
166 ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ | |
167 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ | |
168 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ | |
169 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ | |
170 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ | |
171 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ | |
172 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ | |
173 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ | |
174 | |
175 static double zero = 0.0; | |
176 | |
177 static double __ieee754_log(double x) { | |
178 double hfsq,f,s,z,R,w,t1,t2,dk; | |
179 int k,hx,i,j; | |
180 unsigned lx; | |
181 | |
182 hx = __HI(x); /* high word of x */ | |
183 lx = __LO(x); /* low word of x */ | |
184 | |
185 k=0; | |
186 if (hx < 0x00100000) { /* x < 2**-1022 */ | |
187 if (((hx&0x7fffffff)|lx)==0) | |
188 return -two54/zero; /* log(+-0)=-inf */ | |
189 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ | |
190 k -= 54; x *= two54; /* subnormal number, scale up x */ | |
191 hx = __HI(x); /* high word of x */ | |
192 } | |
193 if (hx >= 0x7ff00000) return x+x; | |
194 k += (hx>>20)-1023; | |
195 hx &= 0x000fffff; | |
196 i = (hx+0x95f64)&0x100000; | |
197 __HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */ | |
198 k += (i>>20); | |
199 f = x-1.0; | |
200 if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */ | |
201 if(f==zero) { | |
202 if (k==0) return zero; | |
203 else {dk=(double)k; return dk*ln2_hi+dk*ln2_lo;} | |
204 } | |
205 R = f*f*(0.5-0.33333333333333333*f); | |
206 if(k==0) return f-R; else {dk=(double)k; | |
207 return dk*ln2_hi-((R-dk*ln2_lo)-f);} | |
208 } | |
209 s = f/(2.0+f); | |
210 dk = (double)k; | |
211 z = s*s; | |
212 i = hx-0x6147a; | |
213 w = z*z; | |
214 j = 0x6b851-hx; | |
215 t1= w*(Lg2+w*(Lg4+w*Lg6)); | |
216 t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); | |
217 i |= j; | |
218 R = t2+t1; | |
219 if(i>0) { | |
220 hfsq=0.5*f*f; | |
221 if(k==0) return f-(hfsq-s*(hfsq+R)); else | |
222 return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); | |
223 } else { | |
224 if(k==0) return f-s*(f-R); else | |
225 return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f); | |
226 } | |
227 } | |
228 | |
229 JRT_LEAF(jdouble, SharedRuntime::dlog(jdouble x)) | |
230 return __ieee754_log(x); | |
231 JRT_END | |
232 | |
233 /* __ieee754_log10(x) | |
234 * Return the base 10 logarithm of x | |
235 * | |
236 * Method : | |
237 * Let log10_2hi = leading 40 bits of log10(2) and | |
238 * log10_2lo = log10(2) - log10_2hi, | |
239 * ivln10 = 1/log(10) rounded. | |
240 * Then | |
241 * n = ilogb(x), | |
242 * if(n<0) n = n+1; | |
243 * x = scalbn(x,-n); | |
244 * log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x)) | |
245 * | |
246 * Note 1: | |
247 * To guarantee log10(10**n)=n, where 10**n is normal, the rounding | |
248 * mode must set to Round-to-Nearest. | |
249 * Note 2: | |
250 * [1/log(10)] rounded to 53 bits has error .198 ulps; | |
251 * log10 is monotonic at all binary break points. | |
252 * | |
253 * Special cases: | |
254 * log10(x) is NaN with signal if x < 0; | |
255 * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; | |
256 * log10(NaN) is that NaN with no signal; | |
257 * log10(10**N) = N for N=0,1,...,22. | |
258 * | |
259 * Constants: | |
260 * The hexadecimal values are the intended ones for the following constants. | |
261 * The decimal values may be used, provided that the compiler will convert | |
262 * from decimal to binary accurately enough to produce the hexadecimal values | |
263 * shown. | |
264 */ | |
265 | |
266 static const double | |
267 ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */ | |
268 log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ | |
269 log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ | |
270 | |
271 static double __ieee754_log10(double x) { | |
272 double y,z; | |
273 int i,k,hx; | |
274 unsigned lx; | |
275 | |
276 hx = __HI(x); /* high word of x */ | |
277 lx = __LO(x); /* low word of x */ | |
278 | |
279 k=0; | |
280 if (hx < 0x00100000) { /* x < 2**-1022 */ | |
281 if (((hx&0x7fffffff)|lx)==0) | |
282 return -two54/zero; /* log(+-0)=-inf */ | |
283 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ | |
284 k -= 54; x *= two54; /* subnormal number, scale up x */ | |
285 hx = __HI(x); /* high word of x */ | |
286 } | |
287 if (hx >= 0x7ff00000) return x+x; | |
288 k += (hx>>20)-1023; | |
289 i = ((unsigned)k&0x80000000)>>31; | |
290 hx = (hx&0x000fffff)|((0x3ff-i)<<20); | |
291 y = (double)(k+i); | |
292 __HI(x) = hx; | |
293 z = y*log10_2lo + ivln10*__ieee754_log(x); | |
294 return z+y*log10_2hi; | |
295 } | |
296 | |
297 JRT_LEAF(jdouble, SharedRuntime::dlog10(jdouble x)) | |
298 return __ieee754_log10(x); | |
299 JRT_END | |
300 | |
301 | |
302 /* __ieee754_exp(x) | |
303 * Returns the exponential of x. | |
304 * | |
305 * Method | |
306 * 1. Argument reduction: | |
307 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. | |
308 * Given x, find r and integer k such that | |
309 * | |
310 * x = k*ln2 + r, |r| <= 0.5*ln2. | |
311 * | |
312 * Here r will be represented as r = hi-lo for better | |
313 * accuracy. | |
314 * | |
315 * 2. Approximation of exp(r) by a special rational function on | |
316 * the interval [0,0.34658]: | |
317 * Write | |
318 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... | |
319 * We use a special Reme algorithm on [0,0.34658] to generate | |
320 * a polynomial of degree 5 to approximate R. The maximum error | |
321 * of this polynomial approximation is bounded by 2**-59. In | |
322 * other words, | |
323 * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 | |
324 * (where z=r*r, and the values of P1 to P5 are listed below) | |
325 * and | |
326 * | 5 | -59 | |
327 * | 2.0+P1*z+...+P5*z - R(z) | <= 2 | |
328 * | | | |
329 * The computation of exp(r) thus becomes | |
330 * 2*r | |
331 * exp(r) = 1 + ------- | |
332 * R - r | |
333 * r*R1(r) | |
334 * = 1 + r + ----------- (for better accuracy) | |
335 * 2 - R1(r) | |
336 * where | |
337 * 2 4 10 | |
338 * R1(r) = r - (P1*r + P2*r + ... + P5*r ). | |
339 * | |
340 * 3. Scale back to obtain exp(x): | |
341 * From step 1, we have | |
342 * exp(x) = 2^k * exp(r) | |
343 * | |
344 * Special cases: | |
345 * exp(INF) is INF, exp(NaN) is NaN; | |
346 * exp(-INF) is 0, and | |
347 * for finite argument, only exp(0)=1 is exact. | |
348 * | |
349 * Accuracy: | |
350 * according to an error analysis, the error is always less than | |
351 * 1 ulp (unit in the last place). | |
352 * | |
353 * Misc. info. | |
354 * For IEEE double | |
355 * if x > 7.09782712893383973096e+02 then exp(x) overflow | |
356 * if x < -7.45133219101941108420e+02 then exp(x) underflow | |
357 * | |
358 * Constants: | |
359 * The hexadecimal values are the intended ones for the following | |
360 * constants. The decimal values may be used, provided that the | |
361 * compiler will convert from decimal to binary accurately enough | |
362 * to produce the hexadecimal values shown. | |
363 */ | |
364 | |
365 static const double | |
366 one = 1.0, | |
367 halF[2] = {0.5,-0.5,}, | |
368 twom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/ | |
369 o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ | |
370 u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ | |
371 ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ | |
372 -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */ | |
373 ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ | |
374 -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */ | |
375 invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ | |
376 P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ | |
377 P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ | |
378 P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ | |
379 P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ | |
380 P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ | |
381 | |
382 static double __ieee754_exp(double x) { | |
383 double y,hi=0,lo=0,c,t; | |
384 int k=0,xsb; | |
385 unsigned hx; | |
386 | |
387 hx = __HI(x); /* high word of x */ | |
388 xsb = (hx>>31)&1; /* sign bit of x */ | |
389 hx &= 0x7fffffff; /* high word of |x| */ | |
390 | |
391 /* filter out non-finite argument */ | |
392 if(hx >= 0x40862E42) { /* if |x|>=709.78... */ | |
393 if(hx>=0x7ff00000) { | |
394 if(((hx&0xfffff)|__LO(x))!=0) | |
395 return x+x; /* NaN */ | |
396 else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ | |
397 } | |
398 if(x > o_threshold) return hugeX*hugeX; /* overflow */ | |
399 if(x < u_threshold) return twom1000*twom1000; /* underflow */ | |
400 } | |
401 | |
402 /* argument reduction */ | |
403 if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ | |
404 if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ | |
405 hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; | |
406 } else { | |
407 k = (int)(invln2*x+halF[xsb]); | |
408 t = k; | |
409 hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */ | |
410 lo = t*ln2LO[0]; | |
411 } | |
412 x = hi - lo; | |
413 } | |
414 else if(hx < 0x3e300000) { /* when |x|<2**-28 */ | |
415 if(hugeX+x>one) return one+x;/* trigger inexact */ | |
416 } | |
417 else k = 0; | |
418 | |
419 /* x is now in primary range */ | |
420 t = x*x; | |
421 c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); | |
422 if(k==0) return one-((x*c)/(c-2.0)-x); | |
423 else y = one-((lo-(x*c)/(2.0-c))-hi); | |
424 if(k >= -1021) { | |
425 __HI(y) += (k<<20); /* add k to y's exponent */ | |
426 return y; | |
427 } else { | |
428 __HI(y) += ((k+1000)<<20);/* add k to y's exponent */ | |
429 return y*twom1000; | |
430 } | |
431 } | |
432 | |
433 JRT_LEAF(jdouble, SharedRuntime::dexp(jdouble x)) | |
434 return __ieee754_exp(x); | |
435 JRT_END | |
436 | |
437 /* __ieee754_pow(x,y) return x**y | |
438 * | |
439 * n | |
440 * Method: Let x = 2 * (1+f) | |
441 * 1. Compute and return log2(x) in two pieces: | |
442 * log2(x) = w1 + w2, | |
443 * where w1 has 53-24 = 29 bit trailing zeros. | |
444 * 2. Perform y*log2(x) = n+y' by simulating muti-precision | |
445 * arithmetic, where |y'|<=0.5. | |
446 * 3. Return x**y = 2**n*exp(y'*log2) | |
447 * | |
448 * Special cases: | |
449 * 1. (anything) ** 0 is 1 | |
450 * 2. (anything) ** 1 is itself | |
451 * 3. (anything) ** NAN is NAN | |
452 * 4. NAN ** (anything except 0) is NAN | |
453 * 5. +-(|x| > 1) ** +INF is +INF | |
454 * 6. +-(|x| > 1) ** -INF is +0 | |
455 * 7. +-(|x| < 1) ** +INF is +0 | |
456 * 8. +-(|x| < 1) ** -INF is +INF | |
457 * 9. +-1 ** +-INF is NAN | |
458 * 10. +0 ** (+anything except 0, NAN) is +0 | |
459 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 | |
460 * 12. +0 ** (-anything except 0, NAN) is +INF | |
461 * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF | |
462 * 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) | |
463 * 15. +INF ** (+anything except 0,NAN) is +INF | |
464 * 16. +INF ** (-anything except 0,NAN) is +0 | |
465 * 17. -INF ** (anything) = -0 ** (-anything) | |
466 * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) | |
467 * 19. (-anything except 0 and inf) ** (non-integer) is NAN | |
468 * | |
469 * Accuracy: | |
470 * pow(x,y) returns x**y nearly rounded. In particular | |
471 * pow(integer,integer) | |
472 * always returns the correct integer provided it is | |
473 * representable. | |
474 * | |
475 * Constants : | |
476 * The hexadecimal values are the intended ones for the following | |
477 * constants. The decimal values may be used, provided that the | |
478 * compiler will convert from decimal to binary accurately enough | |
479 * to produce the hexadecimal values shown. | |
480 */ | |
481 | |
482 static const double | |
483 bp[] = {1.0, 1.5,}, | |
484 dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ | |
485 dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ | |
486 zeroX = 0.0, | |
487 two = 2.0, | |
488 two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ | |
489 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ | |
490 L1X = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ | |
491 L2X = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ | |
492 L3X = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */ | |
493 L4X = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */ | |
494 L5X = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */ | |
495 L6X = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */ | |
496 lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ | |
497 lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */ | |
498 lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */ | |
499 ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */ | |
500 cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */ | |
501 cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */ | |
502 cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/ | |
503 ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */ | |
504 ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/ | |
505 ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ | |
506 | |
507 double __ieee754_pow(double x, double y) { | |
508 double z,ax,z_h,z_l,p_h,p_l; | |
509 double y1,t1,t2,r,s,t,u,v,w; | |
510 int i0,i1,i,j,k,yisint,n; | |
511 int hx,hy,ix,iy; | |
512 unsigned lx,ly; | |
513 | |
514 i0 = ((*(int*)&one)>>29)^1; i1=1-i0; | |
515 hx = __HI(x); lx = __LO(x); | |
516 hy = __HI(y); ly = __LO(y); | |
517 ix = hx&0x7fffffff; iy = hy&0x7fffffff; | |
518 | |
519 /* y==zero: x**0 = 1 */ | |
520 if((iy|ly)==0) return one; | |
521 | |
522 /* +-NaN return x+y */ | |
523 if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) || | |
524 iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) | |
525 return x+y; | |
526 | |
527 /* determine if y is an odd int when x < 0 | |
528 * yisint = 0 ... y is not an integer | |
529 * yisint = 1 ... y is an odd int | |
530 * yisint = 2 ... y is an even int | |
531 */ | |
532 yisint = 0; | |
533 if(hx<0) { | |
534 if(iy>=0x43400000) yisint = 2; /* even integer y */ | |
535 else if(iy>=0x3ff00000) { | |
536 k = (iy>>20)-0x3ff; /* exponent */ | |
537 if(k>20) { | |
538 j = ly>>(52-k); | |
539 if((unsigned)(j<<(52-k))==ly) yisint = 2-(j&1); | |
540 } else if(ly==0) { | |
541 j = iy>>(20-k); | |
542 if((j<<(20-k))==iy) yisint = 2-(j&1); | |
543 } | |
544 } | |
545 } | |
546 | |
547 /* special value of y */ | |
548 if(ly==0) { | |
549 if (iy==0x7ff00000) { /* y is +-inf */ | |
550 if(((ix-0x3ff00000)|lx)==0) | |
551 return y - y; /* inf**+-1 is NaN */ | |
552 else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */ | |
553 return (hy>=0)? y: zeroX; | |
554 else /* (|x|<1)**-,+inf = inf,0 */ | |
555 return (hy<0)?-y: zeroX; | |
556 } | |
557 if(iy==0x3ff00000) { /* y is +-1 */ | |
558 if(hy<0) return one/x; else return x; | |
559 } | |
560 if(hy==0x40000000) return x*x; /* y is 2 */ | |
561 if(hy==0x3fe00000) { /* y is 0.5 */ | |
562 if(hx>=0) /* x >= +0 */ | |
563 return sqrt(x); | |
564 } | |
565 } | |
566 | |
567 ax = fabsd(x); | |
568 /* special value of x */ | |
569 if(lx==0) { | |
570 if(ix==0x7ff00000||ix==0||ix==0x3ff00000){ | |
571 z = ax; /*x is +-0,+-inf,+-1*/ | |
572 if(hy<0) z = one/z; /* z = (1/|x|) */ | |
573 if(hx<0) { | |
574 if(((ix-0x3ff00000)|yisint)==0) { | |
1681
126ea7725993
6953477: Increase portability and flexibility of building Hotspot
bobv
parents:
1552
diff
changeset
|
575 #ifdef CAN_USE_NAN_DEFINE |
126ea7725993
6953477: Increase portability and flexibility of building Hotspot
bobv
parents:
1552
diff
changeset
|
576 z = NAN; |
126ea7725993
6953477: Increase portability and flexibility of building Hotspot
bobv
parents:
1552
diff
changeset
|
577 #else |
0 | 578 z = (z-z)/(z-z); /* (-1)**non-int is NaN */ |
1681
126ea7725993
6953477: Increase portability and flexibility of building Hotspot
bobv
parents:
1552
diff
changeset
|
579 #endif |
0 | 580 } else if(yisint==1) |
581 z = -1.0*z; /* (x<0)**odd = -(|x|**odd) */ | |
582 } | |
583 return z; | |
584 } | |
585 } | |
586 | |
587 n = (hx>>31)+1; | |
588 | |
589 /* (x<0)**(non-int) is NaN */ | |
1681
126ea7725993
6953477: Increase portability and flexibility of building Hotspot
bobv
parents:
1552
diff
changeset
|
590 if((n|yisint)==0) |
126ea7725993
6953477: Increase portability and flexibility of building Hotspot
bobv
parents:
1552
diff
changeset
|
591 #ifdef CAN_USE_NAN_DEFINE |
126ea7725993
6953477: Increase portability and flexibility of building Hotspot
bobv
parents:
1552
diff
changeset
|
592 return NAN; |
126ea7725993
6953477: Increase portability and flexibility of building Hotspot
bobv
parents:
1552
diff
changeset
|
593 #else |
126ea7725993
6953477: Increase portability and flexibility of building Hotspot
bobv
parents:
1552
diff
changeset
|
594 return (x-x)/(x-x); |
126ea7725993
6953477: Increase portability and flexibility of building Hotspot
bobv
parents:
1552
diff
changeset
|
595 #endif |
0 | 596 |
597 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ | |
598 if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */ | |
599 | |
600 /* |y| is huge */ | |
601 if(iy>0x41e00000) { /* if |y| > 2**31 */ | |
602 if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */ | |
603 if(ix<=0x3fefffff) return (hy<0)? hugeX*hugeX:tiny*tiny; | |
604 if(ix>=0x3ff00000) return (hy>0)? hugeX*hugeX:tiny*tiny; | |
605 } | |
606 /* over/underflow if x is not close to one */ | |
607 if(ix<0x3fefffff) return (hy<0)? s*hugeX*hugeX:s*tiny*tiny; | |
608 if(ix>0x3ff00000) return (hy>0)? s*hugeX*hugeX:s*tiny*tiny; | |
609 /* now |1-x| is tiny <= 2**-20, suffice to compute | |
610 log(x) by x-x^2/2+x^3/3-x^4/4 */ | |
611 t = ax-one; /* t has 20 trailing zeros */ | |
612 w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25)); | |
613 u = ivln2_h*t; /* ivln2_h has 21 sig. bits */ | |
614 v = t*ivln2_l-w*ivln2; | |
615 t1 = u+v; | |
616 __LO(t1) = 0; | |
617 t2 = v-(t1-u); | |
618 } else { | |
619 double ss,s2,s_h,s_l,t_h,t_l; | |
620 n = 0; | |
621 /* take care subnormal number */ | |
622 if(ix<0x00100000) | |
623 {ax *= two53; n -= 53; ix = __HI(ax); } | |
624 n += ((ix)>>20)-0x3ff; | |
625 j = ix&0x000fffff; | |
626 /* determine interval */ | |
627 ix = j|0x3ff00000; /* normalize ix */ | |
628 if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */ | |
629 else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */ | |
630 else {k=0;n+=1;ix -= 0x00100000;} | |
631 __HI(ax) = ix; | |
632 | |
633 /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ | |
634 u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ | |
635 v = one/(ax+bp[k]); | |
636 ss = u*v; | |
637 s_h = ss; | |
638 __LO(s_h) = 0; | |
639 /* t_h=ax+bp[k] High */ | |
640 t_h = zeroX; | |
641 __HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18); | |
642 t_l = ax - (t_h-bp[k]); | |
643 s_l = v*((u-s_h*t_h)-s_h*t_l); | |
644 /* compute log(ax) */ | |
645 s2 = ss*ss; | |
646 r = s2*s2*(L1X+s2*(L2X+s2*(L3X+s2*(L4X+s2*(L5X+s2*L6X))))); | |
647 r += s_l*(s_h+ss); | |
648 s2 = s_h*s_h; | |
649 t_h = 3.0+s2+r; | |
650 __LO(t_h) = 0; | |
651 t_l = r-((t_h-3.0)-s2); | |
652 /* u+v = ss*(1+...) */ | |
653 u = s_h*t_h; | |
654 v = s_l*t_h+t_l*ss; | |
655 /* 2/(3log2)*(ss+...) */ | |
656 p_h = u+v; | |
657 __LO(p_h) = 0; | |
658 p_l = v-(p_h-u); | |
659 z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ | |
660 z_l = cp_l*p_h+p_l*cp+dp_l[k]; | |
661 /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */ | |
662 t = (double)n; | |
663 t1 = (((z_h+z_l)+dp_h[k])+t); | |
664 __LO(t1) = 0; | |
665 t2 = z_l-(((t1-t)-dp_h[k])-z_h); | |
666 } | |
667 | |
668 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ | |
669 y1 = y; | |
670 __LO(y1) = 0; | |
671 p_l = (y-y1)*t1+y*t2; | |
672 p_h = y1*t1; | |
673 z = p_l+p_h; | |
674 j = __HI(z); | |
675 i = __LO(z); | |
676 if (j>=0x40900000) { /* z >= 1024 */ | |
677 if(((j-0x40900000)|i)!=0) /* if z > 1024 */ | |
678 return s*hugeX*hugeX; /* overflow */ | |
679 else { | |
680 if(p_l+ovt>z-p_h) return s*hugeX*hugeX; /* overflow */ | |
681 } | |
682 } else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */ | |
683 if(((j-0xc090cc00)|i)!=0) /* z < -1075 */ | |
684 return s*tiny*tiny; /* underflow */ | |
685 else { | |
686 if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ | |
687 } | |
688 } | |
689 /* | |
690 * compute 2**(p_h+p_l) | |
691 */ | |
692 i = j&0x7fffffff; | |
693 k = (i>>20)-0x3ff; | |
694 n = 0; | |
695 if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */ | |
696 n = j+(0x00100000>>(k+1)); | |
697 k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */ | |
698 t = zeroX; | |
699 __HI(t) = (n&~(0x000fffff>>k)); | |
700 n = ((n&0x000fffff)|0x00100000)>>(20-k); | |
701 if(j<0) n = -n; | |
702 p_h -= t; | |
703 } | |
704 t = p_l+p_h; | |
705 __LO(t) = 0; | |
706 u = t*lg2_h; | |
707 v = (p_l-(t-p_h))*lg2+t*lg2_l; | |
708 z = u+v; | |
709 w = v-(z-u); | |
710 t = z*z; | |
711 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); | |
712 r = (z*t1)/(t1-two)-(w+z*w); | |
713 z = one-(r-z); | |
714 j = __HI(z); | |
715 j += (n<<20); | |
716 if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */ | |
717 else __HI(z) += (n<<20); | |
718 return s*z; | |
719 } | |
720 | |
721 | |
722 JRT_LEAF(jdouble, SharedRuntime::dpow(jdouble x, jdouble y)) | |
723 return __ieee754_pow(x, y); | |
724 JRT_END | |
725 | |
726 #ifdef WIN32 | |
727 # pragma optimize ( "", on ) | |
728 #endif |