Re: amd64/144448: sin() broken in libm on amd64

[Peter Jeremy]

The following reply was made to PR standards/144448; it has been noted by GNATS.

From: Peter Jeremy <peter...@acm.org>
To: "Eugene M. Zheganin" <e...@norma.perm.ru>
Cc: freebsd-gn...@FreeBSD.org
Subject: Re: amd64/144448: sin() broken in libm on amd64
Date: Thu, 4 Mar 2010 09:39:04 +1100

--z4+8/lEcDcG5Ke9S
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline
Content-Transfer-Encoding: quoted-printable

On 2010-Mar-03 13:33:53 +0000, "Eugene M. Zheganin" <e...@norma.perm.ru> wro=
te:
>The following code, being compiled, is printing different results on
>i386 and amd64 platforms. The result on amd64 is invalid.

I agree they are different but why do you think the i386 result is valid?

Between imperfect rounding of trig functions and throwing away high
bits (when intermediate results are outside +/-pi), after iterating
20 times, your result is not much better than noise in either case.

Evaluating the following bc program at both scale=3D100 and scale=3D200
gives the same output when truncated to ~73 digits:
=3D=3D=3D=3D 8-< =3D=3D=3D=3D 8-< =3D=3D=3D=3D 8-< =3D=3D=3D=3D 8-< =3D=3D=
=3D=3D
prevval =3D 734
curval =3D s(734)
print "p",prevval, "\nc", curval, "\n"
for (i =3D 0; i < 19; i++) {
curval =3D 16 * curval
prevval =3D curval
curval =3D s(curval)
print "p",prevval, "\nc", curval, "\n"
}
=3D=3D=3D=3D 8-< =3D=3D=3D=3D 8-< =3D=3D=3D=3D 8-< =3D=3D=3D=3D 8-< =3D=3D=
=3D=3D

To demonstrate, save the above as sin.bc and run
$ echo scale=3D100 | cat - sin.bc | BC_LINE_LENGTH=3D300 bc -l | cut -c 1-7=
5 > sin.100
$ echo scale=3D200 | cat - sin.bc | BC_LINE_LENGTH=3D300 bc -l | cut -c 1-7=
5 > sin.200
$ cmp sin.100 sin.200
and the resultant files are identical - suggesting that the result is
evaluated to sufficient accuracy. I agree that this is not perfect
but don't immediately have access to any other arbitrary-precision
trig functions for cross-checking.

The last few lines of the bc calculation are:
p-3.53505577417867235125724743360716594330021904691334620895981062175595458
c.3833892004923101643681658500607744026623664419888653675724107455906899765
p6.134227207876962629890653600972390442597863071821845881158571929451039624
c-.148407850272539886063602844608156280652304697194479250111748637369605222

This is nothing like the double-precision results from either i386
or amd64.

--=20
Peter Jeremy

--z4+8/lEcDcG5Ke9S
Content-Type: application/pgp-signature
Content-Disposition: inline

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.14 (FreeBSD)

iEYEARECAAYFAkuO5QgACgkQ/opHv/APuIfkmgCdFOsBDRp5zNmXq29EDzWM/9uc
XjgAoLqk9gLQH5f+UpNO2AyzZrOZHjYi
=q5EP
-----END PGP SIGNATURE-----

--z4+8/lEcDcG5Ke9S--