lgamma_r and lgammaf_r return the wrong sign for -0.f

[Steve Kargl]

On Fri, Sep 12, 2014 at 02:15:37PM -0700, enh via freebsd-numerics wrote:
> if I pass -0.f to lgammaf_r, the sign returned in *signgamp is 1. this
> is incorrect --- it should be -1.
> 
> both lgamma_r and lgammaf_r are affected, but the other special cases
> in those functions look fine to me.
> 
> this is fixed in OpenBSD and glibc, but FreeBSD and NetBSD both have
> the same bug.

Are you sure FreeBSD has a bug?  From n1256.pdf (committee draft of C99) 

  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-pointr 
     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.

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.

> btw, it looks like you're missing coshl/sinhl/tanhl for ld128 now?
> (they've been removed from imprecise.c without having ld128
> implementations added afaics.)

The ld80 and ld128 versions are in the same file.

troutmask:fvwm:kargl[226] grep LDBL /usr/src/lib/msun/src/e_coshl.c 
#if LDBL_MAX_EXP != 0x4000
#define BIAS    (LDBL_MAX_EXP - 1)
#if LDBL_MANT_DIG == 64
#elif LDBL_MANT_DIG == 113
#endif /* LDBL_MANT_DIG == 64 */
        GET_LDBL_EXPSIGN(ix,x);
            if (ix<BIAS-(LDBL_MANT_DIG+1)/2)    /* |x| < TINY */
#if LDBL_MANT_DIG == 64
#elif LDBL_MANT_DIG == 113

see also e_sinhl.c and s_tanhl.c.

-- 
Steve