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 _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
