1 files changed, 128 insertions, 0 deletions
diff --git a/libm/float/logf.c b/libm/float/logf.c
new file mode 100644
index 000000000..750138564
--- /dev/null
+++ b/libm/float/logf.c
@@ -0,0 +1,128 @@
+/*							logf.c
+ *
+ *	Natural logarithm
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, logf();
+ *
+ * y = logf( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base e (2.718...) logarithm of x.
+ *
+ * The argument is separated into its exponent and fractional
+ * parts.  If the exponent is between -1 and +1, the logarithm
+ * of the fraction is approximated by
+ *
+ *     log(1+x) = x - 0.5 x**2 + x**3 P(x)
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0.5, 2.0    100000       7.6e-8     2.7e-8
+ *    IEEE      1, MAXNUMF  100000                  2.6e-8
+ *
+ * In the tests over the interval [1, MAXNUM], the logarithms
+ * of the random arguments were uniformly distributed over
+ * [0, MAXLOGF].
+ *
+ * ERROR MESSAGES:
+ *
+ * logf singularity:  x = 0; returns MINLOG
+ * logf domain:       x < 0; returns MINLOG
+ */
+
+/*
+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
+*/
+
+/* Single precision natural logarithm
+ * test interval: [sqrt(2)/2, sqrt(2)]
+ * trials: 10000
+ * peak relative error: 7.1e-8
+ * rms relative error: 2.7e-8
+ */
+
+#include <math.h>
+extern float MINLOGF, SQRTHF;
+
+
+float frexpf( float, int * );
+
+float logf( float xx )
+{
+register float y;
+float x, z, fe;
+int e;
+
+x = xx;
+fe = 0.0;
+/* Test for domain */
+if( x <= 0.0 )
+	{
+	if( x == 0.0 )
+		mtherr( "logf", SING );
+	else
+		mtherr( "logf", DOMAIN );
+	return( MINLOGF );
+	}
+
+x = frexpf( x, &e );
+if( x < SQRTHF )
+	{
+	e -= 1;
+	x = x + x - 1.0; /*  2x - 1  */
+	}	
+else
+	{
+	x = x - 1.0;
+	}
+z = x * x;
+/* 3.4e-9 */
+/*
+p = logfcof;
+y = *p++ * x;
+for( i=0; i<8; i++ )
+	{
+	y += *p++;
+	y *= x;
+	}
+y *= z;
+*/
+
+y =
+(((((((( 7.0376836292E-2 * x
+- 1.1514610310E-1) * x
++ 1.1676998740E-1) * x
+- 1.2420140846E-1) * x
++ 1.4249322787E-1) * x
+- 1.6668057665E-1) * x
++ 2.0000714765E-1) * x
+- 2.4999993993E-1) * x
++ 3.3333331174E-1) * x * z;
+
+if( e )
+	{
+	fe = e;
+	y += -2.12194440e-4 * fe;
+	}
+
+y +=  -0.5 * z;  /* y - 0.5 x^2 */
+z = x + y;   /* ... + x  */
+
+if( e )
+	z += 0.693359375 * fe;
+
+return( z );
+}