[Haskell-cafe] Sinus in Haskell

[Daniel Fischer]

Am Samstag, 10. November 2007 00:36 schrieb Carl Witty:
> 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.

arch says i686
>
> 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.

Ah ha!
So the library sine function uses a better approximation (I've heard about 
excess-precision using 80 bit floating point numbers, but didn't think of 
that). Then it's clear.
>
> 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

Thanks a lot.

Since you seem to know a lot about these things, out of curiosity, do you know 
how these functions are actually implemented? Do they use Taylor series or 
other techniques?

Cheers,
Daniel