Re: [Haskell-cafe] Sinus in Haskell

9 Nov 2007

      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  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