0

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?

Jasmeru
  • 113
  • 2
  • 1
    By default floating points are 4 byte, i.e. `277./14336.` is doing the divison of two 4 byte floats and the result is saved in a `real*8` which is on most systems a 8 byte float. – jack Jan 08 '21 at 22:24
  • 1
    [This answer](https://scicomp.stackexchange.com/a/973) explains how to define literals of specific precision. – jack Jan 08 '21 at 22:27
  • 1
    Does this answer your question? [precision of real variable](https://stackoverflow.com/questions/16992131/precision-of-real-variable) – jack Jan 08 '21 at 22:30
  • *First I thought it's because Fortran was running in 32 bits while Python ran in 64 bits.* The Fortran result is a 32-bit precision result. If you use `np.float32` in Python you'll see same. – wim Jan 08 '21 at 22:30
  • Note this is program is not Fortran77 whatever your advisor might claim - half the lines are not valid Fortran77 (and real*8 is not valid in any version of standard Fortran), and the source format is also not valid before Fortran 90 – Ian Bush Jan 08 '21 at 23:02
  • @quamrana: No, that question is not appropriate. Please do not promiscuously mark floating-point questions as duplicates of that one. It hinders providing specific answers to other issues in floating-point arithmetic. – Eric Postpischil Jan 09 '21 at 01:02

1 Answers1

3

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
Ian Bush
  • 6,996
  • 1
  • 21
  • 27
Barmar
  • 741,623
  • 53
  • 500
  • 612
  • 2
    Both answers are correct - which is possible as they are answers to different questions namely "what is the result of this division using default precision in Fortran", and "what is the result of this division using the default precision in python" – Ian Bush Jan 08 '21 at 22:39
  • Changed it to "more accurate". – Barmar Jan 08 '21 at 22:40
  • Is there a particular reason to introduce intermediate variables with conversion from integer values instead of using "64-bit" constants? – francescalus Jan 09 '21 at 01:37
  • @francescalus Didn't know about them, I just came up with something that should work. – Barmar Jan 09 '21 at 06:25
  • 1
    On a side note: `real*8` is non-standard. Could be better written as `use iso_fortran_env, only: real64; real(real64) :: x, y, z` – Scientist Jan 09 '21 at 20:55
  • 1
    I've edited the code part to bring it inline with the standard and show that is gives the desired result – Ian Bush Jan 10 '21 at 11:35
  • @IanBush As you made the edit, I'm curious: is the naming (wp), use of whitespace (indentation and inside parens) and of uppercase in keywords a part of a specific coding convention? –  Jan 10 '21 at 16:31
  • 1
    @Jean-ClaudeArbaut Not to my knowledge. wp for "working precision" seems reasonably common, but my code is more spaced out than most, and the capitalisation is just what emacs gives me. – Ian Bush Jan 10 '21 at 17:28