[Haskell-cafe] Sinus in Haskell

[Carl Witty]

On Fri, 2007-11-09 at 21:34 +0100, Daniel Fischer wrote:
> Am Freitag, 9. November 2007 21:02 schrieb Hans van Thiel:
> > On Fri, 2007-11-09 at 14:30 -0500, Brent Yorgey wrote:
> > > On Nov 9, 2007 2:08 PM, Hans van Thiel <hthiel.char at zonnet.nl> wrote:
> > >         Hello All,
> > >         Can anybody explain the results for 1.0, 2.0 and 3.0 times pi
> > >         below?
> > >         GHCi yields the same results. I did search the Haskell report
> > >         and my
> > >         text books, but to no avail. Thanks in advance,
> > >         Hans van Thiel
> > >
> > >         Hugs> sin (0.0 * pi)
> > >         0.0
> > >         Hugs> sin (0.5 * pi)
> > >         1.0
> > >         Hugs> sin (1.0 * pi)
> > >         1.22460635382238e-16
> > >         Hugs> sin (1.5 * pi)
> > >         -1.0
> > >         Hugs> sin (2.0 * pi)
> > >         -2.44921270764475e-16
> > >         Hugs> sin ( 2.5 * pi)
> > >         1.0
> > >         Hugs> sin (3.0 * pi)
> > >         3.67381906146713e-16
> > >         Hugs>
> > >
> > > More generally, this is due to the fact that floating-point numbers
> > > can only have finite precision, so a little bit of rounding error is
> > > inevitable when dealing with irrational numbers like pi.   This
> > > problem is in no way specific to Haskell.
> > >
> > > -Brent
> >
> > All right, I'd have guessed that myself, if it hadn't been for the exact
> > computation results for 0, 0.5, 1.5 and 2.5 times pi. So the rounding
> > errors are only manifest for 1.0, 2.0 and 3.0 times pi. But look at the
> > difference between sin (1.0 * pi) and sin (3.0 * pi). That's not a
> > rounding error, but a factor 3 difference.. and sin (as well as cos) are
> > modulo (2 * pi), right?
> >
> > Regards,
> > Hans
> >
> The exact results for 0, 0.5*pi and 2.5*pi aren't surprising. Leaving out the 
> more than obvious case 0, we have two cases with sin x == 1. Now what's the 
> smallest positive Double epsilon such that show (1+epsilon) /= show 1.0?
> I think, on my box it's 2^^(-52), which is roughly 2.22e-16, but the result 
> calculated is closer to 1 than that, so the result is actually represented as 
> 1.0, which by a lucky coincidence is the exact value.
> Regarding sin (1.0*pi) and sin (3.0*pi), I don't know how the sine is 
> implemented,that they return nonzero values is expected, that actually 
> sin (3.0*pi) == 3.0*sin pi does surprise me, too.
> Can anybody knowing more about the implementation of sine enlighten me?

Actually, there are about 95 million floating-point values in the
vicinity of pi/2 such that the best possible floating-point
approximation of sin on those values is exactly 1.0 (this number is
2^(53/2), where 53 is the number of mantissa bits in an IEEE
double-precision float); so sin() would be expected to return exactly
1.0 for even an inaccurate version of pi/2.

I don't know how sin is implemented on your architecture (looks like
x86?); but I can guess at enough to explain these results.

The story hinges on 3 numbers: the exact real number pi, which cannot be
represented in floating-point; the best possible IEEE double-precision
approximation to pi, which I will call Prelude.pi; and a closer
precision to pi, having perhaps 68 mantissa bits instead of 53, which I
will call approx_pi.

The value of sin(pi) is of course 0.  Since Prelude.pi is not equal to
pi, sin(Prelude.pi) is not 0; instead, it is approximately
1.22464679915e-16.  Since the derivative of sin near pi is almost
exactly -1, this is approximately (pi - Prelude.pi).  But that's not the
answer you're getting; instead, you're getting exactly (approx_pi -
Prelude.pi), where approx_pi is a 68-bit approximation to pi.

Then sin(3*Prelude.pi) should be approximately (3*pi - 3*Prelude.pi);
but you're getting exactly (3*approx_pi - 3*Prelude.pi), which is
3*(approx_pi - Prelude.pi), which is 3*sin(Prelude.pi).  (Note that
3*Prelude.pi can be computed exactly -- with no further rounding -- in
IEEE double-precision floating-point.)

The above essay was written after much experimentation using the MPFR
library for correctly-rounded arbitrary-precision floating point, as
exposed in the Sage computer algebra system.

Carl Witty