[BUG] Float>>tan
Richard A. O'Keefe
ok at atlas.otago.ac.nz
Wed Aug 22 23:54:05 UTC 2001
Antonio Barros (ajbn at cin.ufpe.br) wrote:
(Float halfPi) tan ==> 8.16588936419192e15
I don't see a problem here.
The Float>>tan implementation is written as "^self sin / self cos".
That is regrettable.
Let's see what's actually there.
Float>>tan ^self sin / self cos
Float>>cos ^(self + Halfpi) sin
Float>>sin <primitive: 56> ....
It's a good thing that the backup code for Float>>sin is not normally used;
the code for range reduction is pretty horrible, and the taylor series
part isn't up to modern standards. (IEEE arithmetic provides a remainder()
function which is expressly designed for jobs like range reduction, although
the trig functions really require a more subtle "infinite precision PI"
technique which is available in some C libraries.)
It might be an excellent idea to plug in primitives for cos and tan as well.
In this case, "self cos" should have evaluated to zero and we
should have got a division by zero. Am I missing something here?
Yes. Halfpi is *not* pi/2, it is pi/2 + epsilon, for some small but non-zero
epsilon. So what we really have here is
(pi/2 + epsilon) sin
-------------------
(pi + 2.epsilon) sin
which is going to be approximately 1/(2 epsilon). Since 64-bit
floating-point numbers have about 16 digits of precision, we expect
an answer somewhere in the vicinity of 10**16, and lo! that's what we get.
In short, you have forgotten roundoff error.
More information about the Squeak-dev
mailing list
|