Confusing difference

The trivial code below produces what I consider inconsistent results. Specifically there is a difference between results when the mod intrinsic function is presented with a plain argument (c) or with an argument which requires evaluation (a+b). The code works exactly as expected when compiled with gfortran on Linux, but gives confusing results with Silverfrost Fortran. I am using version 8.82.0

  program main
  double precision :: a,b,c,d

  logical :: f1, f2
  a = 1.0d0/3.0d0
  b = 2.0d0/3.0d0
  write(*,'(a,e15.10e1)') 'a = ',a
  write(*,'(a,e15.10e1)') 'b = ',b
  c = a + b
  write(*,'(a,e15.10e1)') 'a + b = ',a+b
  write(*,'(a,e15.10e1)') 'c     = ',c
  f1 = (c > 1.0d0)
  f2 = (c < 1.0d0)
  write(*,*) 'c is smaller than 1.0 = ',f2
  write(*,*) 'c is larger  than 1.0 = ',f1
  f1 = (a+b > 1.0d0)
  f2 = (a+b < 1.0d0)
  write(*,*) 'a+b is smaller than 1.0 = ',f2
  write(*,*) 'a+b is larger  than 1.0 = ',f1 
  c = abs(mod(c,1.0d0))
  d = abs(mod(a+b,1.0d0))
  write(*,'(a,e15.10e1)') 'abs(mod(c,1.0))   = ',c
  write(*,'(a,e15.10e1)') 'abs(mod(a+b,1.0)) = ',d
  end

I changed your program to:

program main
  use iso_fortran_env
   implicit none
   double precision :: a,b,c,d
   logical :: f1, f2

   write (*,*) 'Vern : ',compiler_version ()

    a = 1.0d0/3.0d0
    b = 2.0d0/3.0d0
    write(*,'(a,e18.12)') 'a = ',a
    write(*,'(a,e18.12)') 'b = ',b
    c = a + b
    write(*,'(a,e18.12)') 'a + b = ',a+b
    write(*,'(a,e18.12)') 'c = ',c
    f1 = (c > 1.0d0)
    f2 = (c < 1.0d0)
    write(*,*) 'c is smaller than 1.0 = ',f2
    write(*,*) 'c is larger than 1.0 = ',f1
    f1 = (a+b > 1.0d0)
    f2 = (a+b < 1.0d0)
    write(*,*) 'a+b is smaller than 1.0 = ',f2
    write(*,*) 'a+b is larger than 1.0 = ',f1
    c = abs(mod(c,1.0d0))
    d = abs(mod(a+b,1.0d0))
    write(*,'(a,e18.12)') 'abs(mod(c,1.0)) = ',c
    write(*,'(a,e18.12)') 'abs(mod(a+b,1.0)) = ',d
end

The results are:

 Vern : FTN95 v9.04.0
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
 c is smaller than 1.0 =   F
 c is larger than 1.0 =   F
 a+b is smaller than 1.0 =   F
 a+b is larger than 1.0 =   F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000E+00

 Vern : GCC version 14.1.0
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
 c is smaller than 1.0 =  F
 c is larger than 1.0 =  F
 a+b is smaller than 1.0 =  F
 a+b is larger than 1.0 =  F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000E+00

I can't see the confusion ?

perhaps you could use 1.0d0 in your mod functions ?

Running John's program with FTN95 v9.03.0.

The 32 bit complier returns:

 Vern : FTN95 v9.03.0
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
 c is smaller than 1.0 =   F
 c is larger than 1.0 =   F
 a+b is smaller than 1.0 =   T
 a+b is larger than 1.0 =   F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.100000000000E+01

The 64 bit complier returns:

 Vern : FTN95 v9.03.0
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
 c is smaller than 1.0 =   F
 c is larger than 1.0 =   F
 a+b is smaller than 1.0 =   F
 a+b is larger than 1.0 =   F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000E+00

Also

 Vern : GCC version 12.3.0
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
 c is smaller than 1.0 =  F
 c is larger than 1.0 =  F
 a+b is smaller than 1.0 =  F
 a+b is larger than 1.0 =  F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000E+00


 Vern :
 Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel
 (R) 64, Version 2021.4.0 Build 20210910_000000
a = 0.333333333333E+00
b = 0.666666666667E+00
a + b = 0.100000000000E+01
c = 0.100000000000E+01
 c is smaller than 1.0 =  F
 c is larger than 1.0 =  F
 a+b is smaller than 1.0 =  F
 a+b is larger than 1.0 =  F
abs(mod(c,1.0)) = 0.000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000E+00

OK,

A new program that demonstrates 32-bit has different round-off.

program main
  use iso_fortran_env
   implicit none
   double precision :: a,b,c,d, one = 1
   logical :: f1, f2

   write (*,*) 'Version : ',compiler_version ()
   write (*,*) 'Options : ',compiler_options ()

    a = 1.0d0/3.0d0
    b = 2.0d0/3.0d0
    write(*,11) 'a = ',a
    write(*,11) 'b = ',b
    c = a + b
    write(*,11) 'a + b = ',a+b
    write(*,11) 'c = ',c
    f1 = (c > one)
    f2 = (c < one)
    write(*,*) 'c is smaller than 1.0 = ',f2
    write(*,*) 'c is larger than 1.0 = ',f1
    f1 = (a+b > one)
    f2 = (a+b < one)
    write(*,*) 'a+b is smaller than 1.0 = ',f2
    write(*,*) 'a+b is larger than 1.0 = ',f1
    c = abs(mod(c,one))
    d = abs(mod(a+b,one))
    write(*,11) 'abs(mod(c,1.0)) = ',c
    write(*,11) 'abs(mod(a+b,1.0)) = ',d
    write(*,11) 'one-(a+b) = ',one-(a+b)
    write(*,11) '(one-a)-b = ',(one-a)-b
    write(*,11) 'one-a-b   = ',one-a-b
 11 format (a,e24.18)
end

Results for revised program:

 Version : FTN95 v9.04.0
 Options : CHECKMATE;ECHO_OPTIONS;ERROR_NUMBERS;IMPLICIT_NONE;INTL;LINK;LOGL;NO_BANNER;UNLIMITED_ERRORS;VS7;
a = 0.333333333333333315E+00
b = 0.666666666666666630E+00
a + b = 0.100000000000000000E+01
c = 0.100000000000000000E+01
 c is smaller than 1.0 =   F
 c is larger than 1.0 =   F
 a+b is smaller than 1.0 =   T
 a+b is larger than 1.0 =   F
abs(mod(c,1.0)) = 0.000000000000000000E+00
abs(mod(a+b,1.0)) = 0.100000000000000000E+01
one-(a+b) = 0.555111512312578271E-16
(one-a)-b = 0.555111512312578271E-16
one-a-b   = 0.555111512312578271E-16

 Version : FTN95 v9.04.0
 Options : 64;CHECKMATE;ECHO_OPTIONS;ERROR_NUMBERS;IMPLICIT_NONE;INTL;LINK;LOGL;NO_BANNER;UNLIMITED_ERRORS;VS7;
a = 0.333333333333333303E+00
b = 0.666666666666666606E+00
a + b = 0.100000000000000000E+01
c = 0.100000000000000000E+01
 c is smaller than 1.0 =   F
 c is larger than 1.0 =   F
 a+b is smaller than 1.0 =   F
 a+b is larger than 1.0 =   F
abs(mod(c,1.0)) = 0.000000000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000000000E+00
one-(a+b) = 0.000000000000000000E+00
(one-a)-b = 0.111022302462515654E-15
one-a-b   = 0.111022302462515654E-15

 Version : GCC version 14.1.0
 Options : -cpp -iprefix C:/Program Files (x86)/gcc_eq/gcc_14.1.0/bin/../lib/gcc/x86_64-w64-mingw32/14.1.0/ -U_REENTRANT -mtune=generic -march=x86-64 -fno-underscoring -fdollar-ok
a = 0.333333333333333315E+00
b = 0.666666666666666630E+00
a + b = 0.100000000000000000E+01
c = 0.100000000000000000E+01
 c is smaller than 1.0 =  F
 c is larger than 1.0 =  F
 a+b is smaller than 1.0 =  F
 a+b is larger than 1.0 =  F
abs(mod(c,1.0)) = 0.000000000000000000E+00
abs(mod(a+b,1.0)) = 0.000000000000000000E+00
one-(a+b) = 0.000000000000000000E+00
(one-a)-b = 0.111022302462515654E-15
one-a-b   = 0.111022302462515654E-15

The differing results demonstrate that these type of tests should account for round-off uncertainty.

I am really not sure if 32-bit uses 80-bit registers for f2 = (a+b < one) ? c (= a+b) would be rounded to a 64-bit real.

For maximum precision you can use real*10 with Win32 rather than x64.

options(alt_kinds)
program main
real*10 :: a,b,c,d
logical :: f1, f2
a = 1.0_10/3.0_10
b = 2.0_10/3.0_10
write(*,'(a,e22.17e1)') 'a = ',a
write(*,'(a,e22.17e1)') 'b = ',b
c = a + b
write(*,'(a,e22.17e1)') 'a + b = ',a+b
write(*,'(a,e22.17e1)') 'c = ',c
f1 = (c > 1.0_10)
f2 = (c < 1.0_10)
write(*,*) 'c is smaller than 1.0 = ',f2
write(*,*) 'c is larger than 1.0 = ',f1
f1 = (a+b > 1.0_10)
f2 = (a+b < 1.0_10)
write(*,*) 'a+b is smaller than 1.0 = ',f2
write(*,*) 'a+b is larger than 1.0 = ',f1
c = abs(mod(c,1.0_10))
d = abs(mod(a+b,1.0_10))
write(*,'(a,e22.17e1)') 'abs(mod(c,1.0)) = ',c
write(*,'(a,e22.17e1)') 'abs(mod(a+b,1.0)) = ',d
end

a and b are evaluated at compile time whilst c and d are evaluated at runtime.

The main point is that the values of a and b are inevitably subject to round-off error.

Many thanks for all the comments. I kind of guessed that rounding errors were the cause of the discrepancy, but I didn't realise that 32bit executables were so subject to this issue.

I came across this issue when trying to understand why I was getting seriously differing results from a rather large program compiled under Windows & Linux. I normally have used 32 bit executables and in the past I havent noticed any great issues between the platforms

Any suggestions as the the best way forward? Just use 64 bit executables? Louis

64 bit executables are not generally the answer to issues relating to floating point precision.

For example 32 bit FTN95 provides for 80 bit floating point values as well as 32 and 64. Most Fortran compilers only provide 32 bit and 64 bit floating point values.

It is the programmer's responsibility to be aware/beware of calculations that are prone to round-off error.

This is a contrived example, but it illustrates Paul's point about 64 bit executables not being the 'fix all' answer to floating point issues.

program p
use iso_fortran_env
implicit none
integer, parameter :: n = 5
real*8 :: a(n) = [1.d0,1.d-16,1.d-16,-1.d0,1.d-16]
real*8 :: acc
integer :: i(n)
integer :: j
write (*,*) 'Vern :    ',compiler_version ()
write (*,*) 'Options : ',compiler_options ()
write(*,'(a)') 'Input array'
do j = 1, n, 1
  write(*,'(G16.10)') a(j)
end do
write(*,*)
write(*,'(a,G16.10)') 'Intrinsic sum = ', sum(a)
call dsort@(i,abs(a),5)   ! or qsortd from NSWC Library if not using FTN95
acc = 0.d0
do j = n, 1, -1
  acc = acc + a(i(j))
end do
write(*,*)
write(*,'(a,G16.10)') 'Sorted sum =    ',acc
end program  

32 bit executable:

 Vern :    FTN95 v9.03.0
 Options : ERROR_NUMBERS;LINK;NO_BANNER;UNLIMITED_ERRORS;VS7;
Input array
 1.000000000
0.1000000000E-15
0.1000000000E-15
-1.000000000
0.1000000000E-15

Intrinsic sum = 0.2999268806E-15

Sorted sum =    0.3000000000E-15

64 bit executable

 Vern :    FTN95 v9.03.0
 Options : 64;ERROR_NUMBERS;LINK;NO_BANNER;UNLIMITED_ERRORS;VS7;
Input array
 1.000000000
0.1000000000E-15
0.1000000000E-15
-1.000000000
0.1000000000E-15

Intrinsic sum = 0.1000000000E-15

Sorted sum =    0.3000000000E-15

I am curious why 32bit FTN95 still retain 80bit accuracy while 64bit one has only 64bit one? So both Intel and AMD processors have hardware for 80bit while 64bit FTN95 ignores it (as well as other Fortran compilers?)

Also how about SSE2 and AVX instructions with 64bit compilers. Do they still support 80bit accuracy? I've heard this opinion on the net. Also have heard that it's Windows which does not support double precision 80bit accuracy, Linux supports it (is this true?)

If anyone has some free time to investigate all that mess please try to make short confirming tests

64 bit applications can in theory access the the x87 FPU coprocessor which has 80 bit floating point precision. This coprocessor is separate from the the CPU that uses SSE etc. instructions and which basically does not extend to 80 bits.

The problem is that, at least under Windows, a massive time penalty is incurred when mixing x87 instructions with SSE etc. instructions. So compilers normally avoid this mix.

Quoted from DanRRight Also how about SSE2 and AVX instructions with 64bit compilers. Do they still support 80bit accuracy?

No ! this is the problem with AVX registers, they only support grouped 4/8 byte instructions.

As for 10-byte reals; I do not know if they are hardware supported in Intel or Gfortran x64 compilers. I think the answer could be No or maybe!

Given how simd instructions are so preferred, it is too inefficient to target 10-byte reals.

I think this thread is all about how 80-bit registers are used in the 32-bit FTN95 compiler. If the latest value is retained in an 80-bit register for the next calculation, then the round-off errors are much less.