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>
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)
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.
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
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.
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:
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.