Error in NINT(), 32-bit FTN95

FTN95 (32-bit, many versions up to 8.3) was giving me a result off by 1.0 when evaluating NINT(double precision expression). The numbers are such that each term in the expression, as well as the sum of the terms, is exactly representable in IEEE-64 format. I give below a cut-down reproducer program.

The error does not occur with /64, nor with FTN77 V 4.03. The value of the argument to NINT when evaluating TC1 is 8912.5000D0; in IEEE-64 expressed in hex, Z'40C1684000000000'. The incorrect result given is 8912. The correct result is 8913.

I realise that, depending on whether constant subexpressions are evaluated at compile time and those values used at run time, the results can be different, so that the value actually used may be 8912.499999999999 instead of 8912.5. However, the 32-bit debugger does not show the x87 registers, so I could not look into this aspect.

The Fortran standard defines NINT(x) as INT(x+0.5) for x > 0, so the usual 'nearest or even' rule of arithmetic does not apply.

      program xttbl
      implicit none
      double precision ti, tr, co2, co3, r1, r2, r3
      double precision e1, e2, e3, es, tc1, tc2, tc3
      integer ni
      data ti/77400.0d0/, co2/19100d0/, co3/77400d0/
      data r1/0.1d0/, r2/0.12d0/, r3/0.22d0/
C
      ni=int(ti*2d-2)
      tr=ni*50d0+25d0-co3
      e1=r3*tr
      e2=r2*(co3-co2)
      e3=r1*co2
      es=e1+e2+e3+1.0d0
      tc1=nint(r3*tr+r2*(co3-co2)+r1*co2+1.d0)
      tc2=nint(e1+e2+e3+1.0d0)
      tc3=nint(es)
      write(*,10)'e1,e2,e3,es  ',e1,e2,e3,es
      write(*,20)'TC1  ',tc1
      write(*,20)'TC2  ',tc2
      write(*,20)'TC3  ',tc3
   10 format(1x,A,2x,4F12.4)
   20 format(1x,A,14x,F7.0)
      end

mecej4

Many thanks for the feedback. I will make a note that this needs fixing.

I have reproduced the error using FTN95 Ver 8.4.0, although I am not sure the test proves there is something wrong. I changed the example to attempt to print the input to NINT, although if the value is in the x87 register, write may not show that value. My changed code is

      program xttbl 
      implicit none 
      double precision ti, tr, co2, co3, r1, r2, r3 
      double precision e1, e2, e3, es, tc1, tc2, tc3 
      integer ni 
      data ti/77400.0d0/, co2/19100d0/, co3/77400d0/ 
      data r1/0.1d0/, r2/0.12d0/, r3/0.22d0/ 
!C 
      ni=int(ti*2d-2) 
      tr=ni*50d0+25d0-co3 
      e1=r3*tr 
      e2=r2*(co3-co2) 
      e3=r1*co2 
      es=e1+e2+e3+1.0d0 
      tc1=nint(r3*tr+r2*(co3-co2)+r1*co2+1.d0) 
      tc2=nint(e1+e2+e3+1.0d0) 
      tc3=nint(es) 
      write(*,10)'e1,e2,e3,es  ',e1,e2,e3,es 
      write(*,*)'TC1  ',tc1,  (r3*tr+r2*(co3-co2)+r1*co2+1.d0) ! -tc1
      write(*,*)'TC2  ',tc2,  (e1+e2+e3+1.0d0) ! -tc2
      write(*,*)'TC3  ',tc3,  (es) ! -tc3
   10 format(1x,A,2x,4F12.4) 
   20 format(1x,A,14x,F7.0, es18.10) 
      end

The results I get with 32-bit are: [FTN95/Win32 Ver. 8.40.0 Copyright (c) Silverfrost Ltd 1993-2018]

[Current options] LGO;

0023)    20 format(1x,A,14x,F7.0, es18.10) 
WARNING - Label 20 is declared, but not used
    NO ERRORS, 1 WARNING  [<XTTBL> FTN95 v8.40.0]
Creating executable: c:\temp\forum\mecej4\lgotemp@.exe
Program entered
 e1,e2,e3,es          5.5000   6996.0000   1910.0000   8912.5000
 TC1            8912.00000000              8912.50000000    
 TC2            8913.00000000              8912.50000000    
 TC3            8913.00000000              8912.50000000    

The use of NINT for these values, where the result could potentially change based on round-off in the last bit, is not a reliable algorithm. You could test if nint(a) - a < 0 .or. nint(a) - a >= 1 ) then error condition but I would suspect this test is subject to round-off problems as well. Basically I think it is an unstable or unsuitable algorithm.

( In FTN95 write (,) does not show the full precision available, so your use of hex output may be required, something I am not familiar with decoding )

Thanks for your comments, John.

I agree that we are working with a corner case. However, these findings struck me as strange:

(i) X87 arithmetic (as used by programs compiled with FTN77, FTN95-32 bit) is in general slower but more accurate than SSE2 arithmetic (as generated by FTN95 /64). In this example, the opposite happened.

(ii) FTN77 and FTN95 (32-bit) target IA32+X87, yet in this example the results from the two sibling compilers are different.

Added to this is something that I found frustrating: the inability of SDBG to show X87 register values in hex, even after entering assembly mode. I tried running the EXE produced by FTN95 32-bit using two other debuggers, but they do not understand the debug symbols that Silverfrost compilers put out, and they even failed in some places to distinguish between code and data. I gave up trying to find an explanation for the behaviour, since one debugger (the VC debugger) showed register contents but not all the assembler code, and SDBG, which showed the assembler code but not the X87 register contents in hex.

If you ever wish to just display a real/real*8 number in IEEE hex format, you can use the utility at https://babbage.cs.qc.cuny.edu/IEEE-754/ . Just enter 8912.5 under 'Value to analyze' and press the Enter key, and observe that there are lots of trailing zeros under 'Binary64'.

For what it's worth, the same answer comes out if you use IDNINT instead of NINT.

I did a search of several codes I have written, and discovered that NINT is not something I use, as 'nearest integer' is not (always) the same as 'rounding up', which I do on many occasions. Perhaps the bug has not been apparent because NINT isn't used all that much.

Still ought to work properly, though.

Eddie

IDNINT is the specific name for the generic NINT when the argument is double precision. Either name produces the same inline code when FTN95 is used to compile for X87.

The argument, if an expression, will already be in an X87 register. The sign is checked and 0.5 is added or subtracted, as appropriate. The X87 control word is changed to truncate, an integer store to memory is performed and the control word restored (to round to nearest even).

Thus, it is no cause for surprise that changing from NINT to IDNINT produces identical behaviour.

My modified code does now show that the problem is with x87 calculations program xttbl implicit none double precision ti, tr, co2, co3, r1, r2, r3 double precision e1, e2, e3, es, tc1, tc2, tc3 real10 c1,c2,c3 integer ni integer8 jj data ti/77400.0d0/, co2/19100d0/, co3/77400d0/ data r1/0.1d0/, r2/0.12d0/, r3/0.22d0/ !C ni=int(ti2d-2) tr=ni50d0+25d0-co3 e1=r3tr e2=r2(co3-co2) e3=r1co2 es=e1+e2+e3+1.0d0 tc1=nint(r3tr+r2*(co3-co2)+r1co2+1.d0) tc2=nint(e1+e2+e3+1.0d0) tc3=nint(es) write(,10)'e1,e2,e3,es ',e1,e2,e3,es write(,)'TC1 ',tc1, (r3tr+r2(co3-co2)+r1co2+1.d0) ! -tc1 write(,)'TC2 ',tc2, (e1+e2+e3+1.0d0) ! -tc2 write(,*)'TC3 ',tc3, (es) ! -tc3 10 format(1x,A,2x,4F12.4)

      write(*,20)'TC1  ',tc1,  transfer(r3*tr+r2*(co3-co2)+r1*co2+1.d0,jj) ! -tc1
      write(*,20)'TC2  ',tc2,  transfer(e1+e2+e3+1.0d0,jj) ! -tc2
      write(*,20)'TC3  ',tc3,  transfer(es,jj) ! -tc3
   20 format(1x,A,14x,F7.0, z18.14) 

      c1=(r3*tr+r2*(co3-co2)+r1*co2+1.d0) 
      c2=(e1+e2+e3+1.0d0) 
      c3=(es) 
      tc1=nint(c1)
      tc2=nint(c2)
      tc3=nint(c3) 
      write(*,30)'TC1  ',tc1, c1
      write(*,30)'TC2  ',tc2, c2
      write(*,30)'TC3  ',tc3, c3
   30 format(1x,A,14x,F7.0, es26.18, z18.14) 

      end

This 'faulty' result seems to me to be a consequence of the fact that the expression 'r3tr+r2(co3-co2)+r1*co2+1.d0' is stored in an 80 bit register.

To get the 'correct' result this must be truncated to a 64 bit value before applying NINT.

I am not sure that it makes sense to 'fix' the compiler in this context. In a vast majority of cases the status quo should give more accurate results.

p.s. Edited. I should have used the word 'rounded' rather than 'truncated'.