numeric equality between 10E10 and 100000000000

[psing51]

Hi all ,

I have a code in which numeric comparison against 10E10 is not working on intel platform , where as (as far as i can remember) the same code used to work on IBM platform. Here is the code (in the actual code, rdata is calculated & populated from various subroutines) - 

        real*8   rdata(21)
        rdata(1)=100000000000
        print *,"Values:",rdata(1),10E10
        if(rdata(1).eq.10E10) print *,"equal!"
        end

Now, after compiling the code with ifort , i get - 

user@machine1:~/TESTING/dft> ./a.out
 Values:   100000000000.000       9.9999998E+10
user@machine1:~/TESTING/dft>

as in output , i can see - 100000000000.000, hence i was expecting "equal" on the stdout

I tried printing contents of rdata using "gdb", and i got -

Missing separate debuginfos, use: zypper install libgcc_s1-debuginfo-6.2.1+r239768-2.4.x86_64
(gdb) n
3               print *,"Values:",rdata(1),10E10
(gdb) n
 Values:   100000000000.000       9.9999998E+10
4               if(rdata(1).eq.10E10) print *,"equal!"
(gdb) p rdata
$1 = (100000000000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
(gdb)

Is there a way/method to auto-convert  the 100000000000 to 100000000000.000 so as to achieve following effect (& achieve equality against 10E10). I made following modification in  code - 

        real*8   rdata(21)
        rdata(1)=100000000000.000
        print *,"Values:",rdata(1),10E10
        if(rdata(1).eq.10E10) print *,"equal!"
        end

Now, i can see following output-

user@machine1:~/TESTING/dft>
 ./a.out
 Values:   99999997952.0000       9.9999998E+10
 equal!
user@machine1:~/TESTING/dft>

and using gdb - 

Breakpoint 1, _unnamed_main$$ () at basic.f:2
2               rdata(1)=100000000000.000
Missing separate debuginfos, use: zypper install libgcc_s1-debuginfo-6.2.1+r239768-2.4.x86_64
(gdb) n
3               print *,"Values:",rdata(1),10E10
(gdb) n
 Values:   99999997952.0000       9.9999998E+10
4               if(rdata(1).eq.10E10) print *,"equal!"
(gdb) p rdata
$1 = (99999997952, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
(gdb)

[Juergen_R_R]

Hi Puneet,

in general, floating point numbers do not have exact representations in terms of byte in your computing architecture, see e.g. David Goldberg, "What every computer scientist should know about floating-point arithmetic". Particularly, floating point representations might be different if you change the hardware platform. Also note, that real*8 is not a portable floating point kind type and might  be different on different platforms. But my general answer would be that a test on exact equality for floating point numbers never really makes any sense. E.g. in our code we have a numerical utility module which defines an assert_equal routine that tests equality of floating point numbers up to a smallness parameter that is defined as a fudge factor (e.g. 1E3 or so) times epsilon(x) where x is one of the numbers to be compared. 

Cheers,

     JRR

[TimP · Accepted solution]

You may have used an old IBM compiler which didn't follow the Fortran standard about initializing default real constants to single precision. If you insist on avoiding f90, your double precision constant would be 10d10. There are also pitfalls in use of the ancient nonstandard real*8 notation in mixed precision, not that your example necessarily exhibits one, except that being nonstandard you might argue the standard doesn't apply to the constants.

[FortranFan]

@puneet s.,

As explained in the first comment by Juergen R. and as you will find in any reference on numerical computations, floating-point data comparisons are best done with a difference relative to a tolerance.  Depending on your calculations, you can decide on the tolerance you want to apply.  But you need to careful if you are working with legacy code or non-standard "FORTRAN" constructs (e.g., REAL*8) and across different platforms that the data are represented appropriately.  Toward this, you may want to review Dr Fortran blogs such as

https://software.intel.com/en-us/blogs/2013/12/30/doctor-fortran-in-its-a-modern-fortran-world

https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds

A quick illustration of how one can work defined KIND for floating-point representation and check for equality of data is as follows:

   ! WP below refers to a KIND parameter for use with Fortran REAL type
   ! that is effectively equivalent to IEEE "double" precision
   integer, parameter :: WP = selected_real_kind( p=15, r=307, radix=2 )

   real(kind=WP), parameter :: TOL_WP = epsilon(1.0_wp) ! a suitable tolerance

   real(kind=WP) :: x
   real(kind=WP) :: Calc_Result
   real(kind=WP) :: Expected_Result

   x = 10.0_wp
   Calc_Result = x**(sum([1,2,3,4])) ! silly simulation of a calculation
   Expected_Result = 1.0E10_wp

   if ( abs(Calc_Result-Expected_Result) <= TOL_WP ) then
      print *, "Difference between the calculated and expected result"
      print *, "is less than a tolerance of ", TOL_WP
   else
      print *, "Values differ"
   end if

end

Try running the above code on different platforms and see what you get.

[mecej4]

In addition to the good advice that has already been given, note this point: "10E10" is of type REAL, which occupies 32 bits, of which 23 bits (+ an implicit leading bit set to 1) are used for the significand (or mantissa). There is no exact representation possible for 10E10 in the IEEE single precision format. You may find it instructive to use an online utility to display various floating point numbers in decimal, binary and hexadecimal:

     https://babbage.cs.qc.cuny.edu/IEEE-754.old/Decimal.html .

[JVanB]

Looking at the problem from a low-level perspective:

program P
   implicit none
   integer i
   integer(selected_int_kind(18)) L
   real x
   real(selected_real_kind(15)) y
   L = int(10,kind(L))**11
   write(*,'(B64)') L
   write(*,'(*(i1))') (mod(i,10),i=63,0,-1)
   write(*,'(*(i1))') (i/10,i=63,0,-1)
   write(*,'()')
   i = shiftl(int(Z'3F8')+int(Z'008')*36,20)+ibits(L,13,23)
   x = transfer(i,x)
   write(*,'(B32)') i
   write(*,'(*(i1))') (mod(i,10),i=31,0,-1)
   write(*,'(*(i1))') (i/10,i=31,0,-1)
   write(*,'()')
   write(*,'(F13.1)') x
   i = 100000000000
   write(*,'(i0)') i
   y = L
BLOCK
   integer, parameter :: RK = KIND(y)
   if(y == 10E10_RK) write(*,'(a)') 'Equal'
END BLOCK
end program P

Output with ifort:

                           1011101001000011101101110100000000000
3210987654321098765432109876543210987654321098765432109876543210
6666555555555544444444443333333333222222222211111111110000000000

 1010001101110100100001110110111
10987654321098765432109876543210
33222222222211111111110000000000

99999997952.0
1215752192
Equal

The first line is 10**11 in binary. You can see that the highest set bit is bit 36 and the lowest is bit 11. Single precision REALs only carry 24 bits of precision, so that's bits 13:36. To make a REAL out of this, we remove the high bit because it's implicitly set in the REAL, extract the next 23 bits, bits 13:35, being careful to round bit 13 to nearest if necessary, and make those the low 23 bits of the REAL. Then the exponent is set to (exponent of 1.0)+36 and our REAL is complete. What it looks like in binary is the first short line of output.

Then you can see that it is in fact equal to 10**11-2**11 because the only set bit that was left behind on conversion was bit 11. Another thing I worry about is that 100000000000 is a default integer constant. Couldn't the compiler justifiably reduce it mod 2**32 to 1215752192 as would be the case if a default integer were assigned this value?

The last line of output shows the result if an appropriate KIND number is given for the constant 10E10. The syntax with the BLOCK construct allows us to extract the KIND number of variable y and use it in the literal 10E10_RK -- the KIND number has to be present in the literal because that is the place where the undesired rounding happens.