Date: Thu, 11 Nov 2010 00:14:18 +1100 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Bruce Evans <brde@optusnet.com.au> Cc: Ulrich =?utf-8?B?U3DDtnJsZWlu?= <uqs@spoerlein.net>, freebsd-bugs@freebsd.org Subject: Re: bin/144306: [libm] [patch] Nasty bug in jn(3) Message-ID: <20101110235945.K1737@besplex.bde.org> In-Reply-To: <20101110225059.M1461@besplex.bde.org> References: <201011092010.oA9KABNt076837@freefall.freebsd.org> <20101110225059.M1461@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
This message is in MIME format. The first part should be readable text, while the remaining parts are likely unreadable without MIME-aware tools. --0-419011604-1289394858=:1737 Content-Type: TEXT/PLAIN; charset=X-UNKNOWN; format=flowed Content-Transfer-Encoding: QUOTED-PRINTABLE On Wed, 10 Nov 2010, Bruce Evans wrote: > On Tue, 9 Nov 2010, Ulrich [utf-8] Sp=C3=B6rlein wrote: >> FreeBSD -CURRENT: >>=20 >> 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 accurac= y, > 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 precisio= n 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 =3D 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 =3D 21 3 more than LDBL_DIG =3D 18 on i386...?), and printed 3 extra. I left these out to simplify comparisons. Bruce --0-419011604-1289394858=:1737--
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20101110235945.K1737>