forums.silverfrost.com Forum Index forums.silverfrost.com
Welcome to the Silverfrost forums
 
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

Error in NINT(), 32-bit FTN95

 
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Support
View previous topic :: View next topic  
Author Message
mecej4



Joined: 31 Oct 2006
Posts: 1085

PostPosted: Wed Feb 27, 2019 3:33 pm    Post subject: Error in NINT(), 32-bit FTN95 Reply with quote

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.

Code:
      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
Back to top
View user's profile Send private message
PaulLaidler
Site Admin


Joined: 21 Feb 2005
Posts: 5794
Location: Salford, UK

PostPosted: Thu Feb 28, 2019 10:11 am    Post subject: Reply with quote

mecej4

Many thanks for the feedback. I will make a note that this needs fixing.
Back to top
View user's profile Send private message
JohnCampbell



Joined: 16 Feb 2006
Posts: 2065
Location: Sydney

PostPosted: Thu Feb 28, 2019 11:34 am    Post subject: Reply with quote

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
Code:
      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:
Code:
[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 )
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1085

PostPosted: Thu Feb 28, 2019 12:53 pm    Post subject: Reply with quote

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".


Last edited by mecej4 on Sat Mar 02, 2019 3:01 am; edited 1 time in total
Back to top
View user's profile Send private message
LitusSaxonicum



Joined: 23 Aug 2005
Posts: 1960
Location: Yateley, Hants, UK

PostPosted: Fri Mar 01, 2019 1:40 pm    Post subject: Reply with quote

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
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1085

PostPosted: Sat Mar 02, 2019 3:00 am    Post subject: Reply with quote

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.
Back to top
View user's profile Send private message
JohnCampbell



Joined: 16 Feb 2006
Posts: 2065
Location: Sydney

PostPosted: Sat Mar 02, 2019 6:06 am    Post subject: Reply with quote

My modified code does now show that the problem is with x87 calculations
Code:
      program xttbl
      implicit none
      double precision ti, tr, co2, co3, r1, r2, r3
      double precision e1, e2, e3, es, tc1, tc2, tc3
      real*10 c1,c2,c3
      integer ni
      integer*8 jj
      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)
   
      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
Back to top
View user's profile Send private message
PaulLaidler
Site Admin


Joined: 21 Feb 2005
Posts: 5794
Location: Salford, UK

PostPosted: Thu Mar 07, 2019 1:11 pm    Post subject: Reply with quote

This "faulty" result seems to me to be a consequence of the fact that the expression "r3*tr+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".
Back to top
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Support All times are GMT + 1 Hour
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group