dead code in lgamma_r[f]?

Bruce Evans brde at
Fri Dec 6 01:55:28 UTC 2013

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.


More information about the freebsd-numerics mailing list