newlib/libm/mathfp/sf_exp.c

   1
   2 /* @(#)z_expf.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  * Exponential Function
  11  *
  12  * Input:
  13  *   x - floating point value
  14  *
  15  * Output:
  16  *   e raised to x.
  17  *
  18  * Description:
  19  *   This routine returns e raised to the xth power.
  20  *
  21  *****************************************************************/
  22
  23 #include <float.h>
  24 #include "fdlibm.h"
  25 #include "zmath.h"
  26
  27 static const float INV_LN2 = 1.442695040;
  28 static const float LN2 = 0.693147180;
  29 static const float p[] = { 0.249999999950, 0.00416028863 };
  30 static const float q[] = { 0.5, 0.04998717878 };
  31
  32 float
  33 expf (float x)
  34 {
  35   int N;
  36   float g, z, R, P, Q;
  37
  38   switch (numtestf (x))
  39     {
  40       case NAN:
  41         errno = EDOM;
  42         return (x);
  43       case INF:
  44         errno = ERANGE;
  45         if (isposf (x))
  46           return (z_infinity_f.f);
  47         else
  48           return (0.0);
  49       case 0:
  50         return (1.0);
  51     }
  52
  53   /* Check for out of bounds. */
  54   if (x > BIGX || x < SMALLX)
  55     {
  56       errno = ERANGE;
  57       return (x);
  58     }
  59
  60   /* Check for a value too small to calculate. */
  61   if (-z_rooteps_f < x && x < z_rooteps_f)
  62     {
  63       return (1.0);
  64     }
  65
  66   /* Calculate the exponent. */
  67   if (x < 0.0)
  68     N = (int) (x * INV_LN2 - 0.5);
  69   else
  70     N = (int) (x * INV_LN2 + 0.5);
  71
  72   /* Construct the mantissa. */
  73   g = x - N * LN2;
  74   z = g * g;
  75   P = g * (p[1] * z + p[0]);
  76   Q = q[1] * z + q[0];
  77   R = 0.5 + P / (Q - P);
  78
  79   /* Return the floating point value. */
  80   N++;
  81   return (ldexpf (R, N));
  82 }
  83
  84 #ifdef _DOUBLE_IS_32BITS
  85
  86 double exp (double x)
  87 {
  88   return (double) expf ((float) x);
  89 }
  90
  91 #endif /* _DOUBLE_IS_32BITS */