1 files changed, 122 insertions, 0 deletions
diff --git a/libm/float/expf.c b/libm/float/expf.c
new file mode 100644
index 000000000..073678b99
--- /dev/null
+++ b/libm/float/expf.c
@@ -0,0 +1,122 @@
+/*							expf.c
+ *
+ *	Exponential function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, expf();
+ *
+ * y = expf( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns e (2.71828...) raised to the x power.
+ *
+ * Range reduction is accomplished by separating the argument
+ * into an integer k and fraction f such that
+ *
+ *     x    k  f
+ *    e  = 2  e.
+ *
+ * A polynomial is used to approximate exp(f)
+ * in the basic range [-0.5, 0.5].
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      +- MAXLOG   100000      1.7e-7      2.8e-8
+ *
+ *
+ * Error amplification in the exponential function can be
+ * a serious matter.  The error propagation involves
+ * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
+ * which shows that a 1 lsb error in representing X produces
+ * a relative error of X times 1 lsb in the function.
+ * While the routine gives an accurate result for arguments
+ * that are exactly represented by a double precision
+ * computer number, the result contains amplified roundoff
+ * error for large arguments not exactly represented.
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * expf underflow    x < MINLOGF         0.0
+ * expf overflow     x > MAXLOGF         MAXNUMF
+ *
+ */
+
+/*
+Cephes Math Library Release 2.2:  June, 1992
+Copyright 1984, 1987, 1989 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+/* Single precision exponential function.
+ * test interval: [-0.5, +0.5]
+ * trials: 80000
+ * peak relative error: 7.6e-8
+ * rms relative error: 2.8e-8
+ */
+#include <math.h>
+extern float LOG2EF, MAXLOGF, MINLOGF, MAXNUMF;
+
+static float C1 =   0.693359375;
+static float C2 =  -2.12194440e-4;
+
+
+
+float floorf( float ), ldexpf( float, int );
+
+float expf( float xx )
+{
+float x, z;
+int n;
+
+x = xx;
+
+
+if( x > MAXLOGF)
+	{
+	mtherr( "expf", OVERFLOW );
+	return( MAXNUMF );
+	}
+
+if( x < MINLOGF )
+	{
+	mtherr( "expf", UNDERFLOW );
+	return(0.0);
+	}
+
+/* Express e**x = e**g 2**n
+ *   = e**g e**( n loge(2) )
+ *   = e**( g + n loge(2) )
+ */
+z = floorf( LOG2EF * x + 0.5 ); /* floor() truncates toward -infinity. */
+x -= z * C1;
+x -= z * C2;
+n = z;
+
+z = x * x;
+/* Theoretical peak relative error in [-0.5, +0.5] is 4.2e-9. */
+z =
+((((( 1.9875691500E-4  * x
+   + 1.3981999507E-3) * x
+   + 8.3334519073E-3) * x
+   + 4.1665795894E-2) * x
+   + 1.6666665459E-1) * x
+   + 5.0000001201E-1) * z
+   + x
+   + 1.0;
+
+/* multiply by power of 2 */
+x = ldexpf( z, n );
+
+return( x );
+}