Back to index

glibc  2.9
e_expf.c
Go to the documentation of this file.
00001 /* Single-precision floating point e^x.
00002    Copyright (C) 1997, 1998, 2005, 2006 Free Software Foundation, Inc.
00003    This file is part of the GNU C Library.
00004    Contributed by Geoffrey Keating <geoffk@ozemail.com.au>
00005 
00006    The GNU C Library is free software; you can redistribute it and/or
00007    modify it under the terms of the GNU Lesser General Public
00008    License as published by the Free Software Foundation; either
00009    version 2.1 of the License, or (at your option) any later version.
00010 
00011    The GNU C Library is distributed in the hope that it will be useful,
00012    but WITHOUT ANY WARRANTY; without even the implied warranty of
00013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014    Lesser General Public License for more details.
00015 
00016    You should have received a copy of the GNU Lesser General Public
00017    License along with the GNU C Library; if not, write to the Free
00018    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00019    02111-1307 USA.  */
00020 
00021 /* How this works:
00022 
00023    The input value, x, is written as
00024 
00025    x = n * ln(2) + t/512 + delta[t] + x;
00026 
00027    where:
00028    - n is an integer, 127 >= n >= -150;
00029    - t is an integer, 177 >= t >= -177
00030    - delta is based on a table entry, delta[t] < 2^-28
00031    - x is whatever is left, |x| < 2^-10
00032 
00033    Then e^x is approximated as
00034 
00035    e^x = 2^n ( e^(t/512 + delta[t])
00036                + ( e^(t/512 + delta[t])
00037                    * ( p(x + delta[t] + n * ln(2)) - delta ) ) )
00038 
00039    where
00040    - p(x) is a polynomial approximating e(x)-1;
00041    - e^(t/512 + delta[t]) is obtained from a table.
00042 
00043    The table used is the same one as for the double precision version;
00044    since we have the table, we might as well use it.
00045 
00046    It turns out to be faster to do calculations in double precision than
00047    to perform an 'accurate table method' expf, because of the range reduction
00048    overhead (compare exp2f).
00049    */
00050 #ifndef _GNU_SOURCE
00051 #define _GNU_SOURCE
00052 #endif
00053 #include <float.h>
00054 #include <ieee754.h>
00055 #include <math.h>
00056 #include <fenv.h>
00057 #include <inttypes.h>
00058 #include <math_private.h>
00059 
00060 extern const float __exp_deltatable[178];
00061 extern const double __exp_atable[355] /* __attribute__((mode(DF))) */;
00062 
00063 static const volatile float TWOM100 = 7.88860905e-31;
00064 static const volatile float TWO127 = 1.7014118346e+38;
00065 
00066 float
00067 __ieee754_expf (float x)
00068 {
00069   static const float himark = 88.72283935546875;
00070   static const float lomark = -103.972084045410;
00071   /* Check for usual case.  */
00072   if (isless (x, himark) && isgreater (x, lomark))
00073     {
00074       static const float THREEp42 = 13194139533312.0;
00075       static const float THREEp22 = 12582912.0;
00076       /* 1/ln(2).  */
00077 #undef M_1_LN2
00078       static const float M_1_LN2 = 1.44269502163f;
00079       /* ln(2) */
00080 #undef M_LN2
00081       static const double M_LN2 = .6931471805599452862;
00082 
00083       int tval;
00084       double x22, t, result, dx;
00085       float n, delta;
00086       union ieee754_double ex2_u;
00087       fenv_t oldenv;
00088 
00089       feholdexcept (&oldenv);
00090 #ifdef FE_TONEAREST
00091       fesetround (FE_TONEAREST);
00092 #endif
00093 
00094       /* Calculate n.  */
00095       n = x * M_1_LN2 + THREEp22;
00096       n -= THREEp22;
00097       dx = x - n*M_LN2;
00098 
00099       /* Calculate t/512.  */
00100       t = dx + THREEp42;
00101       t -= THREEp42;
00102       dx -= t;
00103 
00104       /* Compute tval = t.  */
00105       tval = (int) (t * 512.0);
00106 
00107       if (t >= 0)
00108        delta = - __exp_deltatable[tval];
00109       else
00110        delta = __exp_deltatable[-tval];
00111 
00112       /* Compute ex2 = 2^n e^(t/512+delta[t]).  */
00113       ex2_u.d = __exp_atable[tval+177];
00114       ex2_u.ieee.exponent += (int) n;
00115 
00116       /* Approximate e^(dx+delta) - 1, using a second-degree polynomial,
00117         with maximum error in [-2^-10-2^-28,2^-10+2^-28]
00118         less than 5e-11.  */
00119       x22 = (0.5000000496709180453 * dx + 1.0000001192102037084) * dx + delta;
00120 
00121       /* Return result.  */
00122       fesetenv (&oldenv);
00123 
00124       result = x22 * ex2_u.d + ex2_u.d;
00125       return (float) result;
00126     }
00127   /* Exceptional cases:  */
00128   else if (isless (x, himark))
00129     {
00130       if (__isinff (x))
00131        /* e^-inf == 0, with no error.  */
00132        return 0;
00133       else
00134        /* Underflow */
00135        return TWOM100 * TWOM100;
00136     }
00137   else
00138     /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
00139     return TWO127*x;
00140 }