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 <[EMAIL PROTECTED]> 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 _______________________________________________ Haskell-Cafe mailing list [email protected] http://www.haskell.org/mailman/listinfo/haskell-cafe
