Category : C Source Code
Archive   : CEPHES.ZIP
Filename : LOG.C

/* log.c
*
* Natural logarithm
*
*
*
* SYNOPSIS:
*
* double x, y, log();
*
* y = log( x );
*
*
*
* DESCRIPTION:
*
* Returns logarithm to the base e (2.718...) 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 range # trials peak rms
* IEEE 0.5, 2.0 35000 1.8e-16 5.5e-17
* IEEE 1, MAXNUM 10000 1.4e-16 4.8e-17
* DEC 0.5, 2.0 20000 2.0e-17 7.0e-18
* DEC 1, MAXNUM 17700 1.9e-17 6.1e-18
*
* In the tests over the interval [1, MAXNUM], the logarithms
* of the random arguments were uniformly distributed over
* [0, MAXLOG].
*
* ERROR MESSAGES:
*
* log singularity: x = 0; returns MINLOG
* log domain: x < 0; returns MINLOG
*/

/*
Cephes Math Library Release 2.0: April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

#include "mconf.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[] = {
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 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 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 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 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 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 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

/* 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[] = {
-7.90003033195003570794E-1,
1.64470394105188802674E1,
-6.44977287217605054794E1,
};
static double S[] = {
/* 1.00000000000000000000E0, */
-3.57676812896550927595E1,
3.13460384625398562220E2,
-7.73972744661126066616E2,
};
#endif

#ifdef DEC
static short R[] = {
0140112,0036643,0103520,0033473,
0041203,0111611,0063001,0116432,
0141600,0177326,0046214,0075606,
};
static short S[] = {

0141417,0011033,0005503,0043546,
0042234,0135355,0161046,0152602,
0142501,0077101,0071322,0134511,
};
#endif

#ifdef IBMPC
static short R[] = {
0x06e7,0x70ea,0x47b4,0xbfe9,
0x33a3,0x2cc0,0x7271,0x4030,
0x8f71,0xc991,0x1fda,0xc050
};
static short S[] = {

0x68ed,0x6168,0xe243,0xc041,
0xdab0,0xbc44,0x975d,0x4073,
0x5729,0x2e5a,0x2fc8,0xc088
};
#endif

#ifdef MIEEE
static short R[] = {
0xbfe9,0x47b4,0x70ea,0x06e7,
0x4030,0x7271,0x2cc0,0x33a3,
0xc050,0x1fda,0xc991,0x8f71
};
static short S[] = {
0xc041,0xe243,0x6168,0x68ed,
0x4073,0x975d,0xbc44,0xdab0,
0xc088,0x2fc8,0x2e5a,0x5729
};
#endif


#define SQRTH 0.70710678118654752440
extern double LOGE2, SQRT2, MAXLOG, MINLOG;

double log(x)
double x;
{
short e, *q;
double y, z;
double frexp(), ldexp(), polevl(), p1evl();

/* Test for domain */
if( x <= 0.0 )
{
if( x == 0.0 )
mtherr( fname, SING );
else
mtherr( fname, DOMAIN );
return( MINLOG );
}

/* 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
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 + x * ( z * polevl( z, R, 2 ) / p1evl( z, S, 3 ) );

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;
y = x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) );
y = y - ldexp( z, -1 ); /* y - 0.5 * z */
z = x + y;

ldone:

/* recombine with exponent term */
if( e != 0 )
{
y = e;
z = z - y * 2.121944400546905827679e-4;
z = z + y * 0.693359375;
}

return( z );
}