From owner-freebsd-standards@FreeBSD.ORG Sat Jun 25 17:14:57 2005 Return-Path: X-Original-To: freebsd-standards@freebsd.org Delivered-To: freebsd-standards@freebsd.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id D4E0716A41C for ; Sat, 25 Jun 2005 17:14:57 +0000 (GMT) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.208.78.105]) by mx1.FreeBSD.org (Postfix) with ESMTP id B556743D1D for ; Sat, 25 Jun 2005 17:14:57 +0000 (GMT) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.13.4/8.13.4) with ESMTP id j5PHEvYq087058 for ; Sat, 25 Jun 2005 10:14:57 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.13.4/8.13.1/Submit) id j5PHEv24087057 for freebsd-standards@freebsd.org; Sat, 25 Jun 2005 10:14:57 -0700 (PDT) (envelope-from sgk) Date: Sat, 25 Jun 2005 10:14:57 -0700 From: Steve Kargl To: freebsd-standards@freebsd.org Message-ID: <20050625171457.GA86993@troutmask.apl.washington.edu> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.4.2.1i Subject: Implementation of logl(3) X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sat, 25 Jun 2005 17:14:58 -0000 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 #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; }