Implementation of logl(3)

[Steve Kargl]

I've noticed that several of C99's long double math
function are not implemented.  In an effort to start
filling in some of the holes, I've implemented logl(3).
Before I proceed in looking at other functions, I 
would appreciate comments because I'm not a number
theorist by training.

The code uses a Taylor series expansion to evaluate the
log() of the fractional portion of x = f 2**n where
f = [0.5,1).  Although log(f) is very smooth over the
indicated range, my attempts to fit a simple polynomial
to the curve have failed.  The failures are probably
due to my fitting routines which are in double not long
double (actually they use Fortran's double precision).

Note, if the following is teemed okay, then cbrtl, sqrtl,
and the other logrithm functions can be implemented in
a similar manner.

-- 
steve


/*
  ln(x) = ln(f * 2**n)
        = ln(f) + n * ln(2) where f = [0.5,1).

  Use a Taylor series about f0 = 0.75 to compute ln(f).

  ln(f) = ln(f0) + sum (1/n) * (-1)**(n-1) * ((f - f0)/f0)**n 
*/

#include <math.h>

#define LOGL2	 6.931471805599453094e-1L
#define LOGLF0	-2.876820724517809274e-1L

static long double zero = 0.0L;

#define NUM 32

long double logl(long double x) {

	int i, n, s = 1;
	long double f, t, val;

	/* log(x) = sNaN if x < 0. */
	if (x < 0.0L)
		return (x - x) / (x - x);
 
	/* log(+Inf) = +Inf */
	if (isinf(x))
		return x*x+x;

	/* log(+-0) = -Inf with signal */
	if (x == 0.0L)
		return - 1.0L / zero;

	/* log(+-0) = -Inf with signal */
	if (isnan(x))
		return (x+x);

	f = frexpl(x, &n);
	f = t = (f - 0.75L) / 0.75L;
	val = LOGLF0;
	for (i = 1; i < NUM; i++) {
		val += s * t / (long double) i;
		t *= f;
		s = -s;
	} 
	val += n * LOGL2;
	return val;
}