Long story short, I'm trying to rewrite in Python a Fortran77 code my advisor sent me, as Python is more convenient for me. While I was testing my code, I realize my outputs were slighlty different from the Fortran code. And it seems all of this stems from some rounding errors or whatnot in Fortran.
For instance, the fraction 277./14336. in Python returns:
print(277./14336.)
> 0.019321986607142856
^
But in Fortran77 I get:
program foo
implicit none
real*8 x
x=277./14336.
write(*,*) x
end program
> 1.9321987405419350E-002
^
So these numbers are equal up to the 8th significant digit, which in general should be good enough. But there are some finely tuned cancellation my code (precisely of order 1 part in 10^8) when I try to evaluate the numerical accuracy, so the error estimate from Fortran is sometimes twice as much as from the Python code.
What is going on? First I thought it's because Fortran was running in 32 bits while Python ran in 64 bits. But I got the same result when I ran a 32-bit version of Python (although I'm not sure it made a difference since I was still using a 64-bit OS) and I read that real*8 in Fortran means 8-byte precision, or 64 bits. Is there a fundamental difference in the representation of floating point numbers in Python and Fortran?
Python's answer is more accurate.
In Fortran you're assigning the result to a 64-bit float, but the inputs are 32-bit floats. Therefore, the division is done in 32-bit mode, and the result is then widened to 64 bits in the assignment.
Use 64-bit floats in the calculation and you should get the proper result.
$ cat div.f90
Program div
Use, Intrinsic :: iso_fortran_env, Only : wp => real64
Implicit None
Real( wp ) :: x
x = 277.0_wp / 14336.0_wp
Write( *, * ) x
End Program div
$ gfortran -Wall -Wextra -std=f2008 -fcheck=all div.f90
$ ./a.out
1.9321986607142856E-002