libm/float/expf.c

   1 /*                                                      expf.c
   2  *
   3  *      Exponential function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, expf();
  10  *
  11  * y = expf( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns e (2.71828...) raised to the x power.
  18  *
  19  * Range reduction is accomplished by separating the argument
  20  * into an integer k and fraction f such that
  21  *
  22  *     x    k  f
  23  *    e  = 2  e.
  24  *
  25  * A polynomial is used to approximate exp(f)
  26  * in the basic range [-0.5, 0.5].
  27  *
  28  *
  29  * ACCURACY:
  30  *
  31  *                      Relative error:
  32  * arithmetic   domain     # trials      peak         rms
  33  *    IEEE      +- MAXLOG   100000      1.7e-7      2.8e-8
  34  *
  35  *
  36  * Error amplification in the exponential function can be
  37  * a serious matter.  The error propagation involves
  38  * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
  39  * which shows that a 1 lsb error in representing X produces
  40  * a relative error of X times 1 lsb in the function.
  41  * While the routine gives an accurate result for arguments
  42  * that are exactly represented by a double precision
  43  * computer number, the result contains amplified roundoff
  44  * error for large arguments not exactly represented.
  45  *
  46  *
  47  * ERROR MESSAGES:
  48  *
  49  *   message         condition      value returned
  50  * expf underflow    x < MINLOGF         0.0
  51  * expf overflow     x > MAXLOGF         MAXNUMF
  52  *
  53  */
  54 \f
  55 /*
  56 Cephes Math Library Release 2.2:  June, 1992
  57 Copyright 1984, 1987, 1989 by Stephen L. Moshier
  58 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  59 */
  60
  61 /* Single precision exponential function.
  62  * test interval: [-0.5, +0.5]
  63  * trials: 80000
  64  * peak relative error: 7.6e-8
  65  * rms relative error: 2.8e-8
  66  */
  67 #include <math.h>
  68 extern float LOG2EF, MAXLOGF, MINLOGF, MAXNUMF;
  69
  70 static float C1 =   0.693359375;
  71 static float C2 =  -2.12194440e-4;
  72
  73
  74
  75 float floorf( float ), ldexpf( float, int );
  76
  77 float expf( float xx )
  78 {
  79 float x, z;
  80 int n;
  81
  82 x = xx;
  83
  84
  85 if( x > MAXLOGF)
  86         {
  87         mtherr( "expf", OVERFLOW );
  88         return( MAXNUMF );
  89         }
  90
  91 if( x < MINLOGF )
  92         {
  93         mtherr( "expf", UNDERFLOW );
  94         return(0.0);
  95         }
  96
  97 /* Express e**x = e**g 2**n
  98  *   = e**g e**( n loge(2) )
  99  *   = e**( g + n loge(2) )
 100  */
 101 z = floorf( LOG2EF * x + 0.5 ); /* floor() truncates toward -infinity. */
 102 x -= z * C1;
 103 x -= z * C2;
 104 n = z;
 105
 106 z = x * x;
 107 /* Theoretical peak relative error in [-0.5, +0.5] is 4.2e-9. */
 108 z =
 109 ((((( 1.9875691500E-4  * x
 110    + 1.3981999507E-3) * x
 111    + 8.3334519073E-3) * x
 112    + 4.1665795894E-2) * x
 113    + 1.6666665459E-1) * x
 114    + 5.0000001201E-1) * z
 115    + x
 116    + 1.0;
 117
 118 /* multiply by power of 2 */
 119 x = ldexpf( z, n );
 120
 121 return( x );
 122 }