summaryrefslogtreecommitdiff
path: root/libm/float/log10f.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/float/log10f.c')
-rw-r--r--libm/float/log10f.c129
1 files changed, 129 insertions, 0 deletions
diff --git a/libm/float/log10f.c b/libm/float/log10f.c
new file mode 100644
index 000000000..6cb2e4d87
--- /dev/null
+++ b/libm/float/log10f.c
@@ -0,0 +1,129 @@
+/* log10f.c
+ *
+ * Common logarithm
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, log10f();
+ *
+ * y = log10f( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns logarithm to the base 10 of x.
+ *
+ * The argument is separated into its exponent and fractional
+ * parts. 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 1.3e-7 3.4e-8
+ * IEEE 0, MAXNUMF 100000 1.3e-7 2.6e-8
+ *
+ * In the tests over the interval [0, MAXNUM], the logarithms
+ * of the random arguments were uniformly distributed over
+ * [-MAXL10, MAXL10].
+ *
+ * ERROR MESSAGES:
+ *
+ * log10f singularity: x = 0; returns -MAXL10
+ * log10f domain: x < 0; returns -MAXL10
+ * MAXL10 = 38.230809449325611792
+ */
+
+/*
+Cephes Math Library Release 2.1: December, 1988
+Copyright 1984, 1987, 1988 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+#include <math.h>
+static char fname[] = {"log10"};
+
+/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
+ * 1/sqrt(2) <= x < sqrt(2)
+ */
+static float P[] = {
+ 7.0376836292E-2,
+-1.1514610310E-1,
+ 1.1676998740E-1,
+-1.2420140846E-1,
+ 1.4249322787E-1,
+-1.6668057665E-1,
+ 2.0000714765E-1,
+-2.4999993993E-1,
+ 3.3333331174E-1
+};
+
+
+#define SQRTH 0.70710678118654752440
+#define L102A 3.0078125E-1
+#define L102B 2.48745663981195213739E-4
+#define L10EA 4.3359375E-1
+#define L10EB 7.00731903251827651129E-4
+
+static float MAXL10 = 38.230809449325611792;
+
+float frexpf(float, int *), polevlf(float, float *, int);
+
+float log10f(float xx)
+{
+float x, y, z;
+int e;
+
+x = xx;
+/* Test for domain */
+if( x <= 0.0 )
+ {
+ if( x == 0.0 )
+ mtherr( fname, SING );
+ else
+ mtherr( fname, DOMAIN );
+ return( -MAXL10 );
+ }
+
+/* separate mantissa from exponent */
+
+x = frexpf( x, &e );
+
+/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x) */
+
+if( x < SQRTH )
+ {
+ e -= 1;
+ x = 2.0*x - 1.0;
+ }
+else
+ {
+ x = x - 1.0;
+ }
+
+
+/* rational form */
+z = x*x;
+y = x * ( z * polevlf( x, P, 8 ) );
+y = y - 0.5 * z; /* y - 0.5 * x**2 */
+
+/* multiply log of fraction by log10(e)
+ * and base 2 exponent by log10(2)
+ */
+z = (x + y) * L10EB; /* accumulate terms in order of size */
+z += y * L10EA;
+z += x * L10EA;
+x = e;
+z += x * L102B;
+z += x * L102A;
+
+
+return( z );
+}