1 files changed, 115 insertions, 0 deletions
diff --git a/libm/float/exp10f.c b/libm/float/exp10f.c
new file mode 100644
index 000000000..c7c62c567
--- /dev/null
+++ b/libm/float/exp10f.c
@@ -0,0 +1,115 @@
+/*							exp10f.c
+ *
+ *	Base 10 exponential function
+ *      (Common antilogarithm)
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, exp10f();
+ *
+ * y = exp10f( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns 10 raised to the x power.
+ *
+ * Range reduction is accomplished by expressing the argument
+ * as 10**x = 2**n 10**f, with |f| < 0.5 log10(2).
+ * A polynomial approximates 10**f.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -38,+38     100000      9.8e-8      2.8e-8
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * exp10 underflow    x < -MAXL10        0.0
+ * exp10 overflow     x > MAXL10       MAXNUM
+ *
+ * IEEE single arithmetic: MAXL10 = 38.230809449325611792.
+ *
+ */
+
+/*
+Cephes Math Library Release 2.2:  June, 1992
+Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+
+#include <math.h>
+
+static float P[] = {
+ 2.063216740311022E-001,
+ 5.420251702225484E-001,
+ 1.171292686296281E+000,
+ 2.034649854009453E+000,
+ 2.650948748208892E+000,
+ 2.302585167056758E+000
+};
+
+/*static float LOG102 = 3.01029995663981195214e-1;*/
+static float LOG210 = 3.32192809488736234787e0;
+static float LG102A = 3.00781250000000000000E-1;
+static float LG102B = 2.48745663981195213739E-4;
+static float MAXL10 = 38.230809449325611792;
+
+
+
+
+extern float MAXNUMF;
+
+float floorf(float), ldexpf(float, int), polevlf(float, float *, int);
+
+float exp10f(float xx)
+{
+float x, px, qx;
+short n;
+
+x = xx;
+if( x > MAXL10 )
+	{
+	mtherr( "exp10f", OVERFLOW );
+	return( MAXNUMF );
+	}
+
+if( x < -MAXL10 )	/* Would like to use MINLOG but can't */
+	{
+	mtherr( "exp10f", UNDERFLOW );
+	return(0.0);
+	}
+
+/* The following is necessary because range reduction blows up: */
+if( x == 0 )
+	return(1.0);
+
+/* Express 10**x = 10**g 2**n
+ *   = 10**g 10**( n log10(2) )
+ *   = 10**( g + n log10(2) )
+ */
+px = x * LOG210;
+qx = floorf( px + 0.5 );
+n = qx;
+x -= qx * LG102A;
+x -= qx * LG102B;
+
+/* rational approximation for exponential
+ * of the fractional part:
+ * 10**x - 1  =  2x P(x**2)/( Q(x**2) - P(x**2) )
+ */
+px = 1.0 + x * polevlf( x, P, 5 );
+
+/* multiply by power of 2 */
+x = ldexpf( px, n );
+
+return(x);
+}