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 }