1 files changed, 92 insertions, 0 deletions
diff --git a/newlib/libm/mathfp/sf_exp.c b/newlib/libm/mathfp/sf_exp.c
new file mode 100644
index 00000000000..45d9c31fc44
--- /dev/null
+++ b/newlib/libm/mathfp/sf_exp.c
@@ -0,0 +1,92 @@
+
+/* @(#)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_32_BITS
+
+double exp (double x)
+{
+  return (double) expf ((float) x);
+}
+
+#endif /* _DOUBLE_IS_32_BITS */