Output of file : LOG.C contained in archive : CEPHES22.ZIP
/* 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 35000 1.8e-16 5.6e-17
* IEEE 1, MAXNUM 10000 1.8e-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.1: December, 1988
Copyright 1984, 1987, 1988 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[3] = {
-7.89580278884799154124E-1,
1.63866645699558079767E1,
-6.41409952958715622951E1,
};
static double S[3] = {
/* 1.00000000000000000000E0,*/
-3.56722798256324312549E1,
3.12093766372244180303E2,
-7.69691943550460008604E2,
};
#endif
#ifdef DEC
static short R[12] = {
0140112,0020756,0161540,0072035,
0041203,0013743,0114023,0155527,
0141600,0044060,0104421,0050400,
};
static short S[12] = {
/*0040200,0000000,0000000,0000000,*/
0141416,0130152,0017543,0064122,
0042234,0006000,0104527,0020155,
0142500,0066110,0146631,0174731,
};
#endif
#ifdef IBMPC
static short R[12] = {
0x0e84,0xdc6c,0x443d,0xbfe9,
0x7b6b,0x7302,0x62fc,0x4030,
0x2a20,0x1122,0x0906,0xc050,
};
static short S[12] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x6d0a,0x43ec,0xd60d,0xc041,
0xe40e,0x112a,0x8180,0x4073,
0x3f3b,0x19b3,0x0d89,0xc088,
};
#endif
#ifdef MIEEE
static short R[12] = {
0xbfe9,0x443d,0xdc6c,0x0e84,
0x4030,0x62fc,0x7302,0x7b6b,
0xc050,0x0906,0x1122,0x2a20,
};
static short S[12] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0xc041,0xd60d,0x43ec,0x6d0a,
0x4073,0x8180,0x112a,0xe40e,
0xc088,0x0d89,0x19b3,0x3f3b,
};
#endif

#define SQRTH 0.70710678118654752440
extern double MINLOG;

double log(x)
double x;
{
int e;
short *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

/* 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 + 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 );
}

3 Responses to “Category : C Source CodeArchive   : CEPHES22.ZIPFilename : LOG.C”

1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/