Recently, I got a mail from a user of my Euler Math Toolbox, complaining about the wrong result of sin(10^22) in Euler. The user argued that 10^22 is completely representable in the IEEE floating point system. So the result should be the correct result rounded to the next IEEE floating point number. Instead, Euler delivers a completely different result.

Euler, like almost all other numerical programs, use the floating point co-processor via the C library. This is fast, but not always transparent. However, we should develop some feeling for that arithmetic.

If we try to do the reduction of the argument to [0,2pi] within the IEEE floating point arithmetic, we will fail, since 10^22/2pi is an integer there (computed correctly). So the result of the sine function will be 0, which is just as wrong as the result by Euler or by the floating point arithmetic of Maxima. I cannot explain this result. However, it could be argued that Euler should issue an error message for such ranges.

In general, the IEEE results of functions are not guaranteed to be the correct rounded results of the precise functions. This is only true as far as it is possible without too much effort. For the normal ranges, the results should only be two significant bits (~2^(-54)) away from the precise results.

By the way, Euler’s interval arithmetic delivers sin(~10^22~) as ~-1,1~, which is the best approximation one can get. Note that 10^22 might be exact in the computer, but its error is about 10^6. The sine function does not make sense at all there.

On the other hand, fpprec:=100 in Maxima delivers a „correct“ result using the command sin(1b22).