toolchain/gcc/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 _DEFUN (expf, (float),
  34         float x)
  35 {
  36   int N;
  37   float g, z, R, P, Q;
  38
  39   switch (numtestf (x))
  40     {
  41       case NAN:
  42         errno = EDOM;
  43         return (x);
  44       case INF:
  45         errno = ERANGE;
  46         if (isposf (x))
  47           return (z_infinity_f.f);
  48         else
  49           return (0.0);
  50       case 0:
  51         return (1.0);
  52     }
  53
  54   /* Check for out of bounds. */
  55   if (x > BIGX || x < SMALLX)
  56     {
  57       errno = ERANGE;
  58       return (x);
  59     }
  60
  61   /* Check for a value too small to calculate. */
  62   if (-z_rooteps_f < x && x < z_rooteps_f)
  63     {
  64       return (1.0);
  65     }
  66
  67   /* Calculate the exponent. */
  68   if (x < 0.0)
  69     N = (int) (x * INV_LN2 - 0.5);
  70   else
  71     N = (int) (x * INV_LN2 + 0.5);
  72
  73   /* Construct the mantissa. */
  74   g = x - N * LN2;
  75   z = g * g;
  76   P = g * (p[1] * z + p[0]);
  77   Q = q[1] * z + q[0];
  78   R = 0.5 + P / (Q - P);
  79
  80   /* Return the floating point value. */
  81   N++;
  82   return (ldexpf (R, N));
  83 }
  84
  85 #ifdef _DOUBLE_IS_32BITS
  86
  87 double exp (double x)
  88 {
  89   return (double) expf ((float) x);
  90 }
  91
  92 #endif /* _DOUBLE_IS_32BITS */