From owner-freebsd-numerics@FreeBSD.ORG Tue Dec 10 05:45:18 2013 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id D0EE835B for ; Tue, 10 Dec 2013 05:45:18 +0000 (UTC) Received: from mail107.syd.optusnet.com.au (mail107.syd.optusnet.com.au [211.29.132.53]) by mx1.freebsd.org (Postfix) with ESMTP id 56B47150C for ; Tue, 10 Dec 2013 05:45:17 +0000 (UTC) Received: from c122-106-156-23.carlnfd1.nsw.optusnet.com.au (c122-106-156-23.carlnfd1.nsw.optusnet.com.au [122.106.156.23]) by mail107.syd.optusnet.com.au (Postfix) with ESMTPS id 7F258D447C0; Tue, 10 Dec 2013 16:45:15 +1100 (EST) Date: Tue, 10 Dec 2013 16:45:14 +1100 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl Subject: Re: dead code in lgamma_r[f]? In-Reply-To: <20131208183339.GA83010@troutmask.apl.washington.edu> Message-ID: <20131210162729.A1022@besplex.bde.org> 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> <20131208012939.GA80145@troutmask.apl.washington.edu> <20131208141627.F1181@besplex.bde.org> <20131208183339.GA83010@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.1 cv=Kejr72oD c=1 sm=1 tr=0 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 a=PO7r1zJSAAAA:8 a=m8NKqXhiB1EA:10 a=kj9zAlcOel0A:10 a=JzwRw_2MAAAA:8 a=ElNGO9NIuDsA:10 a=UvcvclIWdY0omR4vDwoA:9 a=mfOIlMnV6EL5W_DY:21 a=H_VbkGgs1e0FD1jH:21 a=CjuIK1q_8ugA:10 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: Tue, 10 Dec 2013 05:45:19 -0000 On Sun, 8 Dec 2013, Steve Kargl wrote: > On Sun, Dec 08, 2013 at 03:21:04PM +1100, Bruce Evans wrote: >> ... >> gamma() was the first part of libm that I looked closely at (even >> before trig functions in 2004/5). I didn't understand it at first and >> tried to write my own. My first idea was to use the asymptotic >> expansion as a kernel for all cases, especially for complex args. >> ... > Interesting idea. For x in [2,8), fdlibm code (as you probably > know) uses recursion to reduce the evaluation to a log of a > product and lgamma(2+s) with s in [0,1). One can then show that I remember a bit more about why I tried to use the asymptotic expansion. The fdlibm method in [0,8] only works for real args. > lgamma(2+s) = s * (1 - E) + s**2 * P(s) where E is Euler's constant > and P(s) is a polynomial in s. But, instead of using a minimax This works for complex s, but only for s near 0. > polynomial approximation fdlibm code uses a rational approximation > lgamma(2+s) = 0.5*s + p(s) / q(s). This is the origin of the > s0, ..., s6, r1, ..., r6 coefficients. It is indeed hard to see why the code matches the comment. The 0.5*s expansion is presumably to calculate the leading term exactly. (All ration approximations should do this, since p(s) / q(s) has an error of a couple of ulps.) s * (1 - E) would need extra precision. > A quick scan of n1256.pdf shows that C99 does not require > lgamma of a complex argument, but it does reserve the > names clgamma[fl]. Same in C11 (n1570.pdf). >> % lg(z) = >> % { >> % local(bol, rnew); >> % >> % if (abs(z - 1) < 1.e-1, return (lngamma(z));); >> % if (abs(z - 2) < 1.e-1, return (lngamma(z));); >> >> Near 1 and 2, I planned to use a poly approximation, but that is not >> so easy to do in pari so I call the pari builtin lgamma() to fake it. > > fdlibm's lgammaf_r does not appear to have an issue near 1 or 2. > Starting at the lower bound of each interval and using nextafterf > to scan up to the upper bound, I'm seeing I mean that these cases are easy (even for complex z) so I didn't bother testing various implementations of them for accuracy. >> % lgb(z) = (z - 0.5) * log(z) - z + 0.5 * log(2 * Pi) + \ >> % sum(n = 1, N, b[n] / z^(2 * n - 1)) >> >> This is the asymptotic expansion. Pari syntax is only perfect for >> numeric things like this (and probably for some number theory things >> which I don't know well). b[n] is a table of coefficients. > ... > portion of your above code is approximated by a minimax polynomial. > > /* 8.0 <= x < 2**58 */ > } else if (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 This is essentially the same -- the approximation is mostly a minimax poly in 1/z. Another problem for complex args is that it would be surpising if minimax polys worked for them. The polys use delicate cancelations that probably only work for real args. Bruce