sysdeps/ieee754/dbl-64/e_exp.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2015 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:uexp.c                                                     */
  21 /*                                                                         */
  22 /*  FUNCTION:uexp                                                          */
  23 /*           exp1                                                          */
  24 /*                                                                         */
  25 /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h                       */
  26 /*              mpa.c mpexp.x slowexp.c                                    */
  27 /*                                                                         */
  28 /* An ultimate exp routine. Given an IEEE double machine number x          */
  29 /* it computes the correctly rounded (to nearest) value of e^x             */
  30 /* Assumption: Machine arithmetic operations are performed in              */
  31 /* round to nearest mode of IEEE 754 standard.                             */
  32 /*                                                                         */
  33 /***************************************************************************/
  34
  35 #include <math.h>
  36 #include "endian.h"
  37 #include "uexp.h"
  38 #include "mydefs.h"
  39 #include "MathLib.h"
  40 #include "uexp.tbl"
  41 #include <math_private.h>
  42 #include <fenv.h>
  43 #include <float.h>
  44
  45 #ifndef SECTION
  46 # define SECTION
  47 #endif
  48
  49 double __slowexp (double);
  50
  51 /* An ultimate exp routine. Given an IEEE double machine number x it computes
  52    the correctly rounded (to nearest) value of e^x.  */
  53 double
  54 SECTION
  55 __ieee754_exp (double x)
  56 {
  57   double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
  58   mynumber junk1, junk2, binexp = {{0, 0}};
  59   int4 i, j, m, n, ex;
  60   double retval;
  61
  62   {
  63     SET_RESTORE_ROUND (FE_TONEAREST);
  64
  65     junk1.x = x;
  66     m = junk1.i[HIGH_HALF];
  67     n = m & hugeint;
  68
  69     if (n > smallint && n < bigint)
  70       {
  71         y = x * log2e.x + three51.x;
  72         bexp = y - three51.x;   /*  multiply the result by 2**bexp        */
  73
  74         junk1.x = y;
  75
  76         eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps               */
  77         t = x - bexp * ln_two1.x;
  78
  79         y = t + three33.x;
  80         base = y - three33.x;   /* t rounded to a multiple of 2**-18      */
  81         junk2.x = y;
  82         del = (t - base) - eps; /*  x = bexp*ln(2) + base + del           */
  83         eps = del + del * del * (p3.x * del + p2.x);
  84
  85         binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
  86
  87         i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
  88         j = (junk2.i[LOW_HALF] & 511) << 1;
  89
  90         al = coar.x[i] * fine.x[j];
  91         bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
  92                + coar.x[i + 1] * fine.x[j + 1]);
  93
  94         rem = (bet + bet * eps) + al * eps;
  95         res = al + rem;
  96         cor = (al - res) + rem;
  97         if (res == (res + cor * err_0))
  98           {
  99             retval = res * binexp.x;
 100             goto ret;
 101           }
 102         else
 103           {
 104             retval = __slowexp (x);
 105             goto ret;
 106           }                     /*if error is over bound */
 107       }
 108
 109     if (n <= smallint)
 110       {
 111         retval = 1.0;
 112         goto ret;
 113       }
 114
 115     if (n >= badint)
 116       {
 117         if (n > infint)
 118           {
 119             retval = x + x;
 120             goto ret;
 121           }                     /* x is NaN */
 122         if (n < infint)
 123           {
 124             if (x > 0)
 125               goto ret_huge;
 126             else
 127               goto ret_tiny;
 128           }
 129         /* x is finite,  cause either overflow or underflow  */
 130         if (junk1.i[LOW_HALF] != 0)
 131           {
 132             retval = x + x;
 133             goto ret;
 134           }                     /*  x is NaN  */
 135         retval = (x > 0) ? inf.x : zero;        /* |x| = inf;  return either inf or 0 */
 136         goto ret;
 137       }
 138
 139     y = x * log2e.x + three51.x;
 140     bexp = y - three51.x;
 141     junk1.x = y;
 142     eps = bexp * ln_two2.x;
 143     t = x - bexp * ln_two1.x;
 144     y = t + three33.x;
 145     base = y - three33.x;
 146     junk2.x = y;
 147     del = (t - base) - eps;
 148     eps = del + del * del * (p3.x * del + p2.x);
 149     i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
 150     j = (junk2.i[LOW_HALF] & 511) << 1;
 151     al = coar.x[i] * fine.x[j];
 152     bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
 153            + coar.x[i + 1] * fine.x[j + 1]);
 154     rem = (bet + bet * eps) + al * eps;
 155     res = al + rem;
 156     cor = (al - res) + rem;
 157     if (m >> 31)
 158       {
 159         ex = junk1.i[LOW_HALF];
 160         if (res < 1.0)
 161           {
 162             res += res;
 163             cor += cor;
 164             ex -= 1;
 165           }
 166         if (ex >= -1022)
 167           {
 168             binexp.i[HIGH_HALF] = (1023 + ex) << 20;
 169             if (res == (res + cor * err_0))
 170               {
 171                 retval = res * binexp.x;
 172                 goto ret;
 173               }
 174             else
 175               {
 176                 retval = __slowexp (x);
 177                 goto check_uflow_ret;
 178               }                 /*if error is over bound */
 179           }
 180         ex = -(1022 + ex);
 181         binexp.i[HIGH_HALF] = (1023 - ex) << 20;
 182         res *= binexp.x;
 183         cor *= binexp.x;
 184         eps = 1.0000000001 + err_0 * binexp.x;
 185         t = 1.0 + res;
 186         y = ((1.0 - t) + res) + cor;
 187         res = t + y;
 188         cor = (t - res) + y;
 189         if (res == (res + eps * cor))
 190           {
 191             binexp.i[HIGH_HALF] = 0x00100000;
 192             retval = (res - 1.0) * binexp.x;
 193             goto check_uflow_ret;
 194           }
 195         else
 196           {
 197             retval = __slowexp (x);
 198             goto check_uflow_ret;
 199           }                     /*   if error is over bound    */
 200       check_uflow_ret:
 201         if (retval < DBL_MIN)
 202           {
 203 #if FLT_EVAL_METHOD != 0
 204             volatile
 205 #endif
 206               double force_underflow = tiny * tiny;
 207             math_force_eval (force_underflow);
 208           }
 209         if (retval == 0)
 210           goto ret_tiny;
 211         goto ret;
 212       }
 213     else
 214       {
 215         binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
 216         if (res == (res + cor * err_0))
 217           retval = res * binexp.x * t256.x;
 218         else
 219           retval = __slowexp (x);
 220         if (__isinf (retval))
 221           goto ret_huge;
 222         else
 223           goto ret;
 224       }
 225   }
 226 ret:
 227   return retval;
 228
 229  ret_huge:
 230   return hhuge * hhuge;
 231
 232  ret_tiny:
 233   return tiny * tiny;
 234 }
 235 #ifndef __ieee754_exp
 236 strong_alias (__ieee754_exp, __exp_finite)
 237 #endif
 238
 239 /* Compute e^(x+xx).  The routine also receives bound of error of previous
 240    calculation.  If after computing exp the error exceeds the allowed bounds,
 241    the routine returns a non-positive number.  Otherwise it returns the
 242    computed result, which is always positive.  */
 243 double
 244 SECTION
 245 __exp1 (double x, double xx, double error)
 246 {
 247   double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
 248   mynumber junk1, junk2, binexp = {{0, 0}};
 249   int4 i, j, m, n, ex;
 250
 251   junk1.x = x;
 252   m = junk1.i[HIGH_HALF];
 253   n = m & hugeint;              /* no sign */
 254
 255   if (n > smallint && n < bigint)
 256     {
 257       y = x * log2e.x + three51.x;
 258       bexp = y - three51.x;     /*  multiply the result by 2**bexp        */
 259
 260       junk1.x = y;
 261
 262       eps = bexp * ln_two2.x;   /* x = bexp*ln(2) + t - eps               */
 263       t = x - bexp * ln_two1.x;
 264
 265       y = t + three33.x;
 266       base = y - three33.x;     /* t rounded to a multiple of 2**-18      */
 267       junk2.x = y;
 268       del = (t - base) + (xx - eps);    /*  x = bexp*ln(2) + base + del      */
 269       eps = del + del * del * (p3.x * del + p2.x);
 270
 271       binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
 272
 273       i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
 274       j = (junk2.i[LOW_HALF] & 511) << 1;
 275
 276       al = coar.x[i] * fine.x[j];
 277       bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
 278              + coar.x[i + 1] * fine.x[j + 1]);
 279
 280       rem = (bet + bet * eps) + al * eps;
 281       res = al + rem;
 282       cor = (al - res) + rem;
 283       if (res == (res + cor * (1.0 + error + err_1)))
 284         return res * binexp.x;
 285       else
 286         return -10.0;
 287     }
 288
 289   if (n <= smallint)
 290     return 1.0;                 /*  if x->0 e^x=1 */
 291
 292   if (n >= badint)
 293     {
 294       if (n > infint)
 295         return (zero / zero);   /* x is NaN,  return invalid */
 296       if (n < infint)
 297         return ((x > 0) ? (hhuge * hhuge) : (tiny * tiny));
 298       /* x is finite,  cause either overflow or underflow  */
 299       if (junk1.i[LOW_HALF] != 0)
 300         return (zero / zero);   /*  x is NaN  */
 301       return ((x > 0) ? inf.x : zero);  /* |x| = inf;  return either inf or 0 */
 302     }
 303
 304   y = x * log2e.x + three51.x;
 305   bexp = y - three51.x;
 306   junk1.x = y;
 307   eps = bexp * ln_two2.x;
 308   t = x - bexp * ln_two1.x;
 309   y = t + three33.x;
 310   base = y - three33.x;
 311   junk2.x = y;
 312   del = (t - base) + (xx - eps);
 313   eps = del + del * del * (p3.x * del + p2.x);
 314   i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
 315   j = (junk2.i[LOW_HALF] & 511) << 1;
 316   al = coar.x[i] * fine.x[j];
 317   bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
 318          + coar.x[i + 1] * fine.x[j + 1]);
 319   rem = (bet + bet * eps) + al * eps;
 320   res = al + rem;
 321   cor = (al - res) + rem;
 322   if (m >> 31)
 323     {
 324       ex = junk1.i[LOW_HALF];
 325       if (res < 1.0)
 326         {
 327           res += res;
 328           cor += cor;
 329           ex -= 1;
 330         }
 331       if (ex >= -1022)
 332         {
 333           binexp.i[HIGH_HALF] = (1023 + ex) << 20;
 334           if (res == (res + cor * (1.0 + error + err_1)))
 335             return res * binexp.x;
 336           else
 337             return -10.0;
 338         }
 339       ex = -(1022 + ex);
 340       binexp.i[HIGH_HALF] = (1023 - ex) << 20;
 341       res *= binexp.x;
 342       cor *= binexp.x;
 343       eps = 1.00000000001 + (error + err_1) * binexp.x;
 344       t = 1.0 + res;
 345       y = ((1.0 - t) + res) + cor;
 346       res = t + y;
 347       cor = (t - res) + y;
 348       if (res == (res + eps * cor))
 349         {
 350           binexp.i[HIGH_HALF] = 0x00100000;
 351           return (res - 1.0) * binexp.x;
 352         }
 353       else
 354         return -10.0;
 355     }
 356   else
 357     {
 358       binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
 359       if (res == (res + cor * (1.0 + error + err_1)))
 360         return res * binexp.x * t256.x;
 361       else
 362         return -10.0;
 363     }
 364 }