Intel Fortran expression returns incorrect result

Question

I have a code as following:

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.

Yes, that should be 448. What happens if you removed the INT( )? It is an Integer anyway. How about if you compile it without the -autodouble? — Dijkgraaf, Sep 29 '15 at 01:52
Yes, I'm not saying remove it permanently, but just try it for this bit of code on it's own. It may be some sort of binary imprecision bug. Trying those things might help narrow it down. — Dijkgraaf, Sep 29 '15 at 01:55
I did that, without autodouble works correctly, do you have any idea what is going on? — amin bk, Sep 29 '15 at 02:06
Some numbers can't be represented exactly in binary (see http://stackoverflow.com/questions/1089018/why-cant-decimal-numbers-be-represented-exactly-in-binary), this can cause binary imprecision bugs. Try 1) breaking your Na calculation into smaller part e.g. (N1/2)*CH1 on its own etc. or possibly better 2) change it to N1*CH1/2.0 do multiplication before division and the same for the other part, change it to N1*CH2/2.0 — Dijkgraaf, Sep 29 '15 at 02:16
the actual "correct" answer according to my hp is 448. Not going to get into an edit war over it though. In any case `NINT` will probably give what you expect. — agentp, Sep 29 '15 at 02:35
I cannot reproduce this error with ifort 16. With or without your flags, ifort 16 gives Na = 448. — casey, Sep 29 '15 at 04:03
ifort14 precisely reproduces the behavior of the Question. I believe this is not a bug, because the argument of INT() becomes like 447.99999999999994316... (in other cases, the argument becomes 448.0000.... so the difference is only within floating-point roundoff). So anyway, as @agentp suggests, NINT() is recommended if the answer is known a priori to be integers. — roygvib, Sep 29 '15 at 15:07
@Ross Yes I figured though that was unintentional mistake as the rest of your edit I thought was good. — Dijkgraaf, Sep 29 '15 at 19:26

score 2 · Answer 1 · answered Sep 29 '15 at 06:04

This problem can be understood by investigating the numbers upto their full precision.

The number 0.07 is 7.0000000298023224E-02 in single precision and 7.0000000000000007E-02 in double precision. The number 800 can be represented exactly in both models.

However, the product of these numbers is approximately 56.000000238418579 and 56.000000000000006, respectively. Due to the available precision, the first gets rounded towards 56 (in single precision) and the second towards 56.000000000000007 (in double precision).

As a result, the calculation in single precision "gains" some precision by rounding in the "right" direction.

Concerning the different behaviour of ifort 16 mentioned by @casey, I guess there are some differences in rearranging the equation or maybe use of excess-precision for intermediate results. Though, this is only a guess.

When I remove input _file and put variable in the code, I get 448 correct value. Could you please tell why this happens? — amin bk, Sep 29 '15 at 11:21
I don't know. The assembler output (`-S`) should help us here, but I don't have a clue. I guess, that there are some optimizations during compilation which precalculate some results with different precision. Actually, there is now a difference between `(N1/2)*CH1 + (N1/2)*CH2` and `(N1/2)*(CH1+CH2)`. — Stefan, Sep 29 '15 at 12:58

Intel Fortran expression returns incorrect result

1 Answers1