[Haskell-cafe] Sinus in Haskell

[Daniel Fischer]

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?

Cheers,
Daniel