sysdeps/ieee754/flt-32/e_expf.c

   1 /* Single-precision floating point e^x.
   2    Copyright (C) 1997-2014 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Contributed by Geoffrey Keating <geoffk@ozemail.com.au>
   5
   6    The GNU C Library is free software; you can redistribute it and/or
   7    modify it under the terms of the GNU Lesser General Public
   8    License as published by the Free Software Foundation; either
   9    version 2.1 of the License, or (at your option) any later version.
  10
  11    The GNU C Library 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 GNU
  14    Lesser General Public License for more details.
  15
  16    You should have received a copy of the GNU Lesser General Public
  17    License along with the GNU C Library; if not, see
  18    <http://www.gnu.org/licenses/>.  */
  19
  20 /* How this works:
  21
  22    The input value, x, is written as
  23
  24    x = n * ln(2) + t/512 + delta[t] + x;
  25
  26    where:
  27    - n is an integer, 127 >= n >= -150;
  28    - t is an integer, 177 >= t >= -177
  29    - delta is based on a table entry, delta[t] < 2^-28
  30    - x is whatever is left, |x| < 2^-10
  31
  32    Then e^x is approximated as
  33
  34    e^x = 2^n ( e^(t/512 + delta[t])
  35                + ( e^(t/512 + delta[t])
  36                    * ( p(x + delta[t] + n * ln(2)) - delta ) ) )
  37
  38    where
  39    - p(x) is a polynomial approximating e(x)-1;
  40    - e^(t/512 + delta[t]) is obtained from a table.
  41
  42    The table used is the same one as for the double precision version;
  43    since we have the table, we might as well use it.
  44
  45    It turns out to be faster to do calculations in double precision than
  46    to perform an 'accurate table method' expf, because of the range reduction
  47    overhead (compare exp2f).
  48    */
  49 #include <float.h>
  50 #include <ieee754.h>
  51 #include <math.h>
  52 #include <fenv.h>
  53 #include <inttypes.h>
  54 #include <math_private.h>
  55
  56 extern const float __exp_deltatable[178];
  57 extern const double __exp_atable[355] /* __attribute__((mode(DF))) */;
  58
  59 static const float TWOM100 = 7.88860905e-31;
  60 static const float TWO127 = 1.7014118346e+38;
  61
  62 float
  63 __ieee754_expf (float x)
  64 {
  65   static const float himark = 88.72283935546875;
  66   static const float lomark = -103.972084045410;
  67   /* Check for usual case.  */
  68   if (isless (x, himark) && isgreater (x, lomark))
  69     {
  70       static const float THREEp42 = 13194139533312.0;
  71       static const float THREEp22 = 12582912.0;
  72       /* 1/ln(2).  */
  73 #undef M_1_LN2
  74       static const float M_1_LN2 = 1.44269502163f;
  75       /* ln(2) */
  76 #undef M_LN2
  77       static const double M_LN2 = .6931471805599452862;
  78
  79       int tval;
  80       double x22, t, result, dx;
  81       float n, delta;
  82       union ieee754_double ex2_u;
  83
  84       {
  85         SET_RESTORE_ROUND_NOEXF (FE_TONEAREST);
  86
  87         /* Calculate n.  */
  88         n = x * M_1_LN2 + THREEp22;
  89         n -= THREEp22;
  90         dx = x - n*M_LN2;
  91
  92         /* Calculate t/512.  */
  93         t = dx + THREEp42;
  94         t -= THREEp42;
  95         dx -= t;
  96
  97         /* Compute tval = t.  */
  98         tval = (int) (t * 512.0);
  99
 100         if (t >= 0)
 101           delta = - __exp_deltatable[tval];
 102         else
 103           delta = __exp_deltatable[-tval];
 104
 105         /* Compute ex2 = 2^n e^(t/512+delta[t]).  */
 106         ex2_u.d = __exp_atable[tval+177];
 107         ex2_u.ieee.exponent += (int) n;
 108
 109         /* Approximate e^(dx+delta) - 1, using a second-degree polynomial,
 110            with maximum error in [-2^-10-2^-28,2^-10+2^-28]
 111            less than 5e-11.  */
 112         x22 = (0.5000000496709180453 * dx + 1.0000001192102037084) * dx + delta;
 113       }
 114
 115       /* Return result.  */
 116       result = x22 * ex2_u.d + ex2_u.d;
 117       return (float) result;
 118     }
 119   /* Exceptional cases:  */
 120   else if (isless (x, himark))
 121     {
 122       if (__isinff (x))
 123         /* e^-inf == 0, with no error.  */
 124         return 0;
 125       else
 126         /* Underflow */
 127         return TWOM100 * TWOM100;
 128     }
 129   else
 130     /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
 131     return TWO127*x;
 132 }
 133 strong_alias (__ieee754_expf, __expf_finite)