00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
#if defined(LIBM_SCCS) && !defined(lint)
00014
static char rcsid[] =
"$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $";
00015
#endif
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
#include "math.h"
00063
#include "mathP.h"
00064
00065
#ifdef __STDC__
00066
static const double
00067
#else
00068
static double
00069
#endif
00070 bp[] = {1.0, 1.5,},
00071 dp_h[] = { 0.0, 5.84962487220764160156e-01,},
00072 dp_l[] = { 0.0, 1.35003920212974897128e-08,},
00073 zero = 0.0,
00074 one = 1.0,
00075 two = 2.0,
00076 two53 = 9007199254740992.0,
00077 huge = 1.0e300,
00078 tiny = 1.0e-300,
00079
00080 L1 = 5.99999999999994648725e-01,
00081 L2 = 4.28571428578550184252e-01,
00082 L3 = 3.33333329818377432918e-01,
00083 L4 = 2.72728123808534006489e-01,
00084 L5 = 2.30660745775561754067e-01,
00085 L6 = 2.06975017800338417784e-01,
00086 P1 = 1.66666666666666019037e-01,
00087 P2 = -2.77777777770155933842e-03,
00088 P3 = 6.61375632143793436117e-05,
00089 P4 = -1.65339022054652515390e-06,
00090 P5 = 4.13813679705723846039e-08,
00091 lg2 = 6.93147180559945286227e-01,
00092 lg2_h = 6.93147182464599609375e-01,
00093 lg2_l = -1.90465429995776804525e-09,
00094 ovt = 8.0085662595372944372e-0017,
00095 cp = 9.61796693925975554329e-01,
00096 cp_h = 9.61796700954437255859e-01,
00097 cp_l = -7.02846165095275826516e-09,
00098 ivln2 = 1.44269504088896338700e+00,
00099 ivln2_h = 1.44269502162933349609e+00,
00100 ivln2_l = 1.92596299112661746887e-08;
00101
00102
#ifdef __STDC__
00103
double __ieee754_pow(
double x,
double y)
00104 #
else
00105 double __ieee754_pow(x,y)
00106 double x, y;
00107 #endif
00108 {
00109
double z,ax,z_h,z_l,p_h,p_l;
00110
double y1,t1,t2,r,s,t,u,v,w;
00111 int32_t i,j,k,yisint,n;
00112 int32_t hx,hy,ix,iy;
00113 u_int32_t lx,ly;
00114
00115
EXTRACT_WORDS(hx,lx,x);
00116
EXTRACT_WORDS(hy,ly,y);
00117 ix = hx&0x7fffffff; iy = hy&0x7fffffff;
00118
00119
00120
if((iy|ly)==0)
return one;
00121
00122
00123
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
00124 iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
00125
return x+y;
00126
00127
00128
00129
00130
00131
00132 yisint = 0;
00133
if(hx<0) {
00134
if(iy>=0x43400000) yisint = 2;
00135
else if(iy>=0x3ff00000) {
00136 k = (iy>>20)-0x3ff;
00137
if(k>20) {
00138 j = ly>>(52-k);
00139
if((j<<(52-k))==ly) yisint = 2-(j&1);
00140 }
else if(ly==0) {
00141 j = iy>>(20-k);
00142
if((j<<(20-k))==iy) yisint = 2-(j&1);
00143 }
00144 }
00145 }
00146
00147
00148
if(ly==0) {
00149
if (iy==0x7ff00000) {
00150
if(((ix-0x3ff00000)|lx)==0)
00151
return y - y;
00152
else if (ix >= 0x3ff00000)
00153
return (hy>=0)? y: zero;
00154
else
00155
return (hy<0)?-y: zero;
00156 }
00157
if(iy==0x3ff00000) {
00158
if(hy<0)
return one/x;
else return x;
00159 }
00160
if(hy==0x40000000)
return x*x;
00161
if(hy==0x3fe00000) {
00162
if(hx>=0)
00163
return __ieee754_sqrt(x);
00164 }
00165 }
00166
00167 ax =
fabs(x);
00168
00169
if(lx==0) {
00170
if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
00171 z = ax;
00172
if(hy<0) z = one/z;
00173
if(hx<0) {
00174
if(((ix-0x3ff00000)|yisint)==0) {
00175 z = (z-z)/(z-z);
00176 }
else if(yisint==1)
00177 z = -z;
00178 }
00179
return z;
00180 }
00181 }
00182
00183
00184
if(((((u_int32_t)hx>>31)-1)|yisint)==0)
return (x-x)/(x-x);
00185
00186
00187
if(iy>0x41e00000) {
00188
if(iy>0x43f00000){
00189
if(ix<=0x3fefffff)
return (hy<0)? huge*huge:tiny*tiny;
00190
if(ix>=0x3ff00000)
return (hy>0)? huge*huge:tiny*tiny;
00191 }
00192
00193
if(ix<0x3fefffff)
return (hy<0)? huge*huge:tiny*tiny;
00194
if(ix>0x3ff00000)
return (hy>0)? huge*huge:tiny*tiny;
00195
00196
00197 t = x-1;
00198 w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
00199 u = ivln2_h*t;
00200 v = t*ivln2_l-w*ivln2;
00201 t1 = u+v;
00202
SET_LOW_WORD(t1,0);
00203 t2 = v-(t1-u);
00204 }
else {
00205
double s2,s_h,s_l,t_h,t_l;
00206 n = 0;
00207
00208
if(ix<0x00100000)
00209 {ax *= two53; n -= 53;
GET_HIGH_WORD(ix,ax); }
00210 n += ((ix)>>20)-0x3ff;
00211 j = ix&0x000fffff;
00212
00213 ix = j|0x3ff00000;
00214
if(j<=0x3988E) k=0;
00215
else if(j<0xBB67A) k=1;
00216
else {k=0;n+=1;ix -= 0x00100000;}
00217
SET_HIGH_WORD(ax,ix);
00218
00219
00220 u = ax-
bp[k];
00221 v = one/(ax+
bp[k]);
00222 s = u*v;
00223 s_h = s;
00224
SET_LOW_WORD(s_h,0);
00225
00226 t_h = zero;
00227
SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18));
00228 t_l = ax - (t_h-
bp[k]);
00229 s_l = v*((u-s_h*t_h)-s_h*t_l);
00230
00231 s2 = s*s;
00232 r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
00233 r += s_l*(s_h+s);
00234 s2 = s_h*s_h;
00235 t_h = 3.0+s2+r;
00236
SET_LOW_WORD(t_h,0);
00237 t_l = r-((t_h-3.0)-s2);
00238
00239 u = s_h*t_h;
00240 v = s_l*t_h+t_l*s;
00241
00242 p_h = u+v;
00243
SET_LOW_WORD(p_h,0);
00244 p_l = v-(p_h-u);
00245 z_h = cp_h*p_h;
00246 z_l = cp_l*p_h+p_l*cp+dp_l[k];
00247
00248 t = (
double)n;
00249 t1 = (((z_h+z_l)+dp_h[k])+t);
00250
SET_LOW_WORD(t1,0);
00251 t2 = z_l-(((t1-t)-dp_h[k])-z_h);
00252 }
00253
00254 s = one;
00255
if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0)
00256 s = -one;
00257
00258
00259
y1 = y;
00260
SET_LOW_WORD(
y1,0);
00261 p_l = (y-
y1)*t1+y*t2;
00262 p_h =
y1*t1;
00263 z = p_l+p_h;
00264
EXTRACT_WORDS(j,i,z);
00265
if (j>=0x40900000) {
00266
if(((j-0x40900000)|i)!=0)
00267
return s*huge*huge;
00268
else {
00269
if(p_l+ovt>z-p_h)
return s*huge*huge;
00270 }
00271 }
else if((j&0x7fffffff)>=0x4090cc00 ) {
00272
if(((j-0xc090cc00)|i)!=0)
00273
return s*tiny*tiny;
00274
else {
00275
if(p_l<=z-p_h)
return s*tiny*tiny;
00276 }
00277 }
00278
00279
00280
00281 i = j&0x7fffffff;
00282 k = (i>>20)-0x3ff;
00283 n = 0;
00284
if(i>0x3fe00000) {
00285 n = j+(0x00100000>>(k+1));
00286 k = ((n&0x7fffffff)>>20)-0x3ff;
00287 t = zero;
00288
SET_HIGH_WORD(t,n&~(0x000fffff>>k));
00289 n = ((n&0x000fffff)|0x00100000)>>(20-k);
00290
if(j<0) n = -n;
00291 p_h -= t;
00292 }
00293 t = p_l+p_h;
00294
SET_LOW_WORD(t,0);
00295 u = t*lg2_h;
00296 v = (p_l-(t-p_h))*lg2+t*lg2_l;
00297 z = u+v;
00298 w = v-(z-u);
00299 t = z*z;
00300 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
00301 r = (z*t1)/(t1-two)-(w+z*w);
00302 z = one-(r-z);
00303
GET_HIGH_WORD(j,z);
00304 j += (n<<20);
00305
if((j>>20)<=0) z =
scalbn(z,n);
00306
else SET_HIGH_WORD(z,j);
00307
return s*z;
00308 }