Date:Fri, 6 Dec 2013 12:55:16 +1100 (EST)From:Bruce Evans <brde@optusnet.com.au>To:Bruce Evans <brde@optusnet.com.au>Cc:freebsd-numerics@FreeBSD.org, Filipe Maia <filipe.c.maia@gmail.com>, Steve Kargl <sgk@troutmask.apl.washington.edu>Subject:Re: dead code in lgamma_r[f]?Message-ID:<20131206114426.I1329@besplex.bde.org>In-Reply-To:<20131206102724.W1033@besplex.bde.org>References:<20131205173700.GA64575@troutmask.apl.washington.edu> <20131205182324.GA65135@troutmask.apl.washington.edu> <CAN5hRiV6N1arc4RDv=9JbRiXp-J9o3WzAbeZSOpAxks2ZeG+_w@mail.gmail.com> <20131205195225.GA65732@troutmask.apl.washington.edu> <20131206102724.W1033@besplex.bde.org>

Next in thread | Previous in thread | Raw E-Mail | Index | Archive | Help

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". With the limited range, we can improve more things. floor() can be replaced by adding and subtracting 0x1p52 or 0x1.8p52 like we do for trig and exp functions, followed by an adjustment to get floor() if necessary. We can use irint(z) instead of bit magic or (int)z to extract the parity bit of z z == y case. Note that fdlibm can't use (int)z since this can overflow, and overflow is often not benign in practice. irint(z) can overflow too, but it is easy to specify irint() as having benign behaviour. Bruce

Want to link to this message? Use this URL: <http://docs.FreeBSD.org/cgi/mid.cgi?20131206114426.I1329>