Retro68/gcc/newlib/libm/mathfp/sf_exp.c

/* @(#)z_expf.c 1.0 98/08/13 */
/******************************************************************
 * The following routines are coded directly from the algorithms
 * and coefficients given in "Software Manual for the Elementary
 * Functions" by William J. Cody, Jr. and William Waite, Prentice
 * Hall, 1980.
 ******************************************************************/
/******************************************************************
 * Exponential Function
 *
 * Input:
 *   x - floating point value
 *
 * Output:
 *   e raised to x.
 *
 * Description:
 *   This routine returns e raised to the xth power.
 *
 *****************************************************************/

#include <float.h>
#include "fdlibm.h"
#include "zmath.h"

static const float INV_LN2 = 1.442695040;
static const float LN2 = 0.693147180;
static const float p[] = { 0.249999999950, 0.00416028863 };
static const float q[] = { 0.5, 0.04998717878 };

float
_DEFUN (expf, (float),
        float x)
{
  int N;
  float g, z, R, P, Q;

  switch (numtestf (x))
    {
      case NAN:
        errno = EDOM;
        return (x);
      case INF:
        errno = ERANGE;
        if (isposf (x))
          return (z_infinity_f.f);
        else
          return (0.0);
      case 0:
        return (1.0);
    }

  /* Check for out of bounds. */
  if (x > BIGX || x < SMALLX)
    {
      errno = ERANGE;
      return (x);
    }

  /* Check for a value too small to calculate. */
  if (-z_rooteps_f < x && x < z_rooteps_f)
    {
      return (1.0);
    }

  /* Calculate the exponent. */
  if (x < 0.0)
    N = (int) (x * INV_LN2 - 0.5);
  else
    N = (int) (x * INV_LN2 + 0.5);

  /* Construct the mantissa. */
  g = x - N * LN2;
  z = g * g;
  P = g * (p[1] * z + p[0]);
  Q = q[1] * z + q[0];
  R = 0.5 + P / (Q - P);

  /* Return the floating point value. */
  N++;
  return (ldexpf (R, N));
}

#ifdef _DOUBLE_IS_32BITS

double exp (double x)
{
  return (double) expf ((float) x);
}

#endif /* _DOUBLE_IS_32BITS */