
forums.silverfrost.com Welcome to the Silverfrost forums

View previous topic :: View next topic 
Author 
Message 
mecej4
Joined: 31 Oct 2006 Posts: 1215

Posted: Wed Feb 27, 2019 3:33 pm Post subject: Error in NINT(), 32bit FTN95 


FTN95 (32bit, 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 IEEE64 format. I give below a cutdown 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 IEEE64 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 32bit 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*2d2)
tr=ni*50d0+25d0co3
e1=r3*tr
e2=r2*(co3co2)
e3=r1*co2
es=e1+e2+e3+1.0d0
tc1=nint(r3*tr+r2*(co3co2)+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 


PaulLaidler Site Admin
Joined: 21 Feb 2005 Posts: 6092 Location: Salford, UK

Posted: Thu Feb 28, 2019 10:11 am Post subject: 


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

Back to top 


JohnCampbell
Joined: 16 Feb 2006 Posts: 2145 Location: Sydney

Posted: Thu Feb 28, 2019 11:34 am Post subject: 


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*2d2)
tr=ni*50d0+25d0co3
e1=r3*tr
e2=r2*(co3co2)
e3=r1*co2
es=e1+e2+e3+1.0d0
tc1=nint(r3*tr+r2*(co3co2)+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*(co3co2)+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 32bit are: Code:  [FTN95/Win32 Ver. 8.40.0 Copyright (c) Silverfrost Ltd 19932018]
[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 roundoff 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 roundoff 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 


mecej4
Joined: 31 Oct 2006 Posts: 1215

Posted: Thu Feb 28, 2019 12:53 pm Post subject: 


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, FTN9532 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 (32bit) 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 32bit 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/IEEE754/ . 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 


LitusSaxonicum
Joined: 23 Aug 2005 Posts: 2055 Location: Yateley, Hants, UK

Posted: Fri Mar 01, 2019 1:40 pm Post subject: 


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 


mecej4
Joined: 31 Oct 2006 Posts: 1215

Posted: Sat Mar 02, 2019 3:00 am Post subject: 


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 


JohnCampbell
Joined: 16 Feb 2006 Posts: 2145 Location: Sydney

Posted: Sat Mar 02, 2019 6:06 am Post subject: 


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*2d2)
tr=ni*50d0+25d0co3
e1=r3*tr
e2=r2*(co3co2)
e3=r1*co2
es=e1+e2+e3+1.0d0
tc1=nint(r3*tr+r2*(co3co2)+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*(co3co2)+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*(co3co2)+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*(co3co2)+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 


PaulLaidler Site Admin
Joined: 21 Feb 2005 Posts: 6092 Location: Salford, UK

Posted: Thu Mar 07, 2019 1:11 pm Post subject: 


This "faulty" result seems to me to be a consequence of the fact that the expression "r3*tr+r2*(co3co2)+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 




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
