newlib/libm/mathfp/s_exp.c

   1
   2 /* @(#)z_exp.c 1.0 98/08/13 */
   3 /******************************************************************
   4  * The following routines are coded directly from the algorithms
   5  * and coefficients given in "Software Manual for the Elementary
   6  * Functions" by William J. Cody, Jr. and William Waite, Prentice
   7  * Hall, 1980.
   8  ******************************************************************/
   9
  10 /*
  11 FUNCTION
  12         <<exp>>, <<expf>>---exponential
  13 INDEX
  14         exp
  15 INDEX
  16         expf
  17
  18 ANSI_SYNOPSIS
  19         #include <math.h>
  20         double exp(double <[x]>);
  21         float expf(float <[x]>);
  22
  23 TRAD_SYNOPSIS
  24         #include <math.h>
  25         double exp(<[x]>);
  26         double <[x]>;
  27
  28         float expf(<[x]>);
  29         float <[x]>;
  30
  31 DESCRIPTION
  32         <<exp>> and <<expf>> calculate the exponential of <[x]>, that is,
  33         @ifinfo
  34         e raised to the power <[x]> (where e
  35         @end ifinfo
  36         @tex
  37         $e^x$ (where $e$
  38         @end tex
  39         is the base of the natural system of logarithms, approximately 2.71828).
  40
  41 RETURNS
  42         On success, <<exp>> and <<expf>> return the calculated value.
  43         If the result underflows, the returned value is <<0>>.  If the
  44         result overflows, the returned value is <<HUGE_VAL>>.  In
  45         either case, <<errno>> is set to <<ERANGE>>.
  46
  47 PORTABILITY
  48         <<exp>> is ANSI C.  <<expf>> is an extension.
  49
  50 */
  51
  52 /*****************************************************************
  53  * Exponential Function
  54  *
  55  * Input:
  56  *   x - floating point value
  57  *
  58  * Output:
  59  *   e raised to x.
  60  *
  61  * Description:
  62  *   This routine returns e raised to the xth power.
  63  *
  64  *****************************************************************/
  65
  66 #include <float.h>
  67 #include "fdlibm.h"
  68 #include "zmath.h"
  69
  70 #ifndef _DOUBLE_IS_32BITS
  71
  72 static const double INV_LN2 = 1.4426950408889634074;
  73 static const double LN2 = 0.6931471805599453094172321;
  74 static const double p[] = { 0.25, 0.75753180159422776666e-2,
  75                      0.31555192765684646356e-4 };
  76 static const double q[] = { 0.5, 0.56817302698551221787e-1,
  77                      0.63121894374398504557e-3,
  78                      0.75104028399870046114e-6 };
  79
  80 double
  81 _DEFUN (exp, (double),
  82         double x)
  83 {
  84   int N;
  85   double g, z, R, P, Q;
  86
  87   switch (numtest (x))
  88     {
  89       case NAN:
  90         errno = EDOM;
  91         return (x);
  92       case INF:
  93         errno = ERANGE;
  94         if (ispos (x))
  95           return (z_infinity.d);
  96         else
  97           return (0.0);
  98       case 0:
  99         return (1.0);
 100     }
 101
 102   /* Check for out of bounds. */
 103   if (x > BIGX || x < SMALLX)
 104     {
 105       errno = ERANGE;
 106       return (x);
 107     }
 108
 109   /* Check for a value too small to calculate. */
 110   if (-z_rooteps < x && x < z_rooteps)
 111     {
 112       return (1.0);
 113     }
 114
 115   /* Calculate the exponent. */
 116   if (x < 0.0)
 117     N = (int) (x * INV_LN2 - 0.5);
 118   else
 119     N = (int) (x * INV_LN2 + 0.5);
 120
 121   /* Construct the mantissa. */
 122   g = x - N * LN2;
 123   z = g * g;
 124   P = g * ((p[2] * z + p[1]) * z + p[0]);
 125   Q = ((q[3] * z + q[2]) * z + q[1]) * z + q[0];
 126   R = 0.5 + P / (Q - P);
 127
 128   /* Return the floating point value. */
 129   N++;
 130   return (ldexp (R, N));
 131 }
 132
 133 #endif /* _DOUBLE_IS_32BITS */