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