/* logf.c * * Natural logarithm * * * * SYNOPSIS: * * float x, y, logf(); * * y = logf( 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) * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0.5, 2.0 100000 7.6e-8 2.7e-8 * IEEE 1, MAXNUMF 100000 2.6e-8 * * In the tests over the interval [1, MAXNUM], the logarithms * of the random arguments were uniformly distributed over * [0, MAXLOGF]. * * ERROR MESSAGES: * * logf singularity: x = 0; returns MINLOG * logf domain: x < 0; returns MINLOG */ /* Cephes Math Library Release 2.2: June, 1992 Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */ /* Single precision natural logarithm * test interval: [sqrt(2)/2, sqrt(2)] * trials: 10000 * peak relative error: 7.1e-8 * rms relative error: 2.7e-8 */ #include <math.h> extern float MINLOGF, SQRTHF; float frexpf( float, int * ); float logf( float xx ) { register float y; float x, z, fe; int e; x = xx; fe = 0.0; /* Test for domain */ if( x <= 0.0 ) { if( x == 0.0 ) mtherr( "logf", SING ); else mtherr( "logf", DOMAIN ); return( MINLOGF ); } x = frexpf( x, &e ); if( x < SQRTHF ) { e -= 1; x = x + x - 1.0; /* 2x - 1 */ } else { x = x - 1.0; } z = x * x; /* 3.4e-9 */ /* p = logfcof; y = *p++ * x; for( i=0; i<8; i++ ) { y += *p++; y *= x; } y *= z; */ y = (((((((( 7.0376836292E-2 * x - 1.1514610310E-1) * x + 1.1676998740E-1) * x - 1.2420140846E-1) * x + 1.4249322787E-1) * x - 1.6668057665E-1) * x + 2.0000714765E-1) * x - 2.4999993993E-1) * x + 3.3333331174E-1) * x * z; if( e ) { fe = e; y += -2.12194440e-4 * fe; } y += -0.5 * z; /* y - 0.5 x^2 */ z = x + y; /* ... + x */ if( e ) z += 0.693359375 * fe; return( z ); }