Floating point rounding issues

tsbg · ‎04-19-2023

Hi all

According to the Intel documentation, the default I/O rounding mode is 'PROCESSOR_DEFINED'. For ifort, this corresponds to 'NEAREST' which conforms to the ISO/IEC/IEEE 60559:2011 standard specification for 'round to nearest, ties to even'.

The following program shows a different behavior:

program round
  !
  implicit none
  !
  integer, parameter :: sp = kind(0.0)
  integer, parameter :: dp = kind(0.0d0)
  integer, parameter :: qp = selected_real_kind(33)
  !
  character(len=:), allocatable :: inp_sp, inp_dp, inp_qp
  real(sp)                      :: vio_dec_sp, vio_hex_sp, vas_dec_sp
  real(dp)                      :: vio_dec_dp, vio_hex_dp, vas_dec_dp
  real(qp)                      :: vio_dec_qp, vio_hex_qp, vas_dec_qp
  !
  write (*, '(/a/a)') '1. Single precision', &
                      '-------------------'
  ! The number 2**23 + 2**(-1) in decimal and hexadecimal:
  inp_sp = '8388608.5 0x1.000001p+23'
  read (inp_sp, *) vio_dec_sp, vio_hex_sp
  vas_dec_sp = 8388608.5_sp
  !
  write (*, '(a)') 'Values from I/O:'
  write (*, '(a,z8.8,2x,es16.9)') 'decimal    : ', vio_dec_sp, vio_dec_sp
  write (*, '(a,z8.8,2x,es16.9)') 'hexadecimal: ', vio_hex_sp, vio_hex_sp
  !
  write (*, '(/a)') 'Value from assignment:'
  write (*, '(a,z8.8,2x,es16.9)') 'decimal    : ', vas_dec_sp, vas_dec_sp
  !
  !
  write (*, '(/a/a)') '2. Double precision', &
                      '-------------------'
  ! The number 2**52 + 2**(-1) in decimal and hexadecimal:
  inp_dp = '4503599627370496.5 0x1.00000000000008p+52'
  read (inp_dp, *) vio_dec_dp, vio_hex_dp
  vas_dec_dp = 4503599627370496.5_dp
  !
  write (*, '(/a)') 'Values from I/O:'
  write (*, '(a,z16.16,2x,es23.16)') 'decimal    : ', vio_dec_dp, vio_dec_dp
  write (*, '(a,z16.16,2x,es23.16)') 'hexadecimal: ', vio_hex_dp, vio_hex_dp
  !
  write (*, '(/a)') 'Value from assignment:'
  write (*, '(a,z16.16,2x,es23.16)') 'decimal    : ', vas_dec_dp, vas_dec_dp
  !
  !
  write (*, '(/a/a)') '3. Quad precision', &
                      '-----------------'
  ! The number 2**112 + 2**(-1) in decimal and hexadecimal:
  inp_qp = '5192296858534827628530496329220096.5 &
           &0x1.00000000000000000000000000008p+112'
  read (inp_qp, *) vio_dec_qp, vio_hex_qp
  vas_dec_qp = 5192296858534827628530496329220096.5_qp
  !
  write (*, '(/a)') 'Values from I/O:'
  write (*, '(a,z32.32,2x,es40.33)') 'decimal    : ', vio_dec_qp, vio_dec_qp
  write (*, '(a,z32.32,2x,es40.33)') 'hexadecimal: ', vio_hex_qp, vio_hex_qp
  !
  write (*, '(/a)') 'Value from assignment:'
  write (*, '(a,z32.32,2x,es40.33)') 'decimal    : ', vas_dec_qp, vas_dec_qp
  !
end program round

Output from running the program on Windows 10 64bit with the ifort (2021.9.0) and the ifx (2023.1.0) compiler:

1. Single precision
-------------------
Values from I/O:
decimal    : 4B000000   8.388608000E+06
hexadecimal: 4B000000   8.388608000E+06

Value from assignment:
decimal    : 4B000001   8.388609000E+06

2. Double precision
-------------------

Values from I/O:
decimal    : 4330000000000001   4.5035996273704970E+15
hexadecimal: 4330000000000000   4.5035996273704960E+15

Value from assignment:
decimal    : 4330000000000001   4.5035996273704970E+15

3. Quad precision
-----------------

Values from I/O:
decimal    : 406F0000000000000000000000000001   5.192296858534827628530496329220097E+33
hexadecimal: 406F0000000000000000000000000000   5.192296858534827628530496329220096E+33

Value from assignment:
decimal    : 406F0000000000000000000000000001   5.192296858534827628530496329220097E+33

As you can see, the hexadecimal values are all correctly rounded. From the decimal values, only the single precision value is correctly rounded. The other precisions do not honor the 'ties to even' rule.

Another point is rounding on assignment. None of the assigned decimal values is correctly rounded according to the 'round to nearest, ties to even' rule.

I have not found any information about the rounding rules for assignment in the Fortran standard or in the Intel documentation. But I think the rule used by most compilers is 'round to nearest, ties to even'.

Questions:

Is this behavior intended, or is it eventually a bug?
Are there any compiler options to change the behavior?

Thanks for any insight.

--
Thomas

Ron_Green · ‎04-19-2023

I have not studied this carefully enough, so help me out.

I am a little skeptical of the internal reads using list-directed IO instead of explicit format edit descriptors.

Perhaps you can save me some time and show me what you EXPECTED in output. I assume it's 'round to nearest, ties to even' for the values:

8388608.5_sp

4503599627370496.5_dp

5192296858534827628530496329220096.5_qp

so you expected

8.388608000E+06

4.5035996273704960E+15

5.192296858534827628530496329220096E+33

for these 3 cases?

as an aside, gfortran 12.1 crashes on the internal read. perhaps you can report this to them.

gfortran -O0 -g wnew.f90
rwgreen@orcsle153:~/quad$ ./a.out

1. Single precision
-------------------
At line 18 of file wnew.f90
Fortran runtime error: Bad real number in item 2 of list input

Error termination. Backtrace:
#0  0x7fc2f4c23a12 in ???
#1  0x7fc2f4c24509 in ???
#2  0x7fc2f4c2505f in ???
#3  0x7fc2f4e62c73 in ???
#4  0x7fc2f4e6403d in ???
#5  0x7fc2f4e64c8a in ???
#6  0x4013f2 in round
	at /home/quad/wnew.f90:18
#7  0x402181 in main
	at /home/quad/wnew.f90:59

tsbg · ‎04-20-2023

Thanks for looking into this.

I am a little skeptical of the internal reads using list-directed IO instead of explicit format edit descriptors.

I changed the list-directed IO from

read (inp_sp, *) vio_dec_sp, vio_hex_sp
read (inp_dp, *) vio_dec_dp, vio_hex_dp
read (inp_qp, *) vio_dec_qp, vio_hex_qp

to explicit format edit descriptors

read (inp_sp, '(f9.0,1x,ex14.0)') vio_dec_sp, vio_hex_sp
read (inp_dp, '(f18.0,1x,ex22.0)') vio_dec_dp, vio_hex_dp
read (inp_qp, '(f36.0,1x,ex38.0)') vio_dec_qp, vio_hex_qp

but it didn't change anything about the results.

Perhaps you can save me some time and show me what you EXPECTED in output. I assume it's 'round to nearest, ties to even' for the values:
8388608.5_sp

4503599627370496.5_dp

5192296858534827628530496329220096.5_qp

so you expected

8.388608000E+06

4.5035996273704960E+15

5.192296858534827628530496329220096E+33

for these 3 cases?

Correct.

I expect that the output from the decimal floating point input values is the same as the output from the hexadecimal floating point input values which are correctly rounded according to the 'round to nearest, ties to even' rule.

as an aside, gfortran 12.1 crashes on the internal read. perhaps you can report this to them.

Sorry, I should have mentioned that the gfortran compiler does not yet support hexadecimal floating point numbers. I hope this will be the case in the upcoming version 13.

tsbg · ‎04-20-2023

Thanks for looking into this.

I am a little skeptical of the internal reads using list-directed IO instead of explicit format edit descriptors.

I changed the list-directed IO from

read (inp_sp, *) vio_dec_sp, vio_hex_sp
read (inp_dp, *) vio_dec_dp, vio_hex_dp
read (inp_qp, *) vio_dec_qp, vio_hex_qp

to explicit format edit descriptors

read (inp_sp, '(f9.0,1x,ex14.0)') vio_dec_sp, vio_hex_sp
read (inp_dp, '(f18.0,1x,ex22.0)') vio_dec_dp, vio_hex_dp
read (inp_qp, '(f36.0,1x,ex38.0)') vio_dec_qp, vio_hex_qp

but it didn't change anything about the results.

Perhaps you can save me some time and show me what you EXPECTED in output. I assume it's 'round to nearest, ties to even' for the values:
8388608.5_sp

4503599627370496.5_dp

5192296858534827628530496329220096.5_qp

so you expected

8.388608000E+06

4.5035996273704960E+15

5.192296858534827628530496329220096E+33

for these 3 cases?

Correct. I expect the output from the decimal floating point input values to be the same as the output from the hexadecimal floating point input values which are correctly rounded according to the 'round to nearest, ties to even' rule.

as an aside, gfortran 12.1 crashes on the internal read. perhaps you can report this to them.

Sorry, I should have mentioned that the gfortran compiler does not yet support hexadecimal floating point numbers. I hope this will be the case in the upcoming version 13.

tsbg · ‎04-20-2023

Thanks for looking into this.

I am a little skeptical of the internal reads using list-directed IO instead of explicit format edit descriptors.

I changed the list-directed IO from

read (inp_sp, *) vio_dec_sp, vio_hex_sp
read (inp_dp, *) vio_dec_dp, vio_hex_dp
read (inp_qp, *) vio_dec_qp, vio_hex_qp

to explicit format edit descriptors

read (inp_sp, '(f9.0,1x,ex14.0)') vio_dec_sp, vio_hex_sp
read (inp_dp, '(f18.0,1x,ex22.0)') vio_dec_dp, vio_hex_dp
read (inp_qp, '(f36.0,1x,ex38.0)') vio_dec_qp, vio_hex_qp

but it didn't change anything about the results.

Perhaps you can save me some time and show me what you EXPECTED in output. I assume it's 'round to nearest, ties to even' for the values:
8388608.5_sp

4503599627370496.5_dp

5192296858534827628530496329220096.5_qp

so you expected

8.388608000E+06

4.5035996273704960E+15

5.192296858534827628530496329220096E+33

for these 3 cases?

Correct. I expect the same output from the decimal floating point input values as from from the hexadecimal floating point input values which are correctly rounded according to the 'round to nearest, ties to even' rule.

as an aside, gfortran 12.1 crashes on the internal read. perhaps you can report this to them.

Sorry, I should have mentioned that the gfortran compiler does not yet support hexadecimal floating point numbers. I hope this will be the case in the upcoming version 13.

--

Thomas

tsbg · ‎04-24-2023

Please excuse the multiple replies. The first two replies where marked as spam and didn't showed up in the forum first. I observed that my Pi-Hole blocked some URLs referenced from the Intel forum. I put them on the whitelist and replied a third time and this time it was successful. Unfortunately the first two replies showed now up too.

Never less, I made a different example to show what I think is a problem:

program round_c
  !
  use, intrinsic :: iso_c_binding
  !
  implicit none
  !
  integer, parameter :: dp = kind(0.0d0)
  !
  interface
    function strtod(str, endptr) bind(c)
      use, intrinsic :: iso_c_binding
      implicit none
      real(c_double)         :: strtod
      character(kind=c_char) :: str(*)
      type(c_ptr), optional  :: endptr
    end function strtod
  end interface
  !
  character(len=:), allocatable :: inp
  real(dp)                      :: v1, v2
  !
  ! The number 2**52 + 2**(-1) in decimal:
  inp = '4503599627370496.5'
  !
  read (inp, '(f18.0)') v1
  !
  v2 = strtod(inp//c_null_char)
  !
  write (*, '(/a/a)') 'Results', &
                      '-------'
  write (*, '(a,z16.16,2x,es23.16)') 'internal read   : ', v1, v1
  write (*, '(a,z16.16,2x,es23.16)') 'C strtod()      : ', v2, v2
  write (*, '(a,a)') 'Are they equal? : ', trim(merge('Yes', 'No ', v1 == v2))
  !
end program round_c

Output from ifort (2021.9.0) and ifx (2023.1.0):

Results
-------
internal read   : 4330000000000001   4.5035996273704970E+15
C strtod()      : 4330000000000000   4.5035996273704960E+15
Are they equal? : No

Output from gfortran 12.2

Results
-------
internal read   : 4330000000000000   4.5035996273704960E+15
C strtod()      : 4330000000000000   4.5035996273704960E+15
Are they equal? : Yes

I would expect that the Intel compilers produce the same result as the gfortran compiler.

Can someone else reproduce these results?

Thanks

--

Thomas