(LONG and about floating point, so I suspect many Haskellers are not
going to be interested in this message . . .)
Julian Assange wrote:
>
> The precission and or rounding used by hugs/ghc seems strange, to wit:
>
> Prelude> sin(pi)
> -8.74228e-08
> Prelude> pi
> 3.14159
> sin(3.14159265358979323846)
> -8.74228e-08
I presume Hugs is using single precision floats . . .
>
> ghc:
>
> module Main where
> main = do
> print pi
> print (sin pi)
>
> ./a.out
> 3.141592653589793
> 1.2246467991473532e-16
and GHC is using double precisions. Let's summarise what IEEE 754 gives you.
IEEE 754 only defines the basic arithmetic operations (+,*,/,-), squareroot
and I think binary-decimal conversion. IEEE 754 says nothing about the trigonometric
functions. IEEE 754 specifies that for these operations, the result (using
the default rounding) must be the nearest representable real. (No, I can't be
bothered to specify the fine details). So that means the results of IEEE 754
operations must be correct within 0.5UDP, where UDP is the difference between
adjacent reals in your representation.
However for the transcendental functions: sin/cos/tan/exp/log/arctan/arcsin/arccos/
etcetera, virtually no-one guarantees this. In fact it is extremely difficult to
do, especially for sin/cos/tan. (I think I heard of a library which did it for
exp and log.) Broadly, you can expect to get a result within 1UDP in most cases.
But there are two serious (and related) problems:
(1) computing sin/cos/tan of extremely large arguments is quite difficult and
also quite hard. For example consider computing tan(10^100). This is
at least as hard as working out 10^100 (mod 2pi). Think about this and you'll
see you need multiprecision arithmetic somewhere. Some libraries actually do
this, some don't. (Try finding 10^100 on a PC and you might well get an answer
of 0, or 10^100.)
(2) For sin/cos/tan, it is awkward when the true result just happens to be near
0 but the argument is large (like 2pi). Because of the nature of
floating point, the difference between adjacent reals for the result will be
much smaller than the difference between adjacent reals for the argument.
It doesn't really make much sense to go for within 1UDP here because it means
doing extra multiprecision when that's not justified by the precision of the
argument.
But ignoring this, lets suppose we have one of the few really good libraries that
actually does do the multiprecision and so on, and guarantees results within
1UDP, or even 0.7UDP (You are unlikely to get much better than this without
really really hard work; 0.5UDP is of course the theoretical limit). Consider
what happens when you evaluate
sin(pi)
There are two approximations here. The first is of pi, which is represented with
a possible error of 0.5UDP. The second is of sin, which may be wrong to within
0.7UDP. So (true sin)(approximate pi) could be up to 0.5UDP. (approximate sin)
(approximate pi) could be even more inaccurate. Without doing the calculations,
it is easy to see that both Hugs and GHC are accurate like this.
Now I think the original suggestor wanted the library to spot that the answer
is "very close to" 0 and actually replace it with 0. Absolutely not! For example,
suppose I am trying to approximate the derivative of sin(x) near x=pi, by
computing sin(x+epsilon) for very small epsilon, then the last thing I need is
for the graph of the computed function to look like
\
\
\
------
\
\
\
! All sorts of physical calculations would be upset. It's really important
for functions that not only should they be close to the real function, but that
their local behaviour should be similar to that of the real function.
This problem also arises when you compute transcendental functions by doing it
for small sections and then gluing the sections together, since you have to be
very careful that at the joins, the functions still look reasonable; for example if
the true function is monotonic there, the approximation should be too.
