2

I've been working on mixed C++/Fortran numerics code that needs to run on Windows and Linux and traced a discrepancy to the LOG10 function. I'm using gcc/gfortran on Linux and MinGW on Windows.

Here's an example:

PROGRAM FP
REAL VAL1, VAL2, ARG

DATA VAR1 / 12.5663710 /
DATA VAR2 / 10.6640625 /
DATA VAR3 / 1.08791232 /

ARG = VAR1 * VAR2 / VAR3
VAL1 = LOG10 (VAR1 * VAR2 / VAR3)
VAL2 = LOG10 (ARG)

WRITE (*,"(F30.25)") ARG
WRITE (*,"(F30.25)") LOG10(ARG)
WRITE (*,"(F30.25)") VAL1
WRITE (*,"(F30.25)") VAL2

END PROGRAM FP

On Linux, I get:

 123.1795578002929687500000000
   2.0905385017395019531250000
   2.0905385017395019531250000
   2.0905385017395019531250000

On Windows, I get

 123.1795578002929687500000000
   2.0905387401580810546875000
   2.0905387401580810546875000
   2.0905387401580810546875000

The same values are going into LOG10, but 2.09053850 is coming out on Linux and 2.09053874 on Windows. This is enough of a difference to cause substantial problems with testing. What can I do to get the same answer on both platforms?

I'm using someone else's Fortran code and am not an expert in its floating-point implementation details but found the problem by tracing the code side-by-side until the values diverged. The LOG10 seems to be the culprit.

As for compiler versions, on Linux I get:

$ gfortran --version
GNU Fortran (Ubuntu 9.2.1-9ubuntu2) 9.2.1 20191008

On Windows:

> gfortran --version
GNU Fortran (x86_64-posix-seh-rev0, Built by MinGW-W64 project) 8.1.0
Vladimir F Героям слава
  • 57,977
  • 4
  • 76
  • 119
Andreas Yankopolus
  • 852
  • 9
  • 22
  • Can be due to different default thread state. Not the same question, but related: https://stackoverflow.com/q/46028336/126995 – Soonts Nov 18 '19 at 23:42

2 Answers2

3

The difference comes by using a different runtime library. The log10 function implementation does not come from the libgfortran runtime library. Instead, on Linux, the standard C library GNU libc (GLIBC) is called https://www.gnu.org/software/libc/manual/html_node/Exponents-and-Logarithms.html (it is in the libm part).

Other compilers will do differently. Intel Fortran has its own runtime library and it actually gives the answer 2.0905387401580810546875000 on Linux.

On Windows it could depend on the GCC distribution you use. If you are using MinGW, then the Microsoft C runtime library is used instead of the GNU C library.

As far as I know, GLIBC is not available on Windows at all (Can I use glibc under windows?). You might try to take the log function from there and link it with your program but you will have to dig into the internals quite deeply.

Vladimir F Героям слава
  • 57,977
  • 4
  • 76
  • 119
  • This is fascinating and raises a whole bunch of questions on writing numerically consistent code across platforms. To further complicate things, My HP calculator gives 2.090538640692396153936101503375651, splitting the difference between the results. – Andreas Yankopolus Nov 19 '19 at 14:47
  • 1
    @AndreasYankopolus The calculator value will be in double or even higher precision. For higher numerical precision in Fortran use double precision and in limited occasions even quadruple precision maybe in order. The kind notation is the best way to go https://stackoverflow.com/questions/838310/fortran-90-kind-parameter I would expect there will still be some difference between GLIBC and Microsoft even in double precision though. Just smaller. – Vladimir F Героям слава Nov 19 '19 at 14:58
1

A consistent answer across platforms should be attained by promoting the argument of the log10 function to double precision and then converting the result back to single. Try the following function

`

real function mylog10(rarg)
integer, parameter :: dp = selected_real_kind(15,9)
real, intent(in) :: rarg
mylog10=log10(real(rarg,kind=dp))
end function mylog10

`

I get 2.0905387401580810546875000 for both Intel and gfortran on Linux. You could even make the function elemental if working on arrays.

RussF
  • 711
  • 3
  • 4
  • Thanks—this is largely how I ended up solving the problem. I didn't want to change the FORTRAN code for a variety of reasons, so I compiled it with `-fdefault-real-8` and did the float->double->float conversion in the calling C++ code. – Andreas Yankopolus Nov 20 '19 at 14:15