1/* 2 * IBM Accurate Mathematical Library 3 * written by International Business Machines Corp. 4 * Copyright (C) 2001, 2002 Free Software Foundation 5 * 6 * This program is free software; you can redistribute it and/or modify 7 * it under the terms of the GNU Lesser General Public License as published by 8 * the Free Software Foundation; either version 2.1 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU Lesser General Public License for more details. 15 * 16 * You should have received a copy of the GNU Lesser General Public License 17 * along with this program; if not, write to the Free Software 18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 19 */ 20/***************************************************************************/ 21/* MODULE_NAME: upow.c */ 22/* */ 23/* FUNCTIONS: upow */ 24/* power1 */ 25/* log2 */ 26/* log1 */ 27/* checkint */ 28/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h */ 29/* halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c */ 30/* uexp.c upow.c */ 31/* root.tbl uexp.tbl upow.tbl */ 32/* An ultimate power routine. Given two IEEE double machine numbers y,x */ 33/* it computes the correctly rounded (to nearest) value of x^y. */ 34/* Assumption: Machine arithmetic operations are performed in */ 35/* round to nearest mode of IEEE 754 standard. */ 36/* */ 37/***************************************************************************/ 38#include "endian.h" 39#include "upow.h" 40#include "dla.h" 41#include "mydefs.h" 42#include "MathLib.h" 43#include "upow.tbl" 44#include "math_private.h" 45 46 47double __exp1(double x, double xx, double error); 48static double log1(double x, double *delta, double *error); 49static double log2(double x, double *delta, double *error); 50double __slowpow(double x, double y,double z); 51static double power1(double x, double y); 52static int checkint(double x); 53 54/***************************************************************************/ 55/* An ultimate power routine. Given two IEEE double machine numbers y,x */ 56/* it computes the correctly rounded (to nearest) value of X^y. */ 57/***************************************************************************/ 58double __ieee754_pow(double x, double y) { 59 double z,a,aa,error, t,a1,a2,y1,y2; 60#if 0 61 double gor=1.0; 62#endif 63 mynumber u,v; 64 int k; 65 int4 qx,qy; 66 v.x=y; 67 u.x=x; 68 if (v.i[LOW_HALF] == 0) { /* of y */ 69 qx = u.i[HIGH_HALF]&0x7fffffff; 70 /* Checking if x is not too small to compute */ 71 if (((qx==0x7ff00000)&&(u.i[LOW_HALF]!=0))||(qx>0x7ff00000)) return NaNQ.x; 72 if (y == 1.0) return x; 73 if (y == 2.0) return x*x; 74 if (y == -1.0) return 1.0/x; 75 if (y == 0) return 1.0; 76 } 77 /* else */ 78 if(((u.i[HIGH_HALF]>0 && u.i[HIGH_HALF]<0x7ff00000)|| /* x>0 and not x->0 */ 79 (u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0)) && 80 /* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */ 81 (v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) { /* if y<-1 or y>1 */ 82 z = log1(x,&aa,&error); /* x^y =e^(y log (X)) */ 83 t = y*134217729.0; 84 y1 = t - (t-y); 85 y2 = y - y1; 86 t = z*134217729.0; 87 a1 = t - (t-z); 88 a2 = (z - a1)+aa; 89 a = y1*a1; 90 aa = y2*a1 + y*a2; 91 a1 = a+aa; 92 a2 = (a-a1)+aa; 93 error = error*ABS(y); 94 t = __exp1(a1,a2,1.9e16*error); /* return -10 or 0 if wasn't computed exactly */ 95 return (t>0)?t:power1(x,y); 96 } 97 98 if (x == 0) { 99 if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0) 100 || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000) 101 return y; 102 if (ABS(y) > 1.0e20) return (y>0)?0:INF.x; 103 k = checkint(y); 104 if (k == -1) 105 return y < 0 ? 1.0/x : x; 106 else 107 return y < 0 ? 1.0/ABS(x) : 0.0; /* return 0 */ 108 } 109 /* if x<0 */ 110 if (u.i[HIGH_HALF] < 0) { 111 k = checkint(y); 112 if (k==0) { 113 if ((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] == 0) { 114 if (x == -1.0) return 1.0; 115 else if (x > -1.0) return v.i[HIGH_HALF] < 0 ? INF.x : 0.0; 116 else return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x; 117 } 118 else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0) 119 return y < 0 ? 0.0 : INF.x; 120 return NaNQ.x; /* y not integer and x<0 */ 121 } 122 else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0) 123 { 124 if (k < 0) 125 return y < 0 ? nZERO.x : nINF.x; 126 else 127 return y < 0 ? 0.0 : INF.x; 128 } 129 return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */ 130 } 131 /* x>0 */ 132 qx = u.i[HIGH_HALF]&0x7fffffff; /* no sign */ 133 qy = v.i[HIGH_HALF]&0x7fffffff; /* no sign */ 134 135 if (qx > 0x7ff00000 || (qx == 0x7ff00000 && u.i[LOW_HALF] != 0)) return NaNQ.x; 136 /* if 0<x<2^-0x7fe */ 137 if (qy > 0x7ff00000 || (qy == 0x7ff00000 && v.i[LOW_HALF] != 0)) 138 return x == 1.0 ? 1.0 : NaNQ.x; 139 /* if y<2^-0x7fe */ 140 141 if (qx == 0x7ff00000) /* x= 2^-0x3ff */ 142 {if (y == 0) return NaNQ.x; 143 return (y>0)?x:0; } 144 145 if (qy > 0x45f00000 && qy < 0x7ff00000) { 146 if (x == 1.0) return 1.0; 147 if (y>0) return (x>1.0)?INF.x:0; 148 if (y<0) return (x<1.0)?INF.x:0; 149 } 150 151 if (x == 1.0) return 1.0; 152 if (y>0) return (x>1.0)?INF.x:0; 153 if (y<0) return (x<1.0)?INF.x:0; 154 return 0; /* unreachable, to make the compiler happy */ 155} 156 157/**************************************************************************/ 158/* Computing x^y using more accurate but more slow log routine */ 159/**************************************************************************/ 160static double power1(double x, double y) { 161 double z,a,aa,error, t,a1,a2,y1,y2; 162 z = log2(x,&aa,&error); 163 t = y*134217729.0; 164 y1 = t - (t-y); 165 y2 = y - y1; 166 t = z*134217729.0; 167 a1 = t - (t-z); 168 a2 = z - a1; 169 a = y*z; 170 aa = ((y1*a1-a)+y1*a2+y2*a1)+y2*a2+aa*y; 171 a1 = a+aa; 172 a2 = (a-a1)+aa; 173 error = error*ABS(y); 174 t = __exp1(a1,a2,1.9e16*error); 175 return (t >= 0)?t:__slowpow(x,y,z); 176} 177 178/****************************************************************************/ 179/* Computing log(x) (x is left argument). The result is the returned double */ 180/* + the parameter delta. */ 181/* The result is bounded by error (rightmost argument) */ 182/****************************************************************************/ 183static double log1(double x, double *delta, double *error) { 184 int i,j,m; 185#if 0 186 int n; 187#endif 188 double uu,vv,eps,nx,e,e1,e2,t,t1,t2,res,add=0; 189#if 0 190 double cor; 191#endif 192 mynumber u,v; 193#ifdef BIG_ENDI 194 mynumber 195/**/ two52 = {{0x43300000, 0x00000000}}; /* 2**52 */ 196#else 197#ifdef LITTLE_ENDI 198 mynumber 199/**/ two52 = {{0x00000000, 0x43300000}}; /* 2**52 */ 200#endif 201#endif 202 203 u.x = x; 204 m = u.i[HIGH_HALF]; 205 *error = 0; 206 *delta = 0; 207 if (m < 0x00100000) /* 1<x<2^-1007 */ 208 { x = x*t52.x; add = -52.0; u.x = x; m = u.i[HIGH_HALF];} 209 210 if ((m&0x000fffff) < 0x0006a09e) 211 {u.i[HIGH_HALF] = (m&0x000fffff)|0x3ff00000; two52.i[LOW_HALF]=(m>>20); } 212 else 213 {u.i[HIGH_HALF] = (m&0x000fffff)|0x3fe00000; two52.i[LOW_HALF]=(m>>20)+1; } 214 215 v.x = u.x + bigu.x; 216 uu = v.x - bigu.x; 217 i = (v.i[LOW_HALF]&0x000003ff)<<2; 218 if (two52.i[LOW_HALF] == 1023) /* nx = 0 */ 219 { 220 if (i > 1192 && i < 1208) /* |x-1| < 1.5*2**-10 */ 221 { 222 t = x - 1.0; 223 t1 = (t+5.0e6)-5.0e6; 224 t2 = t-t1; 225 e1 = t - 0.5*t1*t1; 226 e2 = t*t*t*(r3+t*(r4+t*(r5+t*(r6+t*(r7+t*r8)))))-0.5*t2*(t+t1); 227 res = e1+e2; 228 *error = 1.0e-21*ABS(t); 229 *delta = (e1-res)+e2; 230 return res; 231 } /* |x-1| < 1.5*2**-10 */ 232 else 233 { 234 v.x = u.x*(ui.x[i]+ui.x[i+1])+bigv.x; 235 vv = v.x-bigv.x; 236 j = v.i[LOW_HALF]&0x0007ffff; 237 j = j+j+j; 238 eps = u.x - uu*vv; 239 e1 = eps*ui.x[i]; 240 e2 = eps*(ui.x[i+1]+vj.x[j]*(ui.x[i]+ui.x[i+1])); 241 e = e1+e2; 242 e2 = ((e1-e)+e2); 243 t=ui.x[i+2]+vj.x[j+1]; 244 t1 = t+e; 245 t2 = (((t-t1)+e)+(ui.x[i+3]+vj.x[j+2]))+e2+e*e*(p2+e*(p3+e*p4)); 246 res=t1+t2; 247 *error = 1.0e-24; 248 *delta = (t1-res)+t2; 249 return res; 250 } 251 } /* nx = 0 */ 252 else /* nx != 0 */ 253 { 254 eps = u.x - uu; 255 nx = (two52.x - two52e.x)+add; 256 e1 = eps*ui.x[i]; 257 e2 = eps*ui.x[i+1]; 258 e=e1+e2; 259 e2 = (e1-e)+e2; 260 t=nx*ln2a.x+ui.x[i+2]; 261 t1=t+e; 262 t2=(((t-t1)+e)+nx*ln2b.x+ui.x[i+3]+e2)+e*e*(q2+e*(q3+e*(q4+e*(q5+e*q6)))); 263 res = t1+t2; 264 *error = 1.0e-21; 265 *delta = (t1-res)+t2; 266 return res; 267 } /* nx != 0 */ 268} 269 270/****************************************************************************/ 271/* More slow but more accurate routine of log */ 272/* Computing log(x)(x is left argument).The result is return double + delta.*/ 273/* The result is bounded by error (right argument) */ 274/****************************************************************************/ 275static double log2(double x, double *delta, double *error) { 276 int i,j,m; 277#if 0 278 int n; 279#endif 280 double uu,vv,eps,nx,e,e1,e2,t,t1,t2,res,add=0; 281#if 0 282 double cor; 283#endif 284 double ou1,ou2,lu1,lu2,ov,lv1,lv2,a,a1,a2; 285 double y,yy,z,zz,j1,j2,j3,j4,j5,j6,j7,j8; 286 mynumber u,v; 287#ifdef BIG_ENDI 288 mynumber 289/**/ two52 = {{0x43300000, 0x00000000}}; /* 2**52 */ 290#else 291#ifdef LITTLE_ENDI 292 mynumber 293/**/ two52 = {{0x00000000, 0x43300000}}; /* 2**52 */ 294#endif 295#endif 296 297 u.x = x; 298 m = u.i[HIGH_HALF]; 299 *error = 0; 300 *delta = 0; 301 add=0; 302 if (m<0x00100000) { /* x < 2^-1022 */ 303 x = x*t52.x; add = -52.0; u.x = x; m = u.i[HIGH_HALF]; } 304 305 if ((m&0x000fffff) < 0x0006a09e) 306 {u.i[HIGH_HALF] = (m&0x000fffff)|0x3ff00000; two52.i[LOW_HALF]=(m>>20); } 307 else 308 {u.i[HIGH_HALF] = (m&0x000fffff)|0x3fe00000; two52.i[LOW_HALF]=(m>>20)+1; } 309 310 v.x = u.x + bigu.x; 311 uu = v.x - bigu.x; 312 i = (v.i[LOW_HALF]&0x000003ff)<<2; 313 /*------------------------------------- |x-1| < 2**-11------------------------------- */ 314 if ((two52.i[LOW_HALF] == 1023) && (i == 1200)) 315 { 316 t = x - 1.0; 317 EMULV(t,s3,y,yy,j1,j2,j3,j4,j5); 318 ADD2(-0.5,0,y,yy,z,zz,j1,j2); 319 MUL2(t,0,z,zz,y,yy,j1,j2,j3,j4,j5,j6,j7,j8); 320 MUL2(t,0,y,yy,z,zz,j1,j2,j3,j4,j5,j6,j7,j8); 321 322 e1 = t+z; 323 e2 = (((t-e1)+z)+zz)+t*t*t*(ss3+t*(s4+t*(s5+t*(s6+t*(s7+t*s8))))); 324 res = e1+e2; 325 *error = 1.0e-25*ABS(t); 326 *delta = (e1-res)+e2; 327 return res; 328 } 329 /*----------------------------- |x-1| > 2**-11 -------------------------- */ 330 else 331 { /*Computing log(x) according to log table */ 332 nx = (two52.x - two52e.x)+add; 333 ou1 = ui.x[i]; 334 ou2 = ui.x[i+1]; 335 lu1 = ui.x[i+2]; 336 lu2 = ui.x[i+3]; 337 v.x = u.x*(ou1+ou2)+bigv.x; 338 vv = v.x-bigv.x; 339 j = v.i[LOW_HALF]&0x0007ffff; 340 j = j+j+j; 341 eps = u.x - uu*vv; 342 ov = vj.x[j]; 343 lv1 = vj.x[j+1]; 344 lv2 = vj.x[j+2]; 345 a = (ou1+ou2)*(1.0+ov); 346 a1 = (a+1.0e10)-1.0e10; 347 a2 = a*(1.0-a1*uu*vv); 348 e1 = eps*a1; 349 e2 = eps*a2; 350 e = e1+e2; 351 e2 = (e1-e)+e2; 352 t=nx*ln2a.x+lu1+lv1; 353 t1 = t+e; 354 t2 = (((t-t1)+e)+(lu2+lv2+nx*ln2b.x+e2))+e*e*(p2+e*(p3+e*p4)); 355 res=t1+t2; 356 *error = 1.0e-27; 357 *delta = (t1-res)+t2; 358 return res; 359 } 360} 361 362/**********************************************************************/ 363/* Routine receives a double x and checks if it is an integer. If not */ 364/* it returns 0, else it returns 1 if even or -1 if odd. */ 365/**********************************************************************/ 366static int checkint(double x) { 367 union {int4 i[2]; double x;} u; 368 int k,m,n; 369#if 0 370 int l; 371#endif 372 u.x = x; 373 m = u.i[HIGH_HALF]&0x7fffffff; /* no sign */ 374 if (m >= 0x7ff00000) return 0; /* x is +/-inf or NaN */ 375 if (m >= 0x43400000) return 1; /* |x| >= 2**53 */ 376 if (m < 0x40000000) return 0; /* |x| < 2, can not be 0 or 1 */ 377 n = u.i[LOW_HALF]; 378 k = (m>>20)-1023; /* 1 <= k <= 52 */ 379 if (k == 52) return (n&1)? -1:1; /* odd or even*/ 380 if (k>20) { 381 if (n<<(k-20)) return 0; /* if not integer */ 382 return (n<<(k-21))?-1:1; 383 } 384 if (n) return 0; /*if not integer*/ 385 if (k == 20) return (m&1)? -1:1; 386 if (m<<(k+12)) return 0; 387 return (m<<(k+11))?-1:1; 388} 389