/*
 *				l o g . c
 */

/*)LIBRARY
*/

#ifdef	DOCUMENTATION
title	log	logarithm function
index	logarithm function
usage
	.s
	double x, f, log();
	.br
	f = log(x);
	.s
description
	.s
	Returns logarithm of x.
	.s
diagnostics
	.s
	If the argument is negative or zero the message 'log arg negative or
	zero', followed by the value of the argument, is written to stderr.
	If the
	argument was zero a value of -HUGE is returned, otherwise the value
	of log(fabs(x)) is returned.
	.s
internal
	.s
	Algorithm from pp. 35-59 of Cody and Waite.
	.s
author
	.s
	Hamish Ross
	.s
date
	.s
	26-Dec-84
#endif

#include <math.h>

#define C1 0.693359375		/* C1 + C2 should represent log 2 to	*/
#define C2 -2.12194440054690583e-4	/* more than machine precision	*/

static double a0 = -0.641249434237455811e+2;
static double a1 =  0.163839435630215342e+2;
static double a2 = -0.789561128874912573e+0;
static double b0 = -0.769499321084948798e+3;
static double b1 =  0.312032220919245328e+3;
static double b2 = -0.356679777390346462e+2;

double log(x)
double x;
{
    double a, b, f, frexp(), r, w, z;
    int n;

    if (x <= 0.0) {
	cmemsg(FP_NLOG, &x);
	if (x == 0.0)
	    return(-HUGE);
	else
	    x = -x;
    }
    f = frexp(x, &n);
    if (f > ROOT_05) {
	a = (f - 0.5) - 0.5;
	b = f * 0.5 + 0.5;
    }
    else {
	n--;
	a = f - 0.5;
	b = a * 0.5 + 0.5;
    }
    z = a / b;
    w = z * z;
    a = (a2 * w + a1) * w + a0;
    b = ((w + b2) * w + b1) * w + b0;
    r = z + z * w * a / b;
    z = n;
    return((z * C2 + r) + z * C1);
}