annotate src/share/vm/runtime/sharedRuntimeTrans.cpp @ 2607:008adfd6d850

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