From owner-freebsd-numerics@FreeBSD.ORG Sun Dec 8 01:29:41 2013 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id 4B08DE6D for ; Sun, 8 Dec 2013 01:29:41 +0000 (UTC) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (No client certificate requested) by mx1.freebsd.org (Postfix) with ESMTPS id 2458B16FF for ; Sun, 8 Dec 2013 01:29:41 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.7/8.14.7) with ESMTP id rB81TeuQ080225; Sat, 7 Dec 2013 17:29:40 -0800 (PST) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.7/8.14.7/Submit) id rB81Te0h080224; Sat, 7 Dec 2013 17:29:40 -0800 (PST) (envelope-from sgk) Date: Sat, 7 Dec 2013 17:29:39 -0800 From: Steve Kargl To: Bruce Evans Subject: Re: dead code in lgamma_r[f]? Message-ID: <20131208012939.GA80145@troutmask.apl.washington.edu> References: <20131205173700.GA64575@troutmask.apl.washington.edu> <20131205182324.GA65135@troutmask.apl.washington.edu> <20131205195225.GA65732@troutmask.apl.washington.edu> <20131206102724.W1033@besplex.bde.org> <20131206114426.I1329@besplex.bde.org> <20131207064602.GA76042@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20131207064602.GA76042@troutmask.apl.washington.edu> User-Agent: Mutt/1.5.21 (2010-09-15) Cc: freebsd-numerics@FreeBSD.org, Filipe Maia X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.17 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 08 Dec 2013 01:29:41 -0000 On Fri, Dec 06, 2013 at 10:46:02PM -0800, Steve Kargl wrote: > On Fri, Dec 06, 2013 at 12:55:16PM +1100, Bruce Evans wrote: > > On Fri, 6 Dec 2013, Bruce Evans wrote: > > > > > On Thu, 5 Dec 2013, Steve Kargl wrote: > > > ... > > >> If we again look at the code from __ieee754_lgamma_r(), we see > > >> that sin_pi() is called if ix < 0x43300000, so by the time we > > >> arrive at the 'if(ix<0x43300000)' statement we already know that > > >> the condition is true. > > > > > > No, only for negative integers. hx<0 classifies negative values, and > > > then ix>=0x43300000 classifies numbers that are so large negative that > > > they must be integers, and the result is a sort of pole error. We > > > are just filtering out this case, perhaps as an optimization. > > > > Oops, sin_pi() is only called for negative integers, so your change > > seems to be correct. Just add a comment about the limited domain > > of sin_pi() (it already has one saying that "x is assumed negative". > > > > I wish to retract my earlier statement that after 2 additional > years of reading fdlibm code that it was easier to work with. > I spent the better part of Friday giving myself a headache trying > to understand the algorithm for lgamma_r(). The code for x in > the interval (0,2) does not match any comment in lgamma_r(). > > I also think, but can't prove yet, that like erff() the > polynomial and rational approximations in lgammaf_r() have > too many terms. > I really need to stop looking at fdlibm code. The threshold(s) in e_lgammaf.c could be chosen better (and actually documented). In the float code /* 8.0 <= x < 2**58 */ } else if (ix < 0x4e800000) { // (ix < 0x5c800000) { t = __ieee754_logf(x); z = one/x; y = z*z; w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); r = (x-half)*(t-one)+w; } else /* 2**58 <= x <= inf */ r = x*(__ieee754_logf(x)-one); the 2**58 is much too large. Asymptotically, we have lgamma(x)~(x-0.5)log(x)-x = x*(log(x)-1) - 0.5*log(x) where the 0.5*log(x) is insignificant when compared to 2**58. I suspect that 2**58 is a sloppy threshold for double, which can be reduced to 2**30 (or smaller) for float. I suppose I could look for a tighter bound, but 2**(p+slop) with p, the precision, probably is sufficient. -- Steve