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