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