uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 223 insertions, 0 deletions

diff --git a/libm/double/exp10.c b/libm/double/exp10.c
new file mode 100644
index 000000000..dd0e5a48f
--- /dev/null
+++ b/libm/double/exp10.c
@@ -0,0 +1,223 @@
+/*							exp10.c
+ *
+ *	Base 10 exponential function
+ *      (Common antilogarithm)
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, exp10();
+ *
+ * y = exp10( 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).
+ * The Pade' form
+ *
+ *    1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
+ *
+ * is used to approximate 10**f.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     -307,+307    30000       2.2e-16     5.5e-17
+ * Test result from an earlier version (2.1):
+ *    DEC       -38,+38     70000       3.1e-17     7.0e-18
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * exp10 underflow    x < -MAXL10        0.0
+ * exp10 overflow     x > MAXL10       MAXNUM
+ *
+ * DEC arithmetic: MAXL10 = 38.230809449325611792.
+ * IEEE arithmetic: MAXL10 = 308.2547155599167.
+ *
+ */
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1991, 2000 by Stephen L. Moshier
+*/
+
+
+#include <math.h>
+
+#ifdef UNK
+static double P[] = {
+ 4.09962519798587023075E-2,
+ 1.17452732554344059015E1,
+ 4.06717289936872725516E2,
+ 2.39423741207388267439E3,
+};
+static double Q[] = {
+/* 1.00000000000000000000E0,*/
+ 8.50936160849306532625E1,
+ 1.27209271178345121210E3,
+ 2.07960819286001865907E3,
+};
+/* static double LOG102 = 3.01029995663981195214e-1; */
+static double LOG210 = 3.32192809488736234787e0;
+static double LG102A = 3.01025390625000000000E-1;
+static double LG102B = 4.60503898119521373889E-6;
+/* static double MAXL10 = 38.230809449325611792; */
+static double MAXL10 = 308.2547155599167;
+#endif
+
+#ifdef DEC
+static unsigned short P[] = {
+0037047,0165657,0114061,0067234,
+0041073,0166243,0123052,0144643,
+0042313,0055720,0024032,0047443,
+0043025,0121714,0070232,0050007,
+};
+static unsigned short Q[] = {
+/*0040200,0000000,0000000,0000000,*/
+0041652,0027756,0071216,0050075,
+0042637,0001367,0077263,0136017,
+0043001,0174673,0024157,0133416,
+};
+/*
+static unsigned short L102[] = {0037632,0020232,0102373,0147770};
+#define LOG102 *(double *)L102
+*/
+static unsigned short L210[] = {0040524,0115170,0045715,0015613};
+#define LOG210 *(double *)L210
+static unsigned short L102A[] = {0037632,0020000,0000000,0000000,};
+#define LG102A *(double *)L102A
+static unsigned short L102B[] = {0033632,0102373,0147767,0114220,};
+#define LG102B *(double *)L102B
+static unsigned short MXL[] = {0041430,0166131,0047761,0154130,};
+#define MAXL10 ( *(double *)MXL )
+#endif
+
+#ifdef IBMPC
+static unsigned short P[] = {
+0x2dd4,0xf306,0xfd75,0x3fa4,
+0x5934,0x74c5,0x7d94,0x4027,
+0x49e4,0x0503,0x6b7a,0x4079,
+0x4a01,0x8e13,0xb479,0x40a2,
+};
+static unsigned short Q[] = {
+/*0x0000,0x0000,0x0000,0x3ff0,*/
+0xca08,0xce51,0x45fd,0x4055,
+0x7782,0xefd6,0xe05e,0x4093,
+0xf6e2,0x650d,0x3f37,0x40a0,
+};
+/*
+static unsigned short L102[] = {0x79ff,0x509f,0x4413,0x3fd3};
+#define LOG102 *(double *)L102
+*/
+static unsigned short L210[] = {0xa371,0x0979,0x934f,0x400a};
+#define LOG210 *(double *)L210
+static unsigned short L102A[] = {0x0000,0x0000,0x4400,0x3fd3,};
+#define LG102A *(double *)L102A
+static unsigned short L102B[] = {0xf312,0x79fe,0x509f,0x3ed3,};
+#define LG102B *(double *)L102B
+static double MAXL10 = 308.2547155599167;
+#endif
+
+#ifdef MIEEE
+static unsigned short P[] = {
+0x3fa4,0xfd75,0xf306,0x2dd4,
+0x4027,0x7d94,0x74c5,0x5934,
+0x4079,0x6b7a,0x0503,0x49e4,
+0x40a2,0xb479,0x8e13,0x4a01,
+};
+static unsigned short Q[] = {
+/*0x3ff0,0x0000,0x0000,0x0000,*/
+0x4055,0x45fd,0xce51,0xca08,
+0x4093,0xe05e,0xefd6,0x7782,
+0x40a0,0x3f37,0x650d,0xf6e2,
+};
+/*
+static unsigned short L102[] = {0x3fd3,0x4413,0x509f,0x79ff};
+#define LOG102 *(double *)L102
+*/
+static unsigned short L210[] = {0x400a,0x934f,0x0979,0xa371};
+#define LOG210 *(double *)L210
+static unsigned short L102A[] = {0x3fd3,0x4400,0x0000,0x0000,};
+#define LG102A *(double *)L102A
+static unsigned short L102B[] = {0x3ed3,0x509f,0x79fe,0xf312,};
+#define LG102B *(double *)L102B
+static double MAXL10 = 308.2547155599167;
+#endif
+
+#ifdef ANSIPROT
+extern double floor ( double );
+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 floor(), ldexp(), polevl(), p1evl();
+int isnan(), isfinite();
+#endif
+extern double MAXNUM;
+#ifdef INFINITIES
+extern double INFINITY;
+#endif
+
+double exp10(x)
+double x;
+{
+double px, xx;
+short n;
+
+#ifdef NANS
+if( isnan(x) )
+	return(x);
+#endif
+if( x > MAXL10 )
+	{
+#ifdef INFINITIES
+	return( INFINITY );
+#else
+	mtherr( "exp10", OVERFLOW );
+	return( MAXNUM );
+#endif
+	}
+
+if( x < -MAXL10 )	/* Would like to use MINLOG but can't */
+	{
+#ifndef INFINITIES
+	mtherr( "exp10", UNDERFLOW );
+#endif
+	return(0.0);
+	}
+
+/* Express 10**x = 10**g 2**n
+ *   = 10**g 10**( n log10(2) )
+ *   = 10**( g + n log10(2) )
+ */
+px = floor( LOG210 * x + 0.5 );
+n = px;
+x -= px * LG102A;
+x -= px * LG102B;
+
+/* rational approximation for exponential
+ * of the fractional part:
+ * 10**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
+ */
+xx = x * x;
+px = x * polevl( xx, P, 3 );
+x =  px/( p1evl( xx, Q, 3 ) - px );
+x = 1.0 + ldexp( x, 1 );
+
+/* multiply by power of 2 */
+x = ldexp( x, n );
+
+return(x);
+}