I have a FORTRAN program, where I assign PI=3.14159265
when I run the program, my numbers could out slightly off. So i recompiled with debug flags set, and started stepping through with gdb. A few lines in, I notice PI is not what i set it to?
(gdb) s
57 PI=3.14159265
(gdb) s
58 TWOPI= 2* PI
(gdb) print pi
$1 = 3.1415927410125732
What is going on here? I would understand if PI was a system value, but if it is, why is the value wrong?? how should i fix this?
It is printing the best binary approximation to the value of pi you entered; I get exactly the same answer in C using a float. You can get more accuracy by going to double precision, but this is as close as you can get to your value using binary with 32 bits and the usual representation of floating point numbers.
I'm not a fortan expert by any means, though F77 was the last version that I used. That does look about correct for the estimation which is single precision floating point representation. The short answer is that all numbers are represented in base 2. Single precision floating point representation for 3.14159265 ends up being 3.14159274... The wikipedia article for Single-Precision is surprisingly good for the subject. If you are interested, you should read Goldberg's "What Every Computer Scientist Should Know About Floating-Point Arithmetic".
Your example is actually represented as something pretty close to:
(1 + (4788187 * (2 ^ -23))) * 2
The reason is the rounding done for the representation as floating-point number.
You can define pi as parameter using
REAL, PARAMETER :: pi = 4. * ATAN(1.) ! set pi = 3.14...
to get a good representation of pi without worrying about the decimal to type in. This approach allows to increase the precession by changing the kind values of pi, 4. and 1.. It might help to get expected results from trigonometric functions, e.g. SIN(pi).
pi=acos(-1.0) is 2 bits more accurate than using atan.
