Skip site navigation (1)Skip section navigation (2)
Date:      Fri, 6 Dec 2013 13:07:57 +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:  <20131206130353.U1549@besplex.bde.org>
In-Reply-To: <20131206114426.I1329@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> <20131206114426.I1329@besplex.bde.org>

Next in thread | Previous in thread | Raw E-Mail | Index | Archive | Help
On Fri, 6 Dec 2013, 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".
>
> 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.

Oops^2.  s/easy/hard/.  When irint(z) overflows on i386, it gives an
overflow exception and a result of INT_MIN instead of benign truncation.
To avoid these problems, it would have to extract the low bits like
lgamma() does now, but that would slow down the usual case of a
severely limited range.

Bruce



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