From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/double/log.c | 341 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 341 insertions(+) create mode 100644 libm/double/log.c (limited to 'libm/double/log.c') 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 +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 ); +} -- cgit v1.2.3