Same number, different results

The following program creates a decimal-signifcand and a hexdecimal-significand number from an integer number. It then reads the two numbers as real64 values. It is expected that the two values are identical and that the program outputs nothing.

program test
  !
  use :: iso_fortran_env
  !
  implicit none
  !
  character(len=32) :: cdec, chex
  real(real64)      :: ddec, dhex
  integer(int64)    :: i, ib, ie
  !
  ib = shiftl(1_int64, 53)   ! 9'007'199'254'740'992 (54 bits)
  ie = ib + 15_int64
  !
  do i = ib, ie
    write (cdec, '(i0,a)') i, '.0'
    write (chex, '(a,z0,a)') '0x', i, 'p0'
    read (cdec, *) ddec
    read (chex, *) dhex
    if (ddec /= dhex) then
      write (*, '(a)') '----- Difference detected -----'
      write (*, '(2a,t46,es23.16,2x,z16.16)') &
            'Decimal-significand    : ', trim(cdec), ddec, ddec
      write (*, '(2a,t46,es23.16,2x,z16.16)') &
            'Hexadecimal-significand: ', trim(chex), dhex, dhex 
    end if
  end do
  !
end program test

Running the program on Windows 10 Pro, compiled with the Intel compilers (ifort: 2021.11.1 / ifx: 2024.0.2), results in this output:

----- Difference detected -----
Decimal-significand    : 9007199254740993.0   9.0071992547409940E+15  4340000000000001
Hexadecimal-significand: 0x20000000000001p0   9.0071992547409920E+15  4340000000000000
----- Difference detected -----
Decimal-significand    : 9007199254740997.0   9.0071992547409980E+15  4340000000000003
Hexadecimal-significand: 0x20000000000005p0   9.0071992547409960E+15  4340000000000002
----- Difference detected -----
Decimal-significand    : 9007199254741001.0   9.0071992547410020E+15  4340000000000005
Hexadecimal-significand: 0x20000000000009p0   9.0071992547410000E+15  4340000000000004
----- Difference detected -----
Decimal-significand    : 9007199254741005.0   9.0071992547410060E+15  4340000000000007
Hexadecimal-significand: 0x2000000000000Dp0   9.0071992547410040E+15  4340000000000006

It seems, the rounding is not correct for all values. I tried several compiler options (/standard-semantics, /fp:strict, /fltconsistency), without success.

Is this behaviour reproducible by others on Windows or on Unix/Linux?

Note:
This program/problem is Intel compiler specific.

1 Like

Welcome to the Fortran Discourse.

A 64-bit IEEE number use 52 bits (plus an implicit leading bit-52 set to 1). Your example, in contrast, uses a number that needs 54 bits to preserve exactness. Change

ib = shiftl(1_int64, 53)

to

ib = shiftl(1_int64, 52)

and run your example again.

Note also that when numbers (integer or real-float-double) are represented in strings and I/O is done with those numbers, errors generally get larger.

1 Like

I’ll also comment that what you’re seeing, in addition to what @mecej4 wrote, is “round to even”, which is the default rounding mode in the Intel compiler. Just to verify, I used the NAG compiler on this source:

integer, parameter :: dp = kind(0.0d0)
real(dp) :: a
a = real(Z'4340000000000001',dp)
write (*,'(F25.4)') a
a = real(Z'4340000000000002',dp)
write (*,'(F25.4)') a
end

and the output was:

    9007199254740994.0000
    9007199254740996.0000
1 Like

The point is that the two values in the output should be equal.
Lets take a closer look:

I think we can agree that the decimal-significand number “9007199254740993.0”
and the hexadecimal-significand number “0x20000000000001p0” are representations
of the same number. The rounding mode for all conversions should be “round to nearest, ties to even”. When we convert these two numbers to real64 values, I expect them to be equal.

This Python script converts both representations to doubles. I choosed Python
because its default rounding mode is also “round to nearest, ties to even”.

cdec = '9007199254740993.0'
chex = '0x20000000000001p0'
#
ddec = float(cdec)
dhex = float.fromhex(chex)
#
print("ddec         : %21.4f" % ddec)
print("dhex         : %21.4f" % dhex)
print("ddec == dhex : %s" % (ddec == dhex))

and the output was:

ddec         : 9007199254740992.0000
dhex         : 9007199254740992.0000
ddec == dhex : True

As we can see both numbers have been correctly rounded and they are equal.
Everything as expected. Now lets do the same with Fortran:

program test
  use :: iso_fortran_env
  implicit none
  !
  character(len=32) :: cdec, chex
  real(real64)      :: ddec, dhex
  !
  cdec = '9007199254740993.0'
  chex = '0x20000000000001p0'
  !
  read (cdec, *) ddec
  read (chex, *) dhex
  !
  write (*, '(a,f21.4)') 'ddec         : ', ddec
  write (*, '(a,f21.4)') 'dhex         : ', dhex
  write (*, '(a,l1)'   ) 'ddec == dhex : ', (ddec == dhex)
  !
end program test

and the output was:

ddec         : 9007199254740994.0000
dhex         : 9007199254740992.0000
ddec == dhex : F

Oops! Something unexpected happened. The decimal number has been rounded up, this is not correct. After experimenting with some more numbers, my conclusion is that the Intel compilers do “round to nearest” instaed of “round to nearest, ties to even”. I would say this is a bug.

1 Like

After doing some more tests, I agree with you. I note that NAG Fortran does the right thing.

program test
implicit none

integer, parameter :: dp = selected_real_kind(15)
integer, parameter :: i8 = selected_int_kind(15)
character(18) :: cval = "9007199254740993.0"

real(dp) :: dpx
integer(i8) :: i8x

write (*,*) "Source constant"
dpx = 9007199254740993.0_dp
i8x = transfer (dpx,i8x)
write (*,'(F30.5,1X,Z16)') dpx,i8x

write (*,*) "Internal read"
read (cval,'(F18.0)') dpx
i8x = transfer(dpx,i8x)
write (*,'(F30.5,1X,Z16)') dpx,i8x

end program test

Intel output:

 Source constant
        9007199254740994.00000 4340000000000001
 Internal read
        9007199254740994.00000 4340000000000001

NAG output:

 Source constant
        9007199254740992.00000 4340000000000000
 Internal read
        9007199254740992.00000 4340000000000000

Even if I explicitly put RN in the format, it does not round so that the LSB is zero. I will report this to Intel. Case number is 06097516.

1 Like

Thank you very much.

1 Like

Intel tells me that this is fixed for the forthcoming (I don’t know when) 2025.0 release.

1 Like

That’s good to hear, thanks for the update.

1 Like

I verified that this is fixed in the 2025.0 release, available now.

2 Likes

Thanks for keeping an eye on this case. I can confirm that the program from the first post now runs correctly. But unfortunately, the following, slightly modified program outputs again some differences:

program test
  !
  use :: iso_fortran_env
  !
  implicit none
  !
  character(len=32) :: cdec, chex
  real(real64)      :: ddec, dhex
  integer(int64)    :: i, ib, ie
  !
  ib = shiftl(1_int64, 52)   ! 4’503’599’627’370’496 (53 bits)
  ie = ib + 15_int64
  !
  do i = ib, ie
    ! Add 0.5 to the integer number.
    write (cdec, '(i0,a)') i, '.5'
    write (chex, '(a,z0,a)') '0x', i, '.8p0'
    read (cdec, *) ddec
    read (chex, *) dhex
    if (ddec /= dhex) then
      write (*, '(/a)') '----- Difference detected -----'
      write (*, '(2a,t50,es23.16,2x,z16.16)') &
            'Decimal-significand    : ', trim(cdec), ddec, ddec
      write (*, '(2a,t50,es23.16,2x,z16.16)') &
            'Hexadecimal-significand: ', trim(chex), dhex, dhex 
    end if
  end do
  !
end program test

Running the program on Windows 10 Pro, compiled with the Intel compiler ifx, version 2025.0.1, results in this output:

----- Difference detected -----
Decimal-significand    : 4503599627370496.5       4.5035996273704970E+15  4330000000000001
Hexadecimal-significand: 0x10000000000000.8p0     4.5035996273704960E+15  4330000000000000

----- Difference detected -----
Decimal-significand    : 4503599627370498.5       4.5035996273704990E+15  4330000000000003
Hexadecimal-significand: 0x10000000000002.8p0     4.5035996273704980E+15  4330000000000002

----- Difference detected -----
Decimal-significand    : 4503599627370500.5       4.5035996273705010E+15  4330000000000005
Hexadecimal-significand: 0x10000000000004.8p0     4.5035996273705000E+15  4330000000000004

----- Difference detected -----
Decimal-significand    : 4503599627370502.5       4.5035996273705030E+15  4330000000000007
Hexadecimal-significand: 0x10000000000006.8p0     4.5035996273705020E+15  4330000000000006

----- Difference detected -----
Decimal-significand    : 4503599627370504.5       4.5035996273705050E+15  4330000000000009
Hexadecimal-significand: 0x10000000000008.8p0     4.5035996273705040E+15  4330000000000008

----- Difference detected -----
Decimal-significand    : 4503599627370506.5       4.5035996273705070E+15  433000000000000B
Hexadecimal-significand: 0x1000000000000A.8p0     4.5035996273705060E+15  433000000000000A

----- Difference detected -----
Decimal-significand    : 4503599627370508.5       4.5035996273705090E+15  433000000000000D
Hexadecimal-significand: 0x1000000000000C.8p0     4.5035996273705080E+15  433000000000000C

----- Difference detected -----
Decimal-significand    : 4503599627370510.5       4.5035996273705110E+15  433000000000000F
Hexadecimal-significand: 0x1000000000000E.8p0     4.5035996273705100E+15  433000000000000E

It seems that the hexadecimal numbers are correctly rounded, but the decimal numbers are not!

Another way to test the correct reading of double precision floating point numbers is to use the data files from this Github project: https://github.com/nigeltao/parse-number-fxx-test-data/. The following program reads a test data file and verifies the conversion. The conversion can be done by the C function strtod() or by internal read:

program verify_dbl
  !
  use :: iso_c_binding
  use :: iso_fortran_env
  !
  implicit none
  !
  interface
    function strtod(string, final) bind(C, name='strtod')
      use :: iso_c_binding
      implicit none
      real(c_double)        :: strtod
      character(c_char)     :: string(*)
      type(c_ptr), optional :: final
    end function strtod
  end interface
  !
  character(len=4096) :: buff
  character(len=512)  :: inp_file
  character(len=10)   :: cmode
  real(real64)        :: dval1, dval2
  integer             :: blen, istat, lun, n_tests, n_success, n_failed
  logical             :: use_strtod
  !
  ! Get the conversion mode.
  if (command_argument_count() < 1) stop 'Error: missing conversion mode.'
  call get_command_argument(1, cmode)
  if (cmode == 'internal') then
    use_strtod = .false.
  else if (cmode == 'strtod') then
    use_strtod = .true.
  else
    stop 'Error: invalid conversion mode, use internal or strtod.'
  end if
  !
  ! Open the input file.
  if (command_argument_count() < 2) stop 'Error: missing input filename.'
  call get_command_argument(2, inp_file)
  open (newunit = lun,            &
        file    = trim(inp_file), &
        form    = 'formatted',    &
        status  = 'old'           )
  !
  ! Perform conversion tests.
  n_tests = 0; n_success = 0; n_failed = 0
  do
    read (lun, '(14x,z16,1x,a)', iostat=istat) dval1, buff
    if (istat /= 0) exit
    blen = len_trim(buff)
    if (use_strtod) then
      buff(blen+1:blen+1) = c_null_char
      dval2 = strtod(buff)
    else
      read (buff(1:blen), *, iostat=istat) dval2
    end if
    n_tests = n_tests + 1
    if (dval1 == dval2) then
      n_success = n_success + 1
    else
      n_failed = n_failed + 1
    end if
  end do
  close (lun)
  !
  write (*, '(a,1x,a)'  ) 'Results for file:', trim(inp_file)
  write (*, '(a,1x,a)'  ) 'Compiler version:', compiler_version()
  write (*, '(a,1x,a)'  ) 'Compiler options:', compiler_options()
  write (*, '(a,1x,a)'  ) 'Conversion mode :', trim(cmode)
  write (*, '(a,1x,i10)') 'Success tests   :', n_success
  write (*, '(a,1x,i10)') 'Failed  tests   :', n_failed
  write (*, '(a,1x,i10)') 'Total   tests   :', n_tests
  !
end program verify_dbl

Running this program with the conversion mode “strtod” and the data file “remyoudompheng-fptest-0.txt” outputs the following:

$ verify_dbl strtod remyoudompheng-fptest-0.txt
Results for file: .remyoudompheng-fptest-0.txt
Compiler version: Intel(R) Fortran Compiler for applications running on Intel(R) 64, Version 2025.0.1 Build 20241113
Compiler options: /warn:all /O2
Conversion mode : strtod
Success tests   :    1000000
Failed  tests   :          0
Total   tests   :    1000000

As you can see, all tests passed! Running the test with conversion mode “internal”, gives the following output:

$ verify_dbl internal remyoudompheng-fptest-0.txt
Results for file: remyoudompheng-fptest-0.txt
Compiler version: Intel(R) Fortran Compiler for applications running on Intel(R) 64, Version 2025.0.1 Build 20241113
Compiler options: /warn:all /O2
Conversion mode : internal
Success tests   :     502582
Failed  tests   :     497418
Total   tests   :    1000000

A lot of tests failed! I think, it shouldn’t matter if the conversion is done by the “strtod()” C function or by the internal Fortran read statement.

The test data files can be found in the “data” folder of the Github project.

It seems, Intel fixed only a part of the problem.