1

I see a discrepancy between python and Fortran when using the sinus function. Could anyone shed light on this, please?

in python:

      import math
      print(math.sin(6.28318530717959))
      >> 3.3077843189710302e-15

in fortran90:

      print*, sin(6.28318530717959d0)
      >> 3.3077720792452914E-15

EDIT: As it seems to be a Fortran compiler issue, I used g95 with

       g95 -O3 test.f90 -o test.exe
Christian
  • 991
  • 2
  • 13
  • 25
  • 1
    probably different math libraries. – Jean-François Fabre Feb 25 '17 at 22:30
  • wolfram alpha is giving `3.523074713233440994231605661193961688584318794 × 10^-15` – juanpa.arrivillaga Feb 25 '17 at 22:38
  • Python uses the C math library. C also reports 3.3077843189710302386e-15. If Wolfram Mathematica is any standard of high-precision computations, it also reports the same answer. – DYZ Feb 25 '17 at 22:38
  • @juanpa.arrivillaga Are you sure? Because Wolfram Mathematica is in agreement with C and Python. – DYZ Feb 25 '17 at 22:39
  • 1
    clearly float32 used for Fortran library: 8 decimal digits only. – B. M. Feb 25 '17 at 22:40
  • @DYZ currently on my phone. I'm surprised that Wolfram Mathematica and Wolfram Alpha would give different results. – juanpa.arrivillaga Feb 25 '17 at 22:40
  • @juanpa.arrivillaga Isn't that odd... – DYZ Feb 25 '17 at 22:41
  • @B.M. using numpy.float32, and `np.sin`, I'm getting `1.7484555e-07` – juanpa.arrivillaga Feb 25 '17 at 22:44
  • My gfortran gives 3.3077843189710302E-015w – tim18 Feb 25 '17 at 22:51
  • @tim18 I used g95 as fortran compiler – Christian Feb 25 '17 at 22:54
  • 1
    My gfortran gives 3.3077843189710302E-015; ifort 3.307784318971030E-015; it seems there are only 2 significant digits correct due to cancellation with range reduction near 2pi. The Oracle compiler math libraries are supposed to have full precision reduction, or you could use gmp et al. (or C long double, Fortran REAL80 or REAL128) It might be more interesting if your argument were the nearest representable value to 2pi. 3.523e-15 seems a consensus. – tim18 Feb 25 '17 at 23:16
  • @B.M. clearly it is double precision. It is in the code (`d0`). I can't imagine how double precision could be 32bit on a modern compiler although it does not have to be the same as C double. – Vladimir F Героям слава Feb 25 '17 at 23:24
  • @Christian You should state which Fortran compiler you are using and which compiler flags. – Vladimir F Героям слава Feb 25 '17 at 23:25
  • gfortran-6.3 on macox10.11 gives 3.3077843189710302E-015, and Julia 0.5.0 gives 3.3077843189710302e-15 (so they are same as python). Because sin(x) takes a value on the order 1, it is probably not very strange to see such difference (with g95) depending on the algorithm/library used? – roygvib Feb 26 '17 at 05:12
  • 1
    Please tell us the compiler version! Also, the answer you've accepted is wrong. – Ross Feb 26 '17 at 09:08

1 Answers1

-1

According to IEEE 754 for float representation:

In [7]: bin(3.3077720792452914e-15.view(np.uint64))
Out[7]: '0b11110011101101110010110011010000000000000000000000000000000000'

shows a truncated mantissa, when 

In [9]: bin(3.3077843189710302e-15.view(np.uint64))
Out[9]: '0b11110011101101110010110011101100111001100111010111010001111111'

shows a plain one.

Probably a type issue, with a float32 in the process, even the origin is mysterious.

B. M.
  • 18,243
  • 2
  • 35
  • 54
  • 1
    Since Fortran 77, the double precision argument automatically promotes the intrinsic to double, so sin() and dsin() (Fortran 66 compatibility) must be the same. I don't know of any Fortran in the last decade where c_double could differ from double precision, but technically they are different data types even if sharing the same kind value. The question is a bit more uncertain when moving to c_long_double and the like (and x87 extended precision mode, not available in all compilers). – tim18 Feb 26 '17 at 00:45
  • 3
    This answer is wrong as mentioned by tim18. SIN is the generic name of all Fortran sinus functions ! The right sinus function is selected depending on the argument type. – Francois Jacq Feb 26 '17 at 08:53