sysdeps/ieee754/dbl-64/e_exp.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2016 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             double force_underflow = tiny * tiny;
 204             math_force_eval (force_underflow);
 205           }
 206         if (retval == 0)
 207           goto ret_tiny;
 208         goto ret;
 209       }
 210     else
 211       {
 212         binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
 213         if (res == (res + cor * err_0))
 214           retval = res * binexp.x * t256.x;
 215         else
 216           retval = __slowexp (x);
 217         if (isinf (retval))
 218           goto ret_huge;
 219         else
 220           goto ret;
 221       }
 222   }
 223 ret:
 224   return retval;
 225
 226  ret_huge:
 227   return hhuge * hhuge;
 228
 229  ret_tiny:
 230   return tiny * tiny;
 231 }
 232 #ifndef __ieee754_exp
 233 strong_alias (__ieee754_exp, __exp_finite)
 234 #endif
 235
 236 /* Compute e^(x+xx).  The routine also receives bound of error of previous
 237    calculation.  If after computing exp the error exceeds the allowed bounds,
 238    the routine returns a non-positive number.  Otherwise it returns the
 239    computed result, which is always positive.  */
 240 double
 241 SECTION
 242 __exp1 (double x, double xx, double error)
 243 {
 244   double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
 245   mynumber junk1, junk2, binexp = {{0, 0}};
 246   int4 i, j, m, n, ex;
 247
 248   junk1.x = x;
 249   m = junk1.i[HIGH_HALF];
 250   n = m & hugeint;              /* no sign */
 251
 252   if (n > smallint && n < bigint)
 253     {
 254       y = x * log2e.x + three51.x;
 255       bexp = y - three51.x;     /*  multiply the result by 2**bexp        */
 256
 257       junk1.x = y;
 258
 259       eps = bexp * ln_two2.x;   /* x = bexp*ln(2) + t - eps               */
 260       t = x - bexp * ln_two1.x;
 261
 262       y = t + three33.x;
 263       base = y - three33.x;     /* t rounded to a multiple of 2**-18      */
 264       junk2.x = y;
 265       del = (t - base) + (xx - eps);    /*  x = bexp*ln(2) + base + del      */
 266       eps = del + del * del * (p3.x * del + p2.x);
 267
 268       binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
 269
 270       i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
 271       j = (junk2.i[LOW_HALF] & 511) << 1;
 272
 273       al = coar.x[i] * fine.x[j];
 274       bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
 275              + coar.x[i + 1] * fine.x[j + 1]);
 276
 277       rem = (bet + bet * eps) + al * eps;
 278       res = al + rem;
 279       cor = (al - res) + rem;
 280       if (res == (res + cor * (1.0 + error + err_1)))
 281         return res * binexp.x;
 282       else
 283         return -10.0;
 284     }
 285
 286   if (n <= smallint)
 287     return 1.0;                 /*  if x->0 e^x=1 */
 288
 289   if (n >= badint)
 290     {
 291       if (n > infint)
 292         return (zero / zero);   /* x is NaN,  return invalid */
 293       if (n < infint)
 294         return ((x > 0) ? (hhuge * hhuge) : (tiny * tiny));
 295       /* x is finite,  cause either overflow or underflow  */
 296       if (junk1.i[LOW_HALF] != 0)
 297         return (zero / zero);   /*  x is NaN  */
 298       return ((x > 0) ? inf.x : zero);  /* |x| = inf;  return either inf or 0 */
 299     }
 300
 301   y = x * log2e.x + three51.x;
 302   bexp = y - three51.x;
 303   junk1.x = y;
 304   eps = bexp * ln_two2.x;
 305   t = x - bexp * ln_two1.x;
 306   y = t + three33.x;
 307   base = y - three33.x;
 308   junk2.x = y;
 309   del = (t - base) + (xx - eps);
 310   eps = del + del * del * (p3.x * del + p2.x);
 311   i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
 312   j = (junk2.i[LOW_HALF] & 511) << 1;
 313   al = coar.x[i] * fine.x[j];
 314   bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
 315          + coar.x[i + 1] * fine.x[j + 1]);
 316   rem = (bet + bet * eps) + al * eps;
 317   res = al + rem;
 318   cor = (al - res) + rem;
 319   if (m >> 31)
 320     {
 321       ex = junk1.i[LOW_HALF];
 322       if (res < 1.0)
 323         {
 324           res += res;
 325           cor += cor;
 326           ex -= 1;
 327         }
 328       if (ex >= -1022)
 329         {
 330           binexp.i[HIGH_HALF] = (1023 + ex) << 20;
 331           if (res == (res + cor * (1.0 + error + err_1)))
 332             return res * binexp.x;
 333           else
 334             return -10.0;
 335         }
 336       ex = -(1022 + ex);
 337       binexp.i[HIGH_HALF] = (1023 - ex) << 20;
 338       res *= binexp.x;
 339       cor *= binexp.x;
 340       eps = 1.00000000001 + (error + err_1) * binexp.x;
 341       t = 1.0 + res;
 342       y = ((1.0 - t) + res) + cor;
 343       res = t + y;
 344       cor = (t - res) + y;
 345       if (res == (res + eps * cor))
 346         {
 347           binexp.i[HIGH_HALF] = 0x00100000;
 348           return (res - 1.0) * binexp.x;
 349         }
 350       else
 351         return -10.0;
 352     }
 353   else
 354     {
 355       binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
 356       if (res == (res + cor * (1.0 + error + err_1)))
 357         return res * binexp.x * t256.x;
 358       else
 359         return -10.0;
 360     }
 361 }