Date: Sat, 25 Jun 2005 10:14:57 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-standards@freebsd.org Subject: Implementation of logl(3) Message-ID: <20050625171457.GA86993@troutmask.apl.washington.edu>
next in thread | raw e-mail | index | archive | help
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; }
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20050625171457.GA86993>