ynl (5, LDBL_MIN) should, like yn (5, DBL_MIN) and ynf (5, FLT_MIN), return -Inf, rather than NaN: #include <float.h> #include <math.h> #include <stdio.h> int main (void) { float f = ynf (10, FLT_MIN); printf ("%f\n", f); double d = yn (10, DBL_MIN); printf ("%f\n", d); long double ld = ynl (10, LDBL_MIN); printf ("%Lf\n", ld); return ! isfinite (ld); }
According to Joseph: "My guess is that the overflow of yn has involved an internal computation such as Inf - Inf or Inf / Inf that generated NaNs. It looks rather like the ynl code that tries to stop this happening is confused about GET_LDOUBLE_WORDS, which appears to set the sign+exponent value to a *sign-extended* value (so 0xffffffff not 0xffff for -Inf)."
This patch seems to fix it: --- sysdeps/ieee754/ldbl-96/e_jnl.c.mp2 2012-05-28 15:14:05.575192340 +0200 +++ sysdeps/ieee754/ldbl-96/e_jnl.c 2012-05-28 15:04:26.493452839 +0200 @@ -361,7 +361,8 @@ b = __ieee754_y1l (x); /* quit if b is -inf */ GET_LDOUBLE_WORDS (se, i0, i1, b); - for (i = 1; i < n && se != 0xffff; i++) + /* Use 0xffffffff since GET_LDOUBLE_WORDS sign-extends SE. */ + for (i = 1; i < n && se != 0xffffffff; i++) { temp = b; b = ((long double) (i + i) / x) * b - a;
The proposed patch may "work", but I think its wrong. If se is a signed type, then it's never equal to 0xffffffff; only reliance on confusing promotion rules makes it work, and it would not work with different type sizes. It should be compared against -1, or (better) the types should be fixed to be unsigned.
se, i0, i1 are all unsigned.
I suppose I should have checked... but that makes it even more disturbing that sign extension is happening.
Because the bitfield is int:16. Note that it is implementation defined whether an int bitfield is signed, though I suppose no GCC target has unsigned bitfields.
Fixed by 541428fecf21cdde271acbd280c53bfe5beaafe2.