lgamma_r and lgammaf_r return the wrong sign for -0.f
Steve Kargl
sgk at troutmask.apl.washington.edu
Fri Sep 12 23:37:50 UTC 2014
On Fri, Sep 12, 2014 at 03:51:52PM -0700, enh wrote:
> On Fri, Sep 12, 2014 at 3:41 PM, Steve Kargl
> <sgk at troutmask.apl.washington.edu> wrote:
> > On Fri, Sep 12, 2014 at 03:17:02PM -0700, enh via freebsd-numerics wrote:
> >> On Fri, Sep 12, 2014 at 3:05 PM, Steve Kargl
> >> > F.9.5.3 The lgamma functions
> >> > -- lgamma(1) returns +0.
> >> > -- lgamma(2) returns +0.
> >> > -- lgamma(x) returns +inf and raises the ``divide-by-zero'' floating-point
> >> > exception for x a negative integer or zero.
> >> > -- lgamma(-inf) returns +inf.
> >> > -- lgamma(+inf) returns +inf.
> >> >
> >> > See the 3rd bullet. -0 is 0 and -0 is a negative integer.
> >
> > Of course, neither POSIX nor ISO C specify lgamma_r() only lgamma.
> > I just spent too much time on the 'divide-by-zero' bug, and
> > conflated that with signgam.
>
> i think the reporter's feeling was that *signgamp should be -1 because
> as x approaches 0 from the negative side, Gamma(x) approaches
> -Infinity. (though i suspect they noticed this when porting code from
> glibc.)
Mathematically, gamma(x) -> +inf as x -> 0 for x > 0 and
gamma(x) -> -inf as x -> 0 for x < 0. Assuming that a processor
supports a signed zero, it happens that x = +-0 is the only
integral floating point value that allows one to infer a sign
for signgam. All other negative integral values, of course,
do not have this possible, so signgam is arbitrarily assigned
1.
> >> > POSIX appears to defer to ISO C. n1570.pdf (committe draft for C11)
> >> > has (almost?) identical text.
> >> >
> >> >> patch below (whitespace mangled courtesy of gmail). i'd prefer to wait
> >> >> for this to be fixed in FreeBSD and pull down the fix rather than just
> >> >> fix it locally.
> >> >
> >> > I have a bigger patch coming with ld80 and ld128 version of lgammal
> >> > and lgammal_r.
> >>
> >> sorry, i haven't even looked at the *l variants.
> >
> > I only just finished writing the *l variants. It took me
> > a long time to unwind the comments in e_lgamma_r.c,
> > so that I could write the long double versions.
>
> is there a freebsd-numerics-commits or equivalent i could subscribe to?
No. I tend to post my patches here. Then Bruce Evans and I
tend to go off on private email threads with code reviews. A few
months later I get around to address issues raised by Bruce, and
submit a new patch. This cycle goes on for a year or so, and
then I finally commit my patches.
There are also a few developers (most MIPS/ARMS) that have been
fixing fenv.h issues. They commit their patches without posting
to freebsd-numerics.
I'll note that my lgamma.diff file is clocking in at 1400 LOC.
mailman may eat my post.
In the end, it is probably prudent to simply look periodically
at FreeBSD's src/lib/msun through the svnweb interface.
--
Steve
More information about the freebsd-numerics
mailing list