diff options
Diffstat (limited to 'libm/double/log10.c')
| -rw-r--r-- | libm/double/log10.c | 250 |
1 files changed, 0 insertions, 250 deletions
diff --git a/libm/double/log10.c b/libm/double/log10.c deleted file mode 100644 index 7dc72e253..000000000 --- a/libm/double/log10.c +++ /dev/null @@ -1,250 +0,0 @@ -/* 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 ); -} |