[Haskell-cafe] Sinus in Haskell
Daniel Fischer
daniel.is.fischer at web.de
Fri Nov 9 19:29:09 EST 2007
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
More information about the Haskell-Cafe
mailing list