bin/144306: [libm] [patch] Nasty bug in jn(3)

[Bruce Evans]

On Wed, 10 Nov 2010, Bruce Evans wrote:

> On Tue, 9 Nov 2010, Ulrich [utf-8] Spörlein wrote:
>> FreeBSD -CURRENT:
>> 
>> freebsd64# cc -o testjn -lm testjn.c && ./testjn
>> 2 4.317548e-01
>> 3 1.989999e-01
>> 4 6.474667e-02
>> 5 1.638924e-02
>> 6 3.404818e-03
>> 7 6.006884e-04
>> 8 9.216579e-05
>> 9 1.251727e-05
>
> FreeBSD agrees perfectly with pari.  This is not surprising, since only
> 7 significant decimal digits are printed and almost 15 should be correct
> -- any error here would be an error of about 100 million decimal ulps
> or a couple of binary gulps.
>
>> freebsd64# cc -o testjnf -lm testjnf.c && ./testjnf
>> 2 4.317548e-01
>> 3 1.989999e-01
>> 4 6.474666e-02
>> 5 1.638924e-02
>> 6 3.404818e-03
>> 7 6.006881e-04
>> 8 9.216577e-05
>> 9 1.251727e-05
>
> This is harder to explain.  Now we should only expect 6 digits of accuracy,
> but we are getting 7.  Anyway, these are all perfectly rounded according to
> pari.

Actually we are getting 6.  The last digit differs from the double precision
case for args 4, 7 and 8.  This is as expected.

> [pari output]
> Single precision with bad decimal field width:
> 2 4.317548e-1
> 3 1.989999e-1
> 4 6.474667e-2
> 5 1.638924e-2
> 6 3.404818e-3
> 7 6.006884e-4
> 8 9.216579e-5
> 9 1.251727e-5

Unfortunately, the correctly rounded (in single precision) result gives
(after further rounding to 1 too many decimal digits) values that agree
with the double precision C output and thus disagree with the single
precision C output.  Therefore, the C results in single precision cannot
be correctly rounded.

I may be confused or have a bug somewhere -- manually rounding the
above single precision results to 6 decimal digits gives the same
values in all cases, but FLT_DIG = 6 is supposed to give enough digits
for preserving a single precision result if it is printed and scanned
back.  When I started replying to this, I thought that extra decimal
digits were needed (why is DECIMAL_DIG = 21 3 more than LDBL_DIG = 18
on i386...?), and printed 3 extra.  I left these out to simplify
comparisons.

Bruce