From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 22 Nov 2001 14:04:29 +0000 Subject: Totally rework the math library, this time based on the MacOs X math library (which is itself based on the math lib from FreeBSD). -Erik --- libm/double/log.c | 341 ------------------------------------------------------ 1 file changed, 341 deletions(-) delete mode 100644 libm/double/log.c (limited to 'libm/double/log.c') diff --git a/libm/double/log.c b/libm/double/log.c deleted file mode 100644 index 2fdea17a7..000000000 --- a/libm/double/log.c +++ /dev/null @@ -1,341 +0,0 @@ -/* 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