INT gives me incorrect value

[Amin_B_]

Hello

I have this code:

program th
  implicit none

  integer N1
  integer maxi,ei,Nc,ns,na
  real CH1,CH2



  OPEN(unit=1,file='input_file',status="old")
  read(1,*) ns                      !!!
  read(1,*) ei                         !!!!!!!!!!!!!!! 
  read(1,*) maxi!!!!!!!!!!!!!!!!!!!    
  read(1,*) N1!!!!!!!!!!!!!!!!  
  close(unit=1)


  CH1 = 0.07
  CH2 = -0.35


  Na = INT(abs(2.*((N1/2)*CH1 + (N1/2)*CH2)))
  write(*,*) Na,abs(2.*((real(N1)/2.)*CH1 + (real(N1)/2.)*CH2));stop
end program  th

and the input file is

1                      !!!!!!!!!!!
1                          !!!!!!!!!! 
1                   !!!   
1600 

Then I compile it with

ifort -O3 -autodouble t1.f90 -o out

but when I execute it I get 447 for na which is not correct. The correct answer is 448

Could you please tell me what is happening? and how I can fix it?

Thanks alot

[jimdempseyatthecove]

The problem you are experiencing is the intermediary components are not completely representable in floating point. (require higher/infinite precision). When taking the INT, you are performing a truncation of those rounding difference bits. To correct for this, you have to determine (specify) a rounding factor

Na = INT(abs(2.*((N1/2)*CH1 + (N1/2)*CH2))+0.5) ! round up within 0.5

or

Na = abs(2.*((N1/2)*CH1 + (N1/2)*CH2))
Na = INT(Na + Na * EPSILON(Na)) ! next larger number within the precision of Na

or

(your choice of rounding method)

Caution, the second method requires that the mathematical integer value of Na be fully representable with in the mantissa portion of the REAL(4) or REAL(8) number less one bit. In your example it does. In the case where it does not, the rounding is effectively adding 0.0).

Please keep in mind that numbers like 0.1 have no exact representation in floating point. Programs such as Calculator return results with the equivalent of the rounding expressed above. You as the programmer have to determine if you want this or do not want this. If you want it then you must add the type of rounding you prefer.

Jim Dempsey

[mecej4]

In general, regardless of whether you are using single, double or even quadruple floating point precision, you should not expect the same results from binary floating point arithmetic as from decimal arithmetic when doing the rounding or truncation of values whose fractional part in decimal notation is exactly half, i.e., 0.5.

Consider the expression (N1/2)*CH1. CH1 is 0.35, which cannot be represented exactly in a floating point register. Thus, the register value could be thought of as 0.349999999... which, when multiplied by (N1/2)=800, gives 279.9999999... instead of 280. If this number is truncated, you would have the integer 279 instead of 280. Apply this kind of reasoning to the rest of your expression, and see what may happen.

If you want exact decimal arithmetic, use a decimal arithmetic package. For a real-world example of mistakes of this nature, please see https://en.wikipedia.org/wiki/Vancouver_Stock_Exchange .

[TimP]

Nint may be useful.

[jimdempseyatthecove]

In response to your private message:

program INT_test
  implicit none

  integer N1
  integer maxi,ei,Nc,ns,na
  real CH1,CH2,temp
  OPEN(unit=1,file='input_file.dat',status="old")
  read(1,*) ns                      !!!
  read(1,*) ei                         !!!!!!!!!!!!!!! 
  read(1,*) maxi!!!!!!!!!!!!!!!!!!!    
  read(1,*) N1!!!!!!!!!!!!!!!!  
  close(unit=1)
  CH1 = 0.07
  CH2 = -0.35
  Na = INT(abs(2.*((N1/2)*CH1 + (N1/2)*CH2)))
  write(*,*) Na,abs(2.*((real(N1)/2.)*CH1 + (real(N1)/2.)*CH2)), "as-was"
  Na = INT(abs(2.*((N1/2)*CH1 + (N1/2)*CH2)) + 0.5)
  write(*,*) Na,abs(2.*((real(N1)/2.)*CH1 + (real(N1)/2.)*CH2)), "+.05"
  temp = abs(2.*((N1/2)*CH1 + (N1/2)*CH2))
  Na = INT(temp + temp*EPSILON(temp))
  write(*,*) Na,abs(2.*((real(N1)/2.)*CH1 + (real(N1)/2.)*CH2)), "+EPSILON"
  stop
end program  INT_test

If you compile a Debug build and break at the first Na =, then look at the Hex values for CH1 and CH2, you will find:

#BFD6666666666666 and #3FB1EB851EB851EC respectively

The CH1 mathematical value will have an infinite number of trailing 666's (it got rounded down because 6 is less than 8)
The CH2 mathematical value will end in an infinite sequence3 of EB851's (#3FB1EB851EB851EB851EB851EB851...) (it got rounded up because hex B is .GE. 8).

Also note that WRITE(*,*) will additionally perform rounding to produce a decimal representation. The numbers you see after WRITE are not necessarily the numbers used internally.

Jim Dempsey

[jimdempseyatthecove]

Additional note.

When compiling in Debug build on 32-bit platform, the default behavior for floating point arithmetic is to use the FPU (internal 80-bit format) as opposed to the SSE/AVX (64-bit with autodouble) performed in optimized build. The INT(expression) would have been evaluated using the larger internal temporaries in the Debug build, and the shorter 64-bit registers in the Release build). You can force the use of SSE/AVX in Debug build to produce more consistent results (at times they will differ due to rearranging expressions in optimization).

Jim Dempsey

[FortranFan]

@Amin,

If you could, can you please provide a bit of your computing background and what you plan to do with Fortran and how much exposure you have to the language capabilities in Fortran?  The reason for this question is your use of so-called "mixed mode arithmetic" in the code example of your original post.  See this link: http://docs.oracle.com/cd/E19957-01/805-4939/z400073a2265/index.html.  

As mecej4 and Jim have alerted you, "mixed mode arithmetic" should preferably be avoided in Fortran.  But it is quite acceptable in many computing environments; in fact, "mixed mode arithmetic" can bring more clarity in most of these other environments, albeit rounding error is something one has to keep in mind in any floating point calculation.

So I was wondering about your background and your needs to use Fortran given that if you want to use Fortran effectively, you will need to delve into Fortran-specific details.

Thanks,

[jimdempseyatthecove]

Amin,

>>  Na = INT(abs(2.*((N1/2)*CH1 + (N1/2)*CH2)))
 

You would perform any of your epsilon code prior to converting to INT

TEMP = abs(2.*((REAL(N1)/2.)*CH1 + (REAL(N1)/2.)*CH2)))
! now fix up TEMP
...
! then
Na = INT(TEMP) ! after you fixed TEMP

Jim Dempsey