1 files changed, 133 insertions, 0 deletions
diff --git a/newlib/libm/mathfp/s_exp.c b/newlib/libm/mathfp/s_exp.c
new file mode 100644
index 00000000000..8c7f723fe13
--- /dev/null
+++ b/newlib/libm/mathfp/s_exp.c
@@ -0,0 +1,133 @@
+
+/* @(#)z_exp.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.
+ ******************************************************************/
+
+/*
+FUNCTION
+        <<exp>>, <<expf>>---exponential
+INDEX
+        exp
+INDEX
+        expf
+
+ANSI_SYNOPSIS
+        #include <math.h>
+        double exp(double <[x]>);
+        float expf(float <[x]>);
+
+TRAD_SYNOPSIS
+        #include <math.h>
+        double exp(<[x]>);
+        double <[x]>;
+
+        float expf(<[x]>);
+        float <[x]>;
+
+DESCRIPTION
+        <<exp>> and <<expf>> calculate the exponential of <[x]>, that is,
+        @ifinfo
+        e raised to the power <[x]> (where e
+        @end ifinfo
+        @tex
+        $e^x$ (where $e$
+        @end tex
+        is the base of the natural system of logarithms, approximately 2.71828).
+
+RETURNS
+        On success, <<exp>> and <<expf>> return the calculated value.
+        If the result underflows, the returned value is <<0>>.  If the
+        result overflows, the returned value is <<HUGE_VAL>>.  In
+        either case, <<errno>> is set to <<ERANGE>>.
+
+PORTABILITY
+        <<exp>> is ANSI C.  <<expf>> is an extension.
+
+*/
+
+/*****************************************************************
+ * 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"
+
+#ifndef _DOUBLE_IS_32BITS
+
+static const double INV_LN2 = 1.4426950408889634074;
+static const double LN2 = 0.6931471805599453094172321;
+static const double p[] = { 0.25, 0.75753180159422776666e-2,
+                     0.31555192765684646356e-4 };
+static const double q[] = { 0.5, 0.56817302698551221787e-1,
+                     0.63121894374398504557e-3,
+                     0.75104028399870046114e-6 };
+
+double
+_DEFUN (exp, (double),
+        double x)
+{
+  int N;
+  double g, z, R, P, Q;
+
+  switch (numtest (x))
+    {
+      case NAN:
+        errno = EDOM;
+        return (x);
+      case INF:
+        errno = ERANGE;
+        if (ispos (x))
+          return (z_infinity.d);
+        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 < x && x < z_rooteps)
+    {
+      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[2] * z + p[1]) * z + p[0]);
+  Q = ((q[3] * z + q[2]) * z + q[1]) * z + q[0];
+  R = 0.5 + P / (Q - P);
+
+  /* Return the floating point value. */
+  N++;
+  return (ldexp (R, N));
+}
+
+#endif /* _DOUBLE_IS_32BITS */