annotate src/share/vm/runtime/sharedRuntimeTrans.cpp @ 1145:e018e6884bd8

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