uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 250 insertions, 0 deletions

diff --git a/libm/double/log10.c b/libm/double/log10.c
new file mode 100644
index 000000000..7dc72e253
--- /dev/null
+++ b/libm/double/log10.c
@@ -0,0 +1,250 @@
+/*							log10.c
+ *
+ *	Common logarithm
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, log10();
+ *
+ * y = log10( 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)/Q(x).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0.5, 2.0     30000      1.5e-16     5.0e-17
+ *    IEEE      0, MAXNUM    30000      1.4e-16     4.8e-17
+ *    DEC       1, MAXNUM    50000      2.5e-17     6.0e-18
+ *
+ * In the tests over the interval [1, MAXNUM], the logarithms
+ * of the random arguments were uniformly distributed over
+ * [0, MAXLOG].
+ *
+ * ERROR MESSAGES:
+ *
+ * log10 singularity:  x = 0; returns -INFINITY
+ * log10 domain:       x < 0; returns NAN
+ */
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1995, 2000 by Stephen L. Moshier
+*/
+
+#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)
+ */
+#ifdef UNK
+static double P[] = {
+  4.58482948458143443514E-5,
+  4.98531067254050724270E-1,
+  6.56312093769992875930E0,
+  2.97877425097986925891E1,
+  6.06127134467767258030E1,
+  5.67349287391754285487E1,
+  1.98892446572874072159E1
+};
+static double Q[] = {
+/* 1.00000000000000000000E0, */
+  1.50314182634250003249E1,
+  8.27410449222435217021E1,
+  2.20664384982121929218E2,
+  3.07254189979530058263E2,
+  2.14955586696422947765E2,
+  5.96677339718622216300E1
+};
+#endif
+
+#ifdef DEC
+static unsigned short P[] = {
+0034500,0046473,0051374,0135174,
+0037777,0037566,0145712,0150321,
+0040722,0002426,0031543,0123107,
+0041356,0046513,0170752,0004346,
+0041562,0071553,0023536,0163343,
+0041542,0170221,0024316,0114216,
+0041237,0016454,0046611,0104602
+};
+static unsigned short Q[] = {
+/*0040200,0000000,0000000,0000000,*/
+0041160,0100260,0067736,0102424,
+0041645,0075552,0036563,0147072,
+0042134,0125025,0021132,0025320,
+0042231,0120211,0046030,0103271,
+0042126,0172241,0052151,0120426,
+0041556,0125702,0072116,0047103
+};
+#endif
+
+#ifdef IBMPC
+static unsigned short P[] = {
+0x974f,0x6a5f,0x09a7,0x3f08,
+0x5a1a,0xd979,0xe7ee,0x3fdf,
+0x74c9,0xc66c,0x40a2,0x401a,
+0x411d,0x7e3d,0xc9a9,0x403d,
+0xdcdc,0x64eb,0x4e6d,0x404e,
+0xd312,0x2519,0x5e12,0x404c,
+0x3130,0x89b1,0xe3a5,0x4033
+};
+static unsigned short Q[] = {
+/*0x0000,0x0000,0x0000,0x3ff0,*/
+0xd0a2,0x0dfb,0x1016,0x402e,
+0x79c7,0x47ae,0xaf6d,0x4054,
+0x455a,0xa44b,0x9542,0x406b,
+0x10d7,0x2983,0x3411,0x4073,
+0x3423,0x2a8d,0xde94,0x406a,
+0xc9c8,0x4e89,0xd578,0x404d
+};
+#endif
+
+#ifdef MIEEE
+static unsigned short P[] = {
+0x3f08,0x09a7,0x6a5f,0x974f,
+0x3fdf,0xe7ee,0xd979,0x5a1a,
+0x401a,0x40a2,0xc66c,0x74c9,
+0x403d,0xc9a9,0x7e3d,0x411d,
+0x404e,0x4e6d,0x64eb,0xdcdc,
+0x404c,0x5e12,0x2519,0xd312,
+0x4033,0xe3a5,0x89b1,0x3130
+};
+static unsigned short Q[] = {
+0x402e,0x1016,0x0dfb,0xd0a2,
+0x4054,0xaf6d,0x47ae,0x79c7,
+0x406b,0x9542,0xa44b,0x455a,
+0x4073,0x3411,0x2983,0x10d7,
+0x406a,0xde94,0x2a8d,0x3423,
+0x404d,0xd578,0x4e89,0xc9c8
+};
+#endif
+
+#define SQRTH 0.70710678118654752440
+#define L102A 3.0078125E-1
+#define L102B 2.48745663981195213739E-4
+#define L10EA 4.3359375E-1
+#define L10EB 7.00731903251827651129E-4
+
+#ifdef ANSIPROT
+extern double frexp ( double, int * );
+extern double ldexp ( double, int );
+extern double polevl ( double, void *, int );
+extern double p1evl ( double, void *, int );
+extern int isnan ( double );
+extern int isfinite ( double );
+#else
+double frexp(), ldexp(), polevl(), p1evl();
+int isnan(), isfinite();
+#endif
+extern double LOGE2, SQRT2, INFINITY, NAN;
+
+double log10(x)
+double x;
+{
+VOLATILE double z;
+double y;
+#ifdef DEC
+short *q;
+#endif
+int e;
+
+#ifdef NANS
+if( isnan(x) )
+	return(x);
+#endif
+#ifdef INFINITIES
+if( x == INFINITY )
+	return(x);
+#endif
+/* Test for domain */
+if( x <= 0.0 )
+	{
+	if( x == 0.0 )
+	        {
+		mtherr( fname, SING );
+		return( -INFINITY );
+	        }
+	else
+	        {
+		mtherr( fname, DOMAIN );
+		return( NAN );
+	        }
+	}
+
+/* separate mantissa from exponent */
+
+#ifdef DEC
+q = (short *)&x;
+e = *q;			/* short containing exponent */
+e = ((e >> 7) & 0377) - 0200;	/* the exponent */
+*q &= 0177;	/* strip exponent from x */
+*q |= 040000;	/* x now between 0.5 and 1 */
+#endif
+
+#ifdef IBMPC
+x = frexp( x, &e );
+/*
+q = (short *)&x;
+q += 3;
+e = *q;
+e = ((e >> 4) & 0x0fff) - 0x3fe;
+*q &= 0x0f;
+*q |= 0x3fe0;
+*/
+#endif
+
+/* Equivalent C language standard library function: */
+#ifdef UNK
+x = frexp( x, &e );
+#endif
+
+#ifdef MIEEE
+x = frexp( x, &e );
+#endif
+
+/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
+
+if( x < SQRTH )
+	{
+	e -= 1;
+	x = ldexp( x, 1 ) - 1.0; /*  2x - 1  */
+	}	
+else
+	{
+	x = x - 1.0;
+	}
+
+
+/* rational form */
+z = x*x;
+y = x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) );
+y = y - ldexp( z, -1 );   /*  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;
+z += e * L102B;
+z += e * L102A;
+
+
+return( z );
+}