1 files changed, 341 insertions, 0 deletions
diff --git a/libm/double/log.c b/libm/double/log.c
new file mode 100644
index 000000000..2fdea17a7
--- /dev/null
+++ b/libm/double/log.c
@@ -0,0 +1,341 @@
+/*							log.c
+ *
+ *	Natural logarithm
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, log();
+ *
+ * y = log( 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)/Q(x).
+ *
+ * Otherwise, setting  z = 2(x-1)/x+1),
+ * 
+ *     log(x) = z + z**3 P(z)/Q(z).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0.5, 2.0    150000      1.44e-16    5.06e-17
+ *    IEEE      +-MAXNUM    30000       1.20e-16    4.78e-17
+ *    DEC       0, 10       170000      1.8e-17     6.3e-18
+ *
+ * In the tests over the interval [+-MAXNUM], the logarithms
+ * of the random arguments were uniformly distributed over
+ * [0, MAXLOG].
+ *
+ * ERROR MESSAGES:
+ *
+ * log singularity:  x = 0; returns -INFINITY
+ * log 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[] = {"log"};
+
+/* 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[] = {
+ 1.01875663804580931796E-4,
+ 4.97494994976747001425E-1,
+ 4.70579119878881725854E0,
+ 1.44989225341610930846E1,
+ 1.79368678507819816313E1,
+ 7.70838733755885391666E0,
+};
+static double Q[] = {
+/* 1.00000000000000000000E0, */
+ 1.12873587189167450590E1,
+ 4.52279145837532221105E1,
+ 8.29875266912776603211E1,
+ 7.11544750618563894466E1,
+ 2.31251620126765340583E1,
+};
+#endif
+
+#ifdef DEC
+static unsigned short P[] = {
+0037777,0127270,0162547,0057274,
+0041001,0054665,0164317,0005341,
+0041451,0034104,0031640,0105773,
+0041677,0011276,0123617,0160135,
+0041701,0126603,0053215,0117250,
+0041420,0115777,0135206,0030232,
+};
+static unsigned short Q[] = {
+/*0040200,0000000,0000000,0000000,*/
+0041220,0144332,0045272,0174241,
+0041742,0164566,0035720,0130431,
+0042246,0126327,0166065,0116357,
+0042372,0033420,0157525,0124560,
+0042271,0167002,0066537,0172303,
+0041730,0164777,0113711,0044407,
+};
+#endif
+
+#ifdef IBMPC
+static unsigned short P[] = {
+0x1bb0,0x93c3,0xb4c2,0x3f1a,
+0x52f2,0x3f56,0xd6f5,0x3fdf,
+0x6911,0xed92,0xd2ba,0x4012,
+0xeb2e,0xc63e,0xff72,0x402c,
+0xc84d,0x924b,0xefd6,0x4031,
+0xdcf8,0x7d7e,0xd563,0x401e,
+};
+static unsigned short Q[] = {
+/*0x0000,0x0000,0x0000,0x3ff0,*/
+0xef8e,0xae97,0x9320,0x4026,
+0xc033,0x4e19,0x9d2c,0x4046,
+0xbdbd,0xa326,0xbf33,0x4054,
+0xae21,0xeb5e,0xc9e2,0x4051,
+0x25b2,0x9e1f,0x200a,0x4037,
+};
+#endif
+
+#ifdef MIEEE
+static unsigned short P[] = {
+0x3f1a,0xb4c2,0x93c3,0x1bb0,
+0x3fdf,0xd6f5,0x3f56,0x52f2,
+0x4012,0xd2ba,0xed92,0x6911,
+0x402c,0xff72,0xc63e,0xeb2e,
+0x4031,0xefd6,0x924b,0xc84d,
+0x401e,0xd563,0x7d7e,0xdcf8,
+};
+static unsigned short Q[] = {
+/*0x3ff0,0x0000,0x0000,0x0000,*/
+0x4026,0x9320,0xae97,0xef8e,
+0x4046,0x9d2c,0x4e19,0xc033,
+0x4054,0xbf33,0xa326,0xbdbd,
+0x4051,0xc9e2,0xeb5e,0xae21,
+0x4037,0x200a,0x9e1f,0x25b2,
+};
+#endif
+
+/* Coefficients for log(x) = z + z**3 P(z)/Q(z),
+ * where z = 2(x-1)/(x+1)
+ * 1/sqrt(2) <= x < sqrt(2)
+ */
+
+#ifdef UNK
+static double R[3] = {
+-7.89580278884799154124E-1,
+ 1.63866645699558079767E1,
+-6.41409952958715622951E1,
+};
+static double S[3] = {
+/* 1.00000000000000000000E0,*/
+-3.56722798256324312549E1,
+ 3.12093766372244180303E2,
+-7.69691943550460008604E2,
+};
+#endif
+#ifdef DEC
+static unsigned short R[12] = {
+0140112,0020756,0161540,0072035,
+0041203,0013743,0114023,0155527,
+0141600,0044060,0104421,0050400,
+};
+static unsigned short S[12] = {
+/*0040200,0000000,0000000,0000000,*/
+0141416,0130152,0017543,0064122,
+0042234,0006000,0104527,0020155,
+0142500,0066110,0146631,0174731,
+};
+#endif
+#ifdef IBMPC
+static unsigned short R[12] = {
+0x0e84,0xdc6c,0x443d,0xbfe9,
+0x7b6b,0x7302,0x62fc,0x4030,
+0x2a20,0x1122,0x0906,0xc050,
+};
+static unsigned short S[12] = {
+/*0x0000,0x0000,0x0000,0x3ff0,*/
+0x6d0a,0x43ec,0xd60d,0xc041,
+0xe40e,0x112a,0x8180,0x4073,
+0x3f3b,0x19b3,0x0d89,0xc088,
+};
+#endif
+#ifdef MIEEE
+static unsigned short R[12] = {
+0xbfe9,0x443d,0xdc6c,0x0e84,
+0x4030,0x62fc,0x7302,0x7b6b,
+0xc050,0x0906,0x1122,0x2a20,
+};
+static unsigned short S[12] = {
+/*0x3ff0,0x0000,0x0000,0x0000,*/
+0xc041,0xd60d,0x43ec,0x6d0a,
+0x4073,0x8180,0x112a,0xe40e,
+0xc088,0x0d89,0x19b3,0x3f3b,
+};
+#endif
+
+#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
+#define SQRTH 0.70710678118654752440
+extern double INFINITY, NAN;
+
+double log(x)
+double x;
+{
+int e;
+#ifdef DEC
+short *q;
+#endif
+double y, z;
+
+#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
+
+/* Note, frexp is used so that denormal numbers
+ * will be handled properly.
+ */
+#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(x) = z + z**3 P(z)/Q(z),
+ * where z = 2(x-1)/x+1)
+ */
+
+if( (e > 2) || (e < -2) )
+{
+if( x < SQRTH )
+	{ /* 2( 2x-1 )/( 2x+1 ) */
+	e -= 1;
+	z = x - 0.5;
+	y = 0.5 * z + 0.5;
+	}	
+else
+	{ /*  2 (x-1)/(x+1)   */
+	z = x - 0.5;
+	z -= 0.5;
+	y = 0.5 * x  + 0.5;
+	}
+
+x = z / y;
+
+
+/* rational form */
+z = x*x;
+z = x * ( z * polevl( z, R, 2 ) / p1evl( z, S, 3 ) );
+y = e;
+z = z - y * 2.121944400546905827679e-4;
+z = z + x;
+z = z + e * 0.693359375;
+goto ldone;
+}
+
+
+
+/* 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;
+#if DEC
+y = x * ( z * polevl( x, P, 5 ) / p1evl( x, Q, 6 ) );
+#else
+y = x * ( z * polevl( x, P, 5 ) / p1evl( x, Q, 5 ) );
+#endif
+if( e )
+	y = y - e * 2.121944400546905827679e-4;
+y = y - ldexp( z, -1 );   /*  y - 0.5 * z  */
+z = x + y;
+if( e )
+	z = z + e * 0.693359375;
+
+ldone:
+
+return( z );
+}