1/* Single-precision floating point e^x.
2   Copyright (C) 1997, 1998 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, write to the Free
18   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
19   02111-1307 USA.  */
20
21/* How this works:
22
23   The input value, x, is written as
24
25   x = n * ln(2) + t/512 + delta[t] + x;
26
27   where:
28   - n is an integer, 127 >= n >= -150;
29   - t is an integer, 177 >= t >= -177
30   - delta is based on a table entry, delta[t] < 2^-28
31   - x is whatever is left, |x| < 2^-10
32
33   Then e^x is approximated as
34
35   e^x = 2^n ( e^(t/512 + delta[t])
36               + ( e^(t/512 + delta[t])
37                   * ( p(x + delta[t] + n * ln(2)) - delta ) ) )
38
39   where
40   - p(x) is a polynomial approximating e(x)-1;
41   - e^(t/512 + delta[t]) is obtained from a table.
42
43   The table used is the same one as for the double precision version;
44   since we have the table, we might as well use it.
45
46   It turns out to be faster to do calculations in double precision than
47   to perform an 'accurate table method' expf, because of the range reduction
48   overhead (compare exp2f).
49   */
50#ifndef _GNU_SOURCE
51#define _GNU_SOURCE
52#endif
53#include <float.h>
54#include <ieee754.h>
55#include <math.h>
56#include <fenv.h>
57#include <inttypes.h>
58#include <math_private.h>
59
60extern const float __exp_deltatable[178];
61extern const double __exp_atable[355] /* __attribute__((mode(DF))) */;
62
63static const volatile float TWOM100 = 7.88860905e-31;
64static const volatile float TWO127 = 1.7014118346e+38;
65
66float
67__ieee754_expf (float x)
68{
69  static const float himark = 88.72283935546875;
70  static const float lomark = -103.972084045410;
71  /* Check for usual case.  */
72  if (isless (x, himark) && isgreater (x, lomark))
73    {
74      static const float THREEp42 = 13194139533312.0;
75      static const float THREEp22 = 12582912.0;
76      /* 1/ln(2).  */
77#undef M_1_LN2
78      static const float M_1_LN2 = 1.44269502163f;
79      /* ln(2) */
80#undef M_LN2
81      static const double M_LN2 = .6931471805599452862;
82
83      int tval;
84      double x22, t, result, dx;
85      float n, delta;
86      union ieee754_double ex2_u;
87      fenv_t oldenv;
88
89      feholdexcept (&oldenv);
90#ifdef FE_TONEAREST
91      fesetround (FE_TONEAREST);
92#endif
93
94      /* Calculate n.  */
95      n = x * M_1_LN2 + THREEp22;
96      n -= THREEp22;
97      dx = x - n*M_LN2;
98
99      /* Calculate t/512.  */
100      t = dx + THREEp42;
101      t -= THREEp42;
102      dx -= t;
103
104      /* Compute tval = t.  */
105      tval = (int) (t * 512.0);
106
107      if (t >= 0)
108	delta = - __exp_deltatable[tval];
109      else
110	delta = __exp_deltatable[-tval];
111
112      /* Compute ex2 = 2^n e^(t/512+delta[t]).  */
113      ex2_u.d = __exp_atable[tval+177];
114      ex2_u.ieee.exponent += (int) n;
115
116      /* Approximate e^(dx+delta) - 1, using a second-degree polynomial,
117	 with maximum error in [-2^-10-2^-28,2^-10+2^-28]
118	 less than 5e-11.  */
119      x22 = (0.5000000496709180453 * dx + 1.0000001192102037084) * dx + delta;
120
121      /* Return result.  */
122      fesetenv (&oldenv);
123
124      result = x22 * ex2_u.d + ex2_u.d;
125      return (float) result;
126    }
127  /* Exceptional cases:  */
128  else if (isless (x, himark))
129    {
130      if (__isinff (x))
131	/* e^-inf == 0, with no error.  */
132	return 0;
133      else
134	/* Underflow */
135	return TWOM100 * TWOM100;
136    }
137  else
138    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
139    return TWO127*x;
140}
141