sysdeps/ieee754/dbl-64/e_pow.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2014 Free Software Foundation, Inc.
   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, see <http://www.gnu.org/licenses/>.
  18  */
  19 /***************************************************************************/
  20 /*  MODULE_NAME: upow.c                                                    */
  21 /*                                                                         */
  22 /*  FUNCTIONS: upow                                                        */
  23 /*             power1                                                      */
  24 /*             my_log2                                                     */
  25 /*             log1                                                        */
  26 /*             checkint                                                    */
  27 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h                             */
  28 /*               halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c       */
  29 /*                          uexp.c  upow.c                                 */
  30 /*               root.tbl uexp.tbl upow.tbl                                */
  31 /* An ultimate power routine. Given two IEEE double machine numbers y,x    */
  32 /* it computes the correctly rounded (to nearest) value of x^y.            */
  33 /* Assumption: Machine arithmetic operations are performed in              */
  34 /* round to nearest mode of IEEE 754 standard.                             */
  35 /*                                                                         */
  36 /***************************************************************************/
  37 #include <math.h>
  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 #include <fenv.h>
  46
  47 #ifndef SECTION
  48 # define SECTION
  49 #endif
  50
  51 static const double huge = 1.0e300, tiny = 1.0e-300;
  52
  53 double __exp1 (double x, double xx, double error);
  54 static double log1 (double x, double *delta, double *error);
  55 static double my_log2 (double x, double *delta, double *error);
  56 double __slowpow (double x, double y, double z);
  57 static double power1 (double x, double y);
  58 static int checkint (double x);
  59
  60 /* An ultimate power routine. Given two IEEE double machine numbers y, x it
  61    computes the correctly rounded (to nearest) value of X^y.  */
  62 double
  63 SECTION
  64 __ieee754_pow (double x, double y)
  65 {
  66   double z, a, aa, error, t, a1, a2, y1, y2;
  67   mynumber u, v;
  68   int k;
  69   int4 qx, qy;
  70   v.x = y;
  71   u.x = x;
  72   if (v.i[LOW_HALF] == 0)
  73     {                           /* of y */
  74       qx = u.i[HIGH_HALF] & 0x7fffffff;
  75       /* Is x a NaN?  */
  76       if (((qx == 0x7ff00000) && (u.i[LOW_HALF] != 0)) || (qx > 0x7ff00000))
  77         return x;
  78       if (y == 1.0)
  79         return x;
  80       if (y == 2.0)
  81         return x * x;
  82       if (y == -1.0)
  83         return 1.0 / x;
  84       if (y == 0)
  85         return 1.0;
  86     }
  87   /* else */
  88   if (((u.i[HIGH_HALF] > 0 && u.i[HIGH_HALF] < 0x7ff00000) ||   /* x>0 and not x->0 */
  89        (u.i[HIGH_HALF] == 0 && u.i[LOW_HALF] != 0)) &&
  90       /*   2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
  91       (v.i[HIGH_HALF] & 0x7fffffff) < 0x4ff00000)
  92     {                           /* if y<-1 or y>1   */
  93       double retval;
  94
  95       {
  96         SET_RESTORE_ROUND (FE_TONEAREST);
  97
  98         /* Avoid internal underflow for tiny y.  The exact value of y does
  99            not matter if |y| <= 2**-64.  */
 100         if (ABS (y) < 0x1p-64)
 101           y = y < 0 ? -0x1p-64 : 0x1p-64;
 102         z = log1 (x, &aa, &error);      /* x^y  =e^(y log (X)) */
 103         t = y * CN;
 104         y1 = t - (t - y);
 105         y2 = y - y1;
 106         t = z * CN;
 107         a1 = t - (t - z);
 108         a2 = (z - a1) + aa;
 109         a = y1 * a1;
 110         aa = y2 * a1 + y * a2;
 111         a1 = a + aa;
 112         a2 = (a - a1) + aa;
 113         error = error * ABS (y);
 114         t = __exp1 (a1, a2, 1.9e16 * error);    /* return -10 or 0 if wasn't computed exactly */
 115         retval = (t > 0) ? t : power1 (x, y);
 116       }
 117
 118       if (__isinf (retval))
 119         retval = huge * huge;
 120       else if (retval == 0)
 121         retval = tiny * tiny;
 122       return retval;
 123     }
 124
 125   if (x == 0)
 126     {
 127       if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
 128           || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)        /* NaN */
 129         return y;
 130       if (ABS (y) > 1.0e20)
 131         return (y > 0) ? 0 : 1.0 / 0.0;
 132       k = checkint (y);
 133       if (k == -1)
 134         return y < 0 ? 1.0 / x : x;
 135       else
 136         return y < 0 ? 1.0 / 0.0 : 0.0; /* return 0 */
 137     }
 138
 139   qx = u.i[HIGH_HALF] & 0x7fffffff;     /*   no sign   */
 140   qy = v.i[HIGH_HALF] & 0x7fffffff;     /*   no sign   */
 141
 142   if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0))      /* NaN */
 143     return x;
 144   if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0))      /* NaN */
 145     return x == 1.0 ? 1.0 : y;
 146
 147   /* if x<0 */
 148   if (u.i[HIGH_HALF] < 0)
 149     {
 150       k = checkint (y);
 151       if (k == 0)
 152         {
 153           if (qy == 0x7ff00000)
 154             {
 155               if (x == -1.0)
 156                 return 1.0;
 157               else if (x > -1.0)
 158                 return v.i[HIGH_HALF] < 0 ? INF.x : 0.0;
 159               else
 160                 return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x;
 161             }
 162           else if (qx == 0x7ff00000)
 163             return y < 0 ? 0.0 : INF.x;
 164           return (x - x) / (x - x);     /* y not integer and x<0 */
 165         }
 166       else if (qx == 0x7ff00000)
 167         {
 168           if (k < 0)
 169             return y < 0 ? nZERO.x : nINF.x;
 170           else
 171             return y < 0 ? 0.0 : INF.x;
 172         }
 173       /* if y even or odd */
 174       if (k == 1)
 175         return __ieee754_pow (-x, y);
 176       else
 177         {
 178           double retval;
 179           {
 180             SET_RESTORE_ROUND (FE_TONEAREST);
 181             retval = -__ieee754_pow (-x, y);
 182           }
 183           if (__isinf (retval))
 184             retval = -huge * huge;
 185           else if (retval == 0)
 186             retval = -tiny * tiny;
 187           return retval;
 188         }
 189     }
 190   /* x>0 */
 191
 192   if (qx == 0x7ff00000)         /* x= 2^-0x3ff */
 193     return y > 0 ? x : 0;
 194
 195   if (qy > 0x45f00000 && qy < 0x7ff00000)
 196     {
 197       if (x == 1.0)
 198         return 1.0;
 199       if (y > 0)
 200         return (x > 1.0) ? huge * huge : tiny * tiny;
 201       if (y < 0)
 202         return (x < 1.0) ? huge * huge : tiny * tiny;
 203     }
 204
 205   if (x == 1.0)
 206     return 1.0;
 207   if (y > 0)
 208     return (x > 1.0) ? INF.x : 0;
 209   if (y < 0)
 210     return (x < 1.0) ? INF.x : 0;
 211   return 0;                     /* unreachable, to make the compiler happy */
 212 }
 213
 214 #ifndef __ieee754_pow
 215 strong_alias (__ieee754_pow, __pow_finite)
 216 #endif
 217
 218 /* Compute x^y using more accurate but more slow log routine.  */
 219 static double
 220 SECTION
 221 power1 (double x, double y)
 222 {
 223   double z, a, aa, error, t, a1, a2, y1, y2;
 224   z = my_log2 (x, &aa, &error);
 225   t = y * CN;
 226   y1 = t - (t - y);
 227   y2 = y - y1;
 228   t = z * CN;
 229   a1 = t - (t - z);
 230   a2 = z - a1;
 231   a = y * z;
 232   aa = ((y1 * a1 - a) + y1 * a2 + y2 * a1) + y2 * a2 + aa * y;
 233   a1 = a + aa;
 234   a2 = (a - a1) + aa;
 235   error = error * ABS (y);
 236   t = __exp1 (a1, a2, 1.9e16 * error);
 237   return (t >= 0) ? t : __slowpow (x, y, z);
 238 }
 239
 240 /* Compute log(x) (x is left argument). The result is the returned double + the
 241    parameter DELTA.  The result is bounded by ERROR.  */
 242 static double
 243 SECTION
 244 log1 (double x, double *delta, double *error)
 245 {
 246   int i, j, m;
 247   double uu, vv, eps, nx, e, e1, e2, t, t1, t2, res, add = 0;
 248   mynumber u, v;
 249 #ifdef BIG_ENDI
 250   mynumber /**/ two52 = {{0x43300000, 0x00000000}};     /* 2**52  */
 251 #else
 252 # ifdef LITTLE_ENDI
 253   mynumber /**/ two52 = {{0x00000000, 0x43300000}};     /* 2**52  */
 254 # endif
 255 #endif
 256
 257   u.x = x;
 258   m = u.i[HIGH_HALF];
 259   *error = 0;
 260   *delta = 0;
 261   if (m < 0x00100000)           /*  1<x<2^-1007 */
 262     {
 263       x = x * t52.x;
 264       add = -52.0;
 265       u.x = x;
 266       m = u.i[HIGH_HALF];
 267     }
 268
 269   if ((m & 0x000fffff) < 0x0006a09e)
 270     {
 271       u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3ff00000;
 272       two52.i[LOW_HALF] = (m >> 20);
 273     }
 274   else
 275     {
 276       u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3fe00000;
 277       two52.i[LOW_HALF] = (m >> 20) + 1;
 278     }
 279
 280   v.x = u.x + bigu.x;
 281   uu = v.x - bigu.x;
 282   i = (v.i[LOW_HALF] & 0x000003ff) << 2;
 283   if (two52.i[LOW_HALF] == 1023)        /* nx = 0              */
 284     {
 285       if (i > 1192 && i < 1208) /* |x-1| < 1.5*2**-10  */
 286         {
 287           t = x - 1.0;
 288           t1 = (t + 5.0e6) - 5.0e6;
 289           t2 = t - t1;
 290           e1 = t - 0.5 * t1 * t1;
 291           e2 = (t * t * t * (r3 + t * (r4 + t * (r5 + t * (r6 + t
 292                                                            * (r7 + t * r8)))))
 293                 - 0.5 * t2 * (t + t1));
 294           res = e1 + e2;
 295           *error = 1.0e-21 * ABS (t);
 296           *delta = (e1 - res) + e2;
 297           return res;
 298         }                       /* |x-1| < 1.5*2**-10  */
 299       else
 300         {
 301           v.x = u.x * (ui.x[i] + ui.x[i + 1]) + bigv.x;
 302           vv = v.x - bigv.x;
 303           j = v.i[LOW_HALF] & 0x0007ffff;
 304           j = j + j + j;
 305           eps = u.x - uu * vv;
 306           e1 = eps * ui.x[i];
 307           e2 = eps * (ui.x[i + 1] + vj.x[j] * (ui.x[i] + ui.x[i + 1]));
 308           e = e1 + e2;
 309           e2 = ((e1 - e) + e2);
 310           t = ui.x[i + 2] + vj.x[j + 1];
 311           t1 = t + e;
 312           t2 = ((((t - t1) + e) + (ui.x[i + 3] + vj.x[j + 2])) + e2 + e * e
 313                 * (p2 + e * (p3 + e * p4)));
 314           res = t1 + t2;
 315           *error = 1.0e-24;
 316           *delta = (t1 - res) + t2;
 317           return res;
 318         }
 319     }                           /* nx = 0 */
 320   else                          /* nx != 0   */
 321     {
 322       eps = u.x - uu;
 323       nx = (two52.x - two52e.x) + add;
 324       e1 = eps * ui.x[i];
 325       e2 = eps * ui.x[i + 1];
 326       e = e1 + e2;
 327       e2 = (e1 - e) + e2;
 328       t = nx * ln2a.x + ui.x[i + 2];
 329       t1 = t + e;
 330       t2 = ((((t - t1) + e) + nx * ln2b.x + ui.x[i + 3] + e2) + e * e
 331             * (q2 + e * (q3 + e * (q4 + e * (q5 + e * q6)))));
 332       res = t1 + t2;
 333       *error = 1.0e-21;
 334       *delta = (t1 - res) + t2;
 335       return res;
 336     }                           /* nx != 0   */
 337 }
 338
 339 /* Slower but more accurate routine of log.  The returned result is double +
 340    DELTA.  The result is bounded by ERROR.  */
 341 static double
 342 SECTION
 343 my_log2 (double x, double *delta, double *error)
 344 {
 345   int i, j, m;
 346   double uu, vv, eps, nx, e, e1, e2, t, t1, t2, res, add = 0;
 347   double ou1, ou2, lu1, lu2, ov, lv1, lv2, a, a1, a2;
 348   double y, yy, z, zz, j1, j2, j7, j8;
 349 #ifndef DLA_FMS
 350   double j3, j4, j5, j6;
 351 #endif
 352   mynumber u, v;
 353 #ifdef BIG_ENDI
 354   mynumber /**/ two52 = {{0x43300000, 0x00000000}};     /* 2**52  */
 355 #else
 356 # ifdef LITTLE_ENDI
 357   mynumber /**/ two52 = {{0x00000000, 0x43300000}};     /* 2**52  */
 358 # endif
 359 #endif
 360
 361   u.x = x;
 362   m = u.i[HIGH_HALF];
 363   *error = 0;
 364   *delta = 0;
 365   add = 0;
 366   if (m < 0x00100000)
 367     {                           /* x < 2^-1022 */
 368       x = x * t52.x;
 369       add = -52.0;
 370       u.x = x;
 371       m = u.i[HIGH_HALF];
 372     }
 373
 374   if ((m & 0x000fffff) < 0x0006a09e)
 375     {
 376       u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3ff00000;
 377       two52.i[LOW_HALF] = (m >> 20);
 378     }
 379   else
 380     {
 381       u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3fe00000;
 382       two52.i[LOW_HALF] = (m >> 20) + 1;
 383     }
 384
 385   v.x = u.x + bigu.x;
 386   uu = v.x - bigu.x;
 387   i = (v.i[LOW_HALF] & 0x000003ff) << 2;
 388   /*------------------------------------- |x-1| < 2**-11-------------------------------  */
 389   if ((two52.i[LOW_HALF] == 1023) && (i == 1200))
 390     {
 391       t = x - 1.0;
 392       EMULV (t, s3, y, yy, j1, j2, j3, j4, j5);
 393       ADD2 (-0.5, 0, y, yy, z, zz, j1, j2);
 394       MUL2 (t, 0, z, zz, y, yy, j1, j2, j3, j4, j5, j6, j7, j8);
 395       MUL2 (t, 0, y, yy, z, zz, j1, j2, j3, j4, j5, j6, j7, j8);
 396
 397       e1 = t + z;
 398       e2 = ((((t - e1) + z) + zz) + t * t * t
 399             * (ss3 + t * (s4 + t * (s5 + t * (s6 + t * (s7 + t * s8))))));
 400       res = e1 + e2;
 401       *error = 1.0e-25 * ABS (t);
 402       *delta = (e1 - res) + e2;
 403       return res;
 404     }
 405   /*----------------------------- |x-1| > 2**-11  --------------------------  */
 406   else
 407     {                           /*Computing log(x) according to log table                        */
 408       nx = (two52.x - two52e.x) + add;
 409       ou1 = ui.x[i];
 410       ou2 = ui.x[i + 1];
 411       lu1 = ui.x[i + 2];
 412       lu2 = ui.x[i + 3];
 413       v.x = u.x * (ou1 + ou2) + bigv.x;
 414       vv = v.x - bigv.x;
 415       j = v.i[LOW_HALF] & 0x0007ffff;
 416       j = j + j + j;
 417       eps = u.x - uu * vv;
 418       ov = vj.x[j];
 419       lv1 = vj.x[j + 1];
 420       lv2 = vj.x[j + 2];
 421       a = (ou1 + ou2) * (1.0 + ov);
 422       a1 = (a + 1.0e10) - 1.0e10;
 423       a2 = a * (1.0 - a1 * uu * vv);
 424       e1 = eps * a1;
 425       e2 = eps * a2;
 426       e = e1 + e2;
 427       e2 = (e1 - e) + e2;
 428       t = nx * ln2a.x + lu1 + lv1;
 429       t1 = t + e;
 430       t2 = ((((t - t1) + e) + (lu2 + lv2 + nx * ln2b.x + e2)) + e * e
 431             * (p2 + e * (p3 + e * p4)));
 432       res = t1 + t2;
 433       *error = 1.0e-27;
 434       *delta = (t1 - res) + t2;
 435       return res;
 436     }
 437 }
 438
 439 /* This function receives a double x and checks if it is an integer.  If not,
 440    it returns 0, else it returns 1 if even or -1 if odd.  */
 441 static int
 442 SECTION
 443 checkint (double x)
 444 {
 445   union
 446   {
 447     int4 i[2];
 448     double x;
 449   } u;
 450   int k, m, n;
 451   u.x = x;
 452   m = u.i[HIGH_HALF] & 0x7fffffff;      /* no sign */
 453   if (m >= 0x7ff00000)
 454     return 0;                   /*  x is +/-inf or NaN  */
 455   if (m >= 0x43400000)
 456     return 1;                   /*  |x| >= 2**53   */
 457   if (m < 0x40000000)
 458     return 0;                   /* |x| < 2,  can not be 0 or 1  */
 459   n = u.i[LOW_HALF];
 460   k = (m >> 20) - 1023;         /*  1 <= k <= 52   */
 461   if (k == 52)
 462     return (n & 1) ? -1 : 1;    /* odd or even */
 463   if (k > 20)
 464     {
 465       if (n << (k - 20))
 466         return 0;               /* if not integer */
 467       return (n << (k - 21)) ? -1 : 1;
 468     }
 469   if (n)
 470     return 0;                   /*if  not integer */
 471   if (k == 20)
 472     return (m & 1) ? -1 : 1;
 473   if (m << (k + 12))
 474     return 0;
 475   return (m << (k + 11)) ? -1 : 1;
 476 }