[BUG] Float>>tan

Richard A. O'Keefe ok at atlas.otago.ac.nz
Thu Aug 23 01:04:07 UTC 2001


Ned Konz <ned at bike-nomad.com> wrote:
	"Float halfPi cos" evaluates to a very small number (the reciprocal of the 
	number above, about 1.225e-16 on my system). Should be 0.
	
No, it should NOT be zero.  The problem is that halfPi is NOT pi/2.
                pi
Recall that cos(-- + delta) = - sin(delta) approx= - delta.
                2
Since pi/2 is a transcendental number, it *cannot* be represented exactly
by a floating-point number.  There *must* be some error.

Float halfPi prints as    1.570796326794897
A better approximation is 1.57079632679489661923
so we see that the error is about 4E-16.

I just ran the following C program:
#include <math.h>
#include <stdio.h>
int main(void) {
    printf("sin: %.17e\n", sin(M_PI_2));
    printf("cos: %.17e\n", cos(M_PI_2));
    printf("tan: %.17e\n", tan(M_PI_2));
    return 0;
}

The output was
sin: 1.00000000000000000e+00
cos: 6.12323399573676604e-17
tan: 1.63312393531953700e+16

NB: this was on a SPARC/Solaris 2.8 system.  The trig library uses the
"infinite precision PI" technique.  What happens if the trig library uses
"machine precision PI" instead?  Well, to start with, the results get worse
and worse and worse as the arguments get bigger.  But *then* you might
expect Float halfPi cos to be zero.  PROVIDED that the "machine precision PI"
used by the math library was identical to the value of PI used by Float.
Float initializes
    Pi := 3.lots and lots of digits.
so we are dependent on the code that converts floating point numbers from
character strings to binary getting its rounding exactly right.

Squeak, at least as of 3.0, takes extra care to print numbers well.
It does not take quite such care when reading floating point numbers.
xxxx.yyyy is computed as
	(xxxx, an exact integer) asFloat +
     (  (yyyy, an exact Integer) asFloat /
        (10 raisedTo: number of digits in yyyy) asFloat )

In a workspace, do this.
    Type "Float halfPi  "
    PrintIt (Cmd-P, Alt-P, whatever)
    Insert a "-" between the two spaces you typed in the first step.
    Select the whole expression "Float halfPi - <number>".
    PrintIt

For me, the final result comes out as
	-4.44089209850063e-16

That is, subtracting halfPi from "itself" does not give zero!

I believe that the cause is roundoff error in Number>>readFromString:.

This means that the value of pi (hence also the value of halfPi) used
by Float is ever-so-slightly *wrong*.  Even if the host math library does
use "machine precision pi", the value stored in the HalfPi class variable
of Float is NOT THE SAME NUMBER.

For the last 50 years, savvy programmers who have wanted a good approximation
of pi in the programs have asked the math library to compute it, e.g.
	HALFPI = 2.0*ATAN(1.0)
in Fortran.  Decimal/binary conversion errors are not the least of the
reasons for doing this.





More information about the Squeak-dev mailing list