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.

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.

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

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.

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.

Thank you very much.