From: kagel AT quasar DOT bloomberg DOT com
Date: Fri, 7 Feb 1997 09:39:26 -0500
Message-Id: <9702071439.AA27934@quasar.bloomberg.com >
To: V DOT O DOT Myskin AT inp DOT nsk DOT su
Cc: djgpp AT delorie DOT com
In-Reply-To: <6CF0313A6D@csd.inp.nsk.su> (V.O.Myskin@inp.nsk.su)
Subject: Re: double-->int: What's wrong here?
Reply-To: kagel AT dg1 DOT bloomberg DOT com

   Errors-To: postmaster AT ns1
   From: "Vyacheslav O. Myskin" <V DOT O DOT Myskin AT inp DOT nsk DOT su>
   Organization: BINP RAS, Novosibirsk, Russia
   Date: Fri, 7 Feb 1997 00:06:30 +0600
   Reply-To: V DOT O DOT Myskin AT inp DOT nsk DOT su
   Cc: djgpp AT delorie DOT com
   Priority: normal
   References: <Pine DOT SUN DOT 3 DOT 91 DOT 970206083831 DOT 4379F-100000 AT is> (message from Eli Zaretskii on Thu, 6 Feb 1997 08:44:11 +0200 (IST))
   X-Mailer: Pegasus Mail for Windows (v2.52)
   Content-Type: text
   Content-Length: 1813

   Hello!


   >    > The output is:
   >    > n1=214  n2=215  d2=215.000000
   >    >
   >    > So why  n1 is not equal to n2? It's not fun because I used n1 as an
   >    > argument to malloc() and kept getting SIGSEGVs later :(
   >
   >    Never assume floating-point computations are exact: they aren't.  The
   >    exact reason for what your example printed are immaterial (it's a long
   >    and quite dull story; you might consider adding -S to gcc command line
   >    and looking at the assembly it generates to see what's going on).  The
   >    real lesson is that you should always round the FP numbers explicitly
   >    before using them as counters of anything.  In your case, say n1 =
   >    .05/d1 + 0.5 and live happily ever after.
   >
   > But Eli, 1.0/4300.0 ==> 0.000232558139535
   > and n1 = .05/d1 ==> 215.000000000000000
   > to the limit of a double's precision (the next 5 digits would 86000 in the
   > 80bit internal x87 representation.  There is no imprecision in this
   > calculation due to binary floating point effects and adding 0.5 will not
   > change the results. Vyacheslav may have indeed found something. 
   > Vyacheslav are you using the emulator or an x87?  Which?

   I ran it on 486DX, so it's x87. The assembler output ( gcc -S ...) shows
   that n1 is stored as an integer (fistp) immediately after calculation; n2 is
   converted after saving the value of d2 and reloading it (as it should be for
   unoptimized output). The ONLY difference is this fstp/fld combination. Is it
   OK that the data saved in memory is not identical to that in the FPU
   registers (can't assume anything else) ? Shame on me, my competence in FPU
   topics is equal to zero, hope haven't misspelled instruction names ;)

OK, it seems that my investigation of Vyacheslav's problem via comparison of
IEEE -vs- 'REAL' arithmetic did not go far enough.  The value of the first
division is the problem:

1.0/4300.0		64-bit IEEE =	0.000232558139535
			80-bit IEEE =	0.00023255813953488370
			'REAL' value=	0.000232558139534883720

This difference in the last representable digit affects the result:

0.05/0.000232558139535		= 214.999999999892500
0.05/0.00023255813953488370	= 215.00000000000001935000
0.05/0.000232558139534883720	= 215.00000000000000086000

My error yesterday was that I used a truncated calculation of the first
division, to 15 digits; for the 64-bit IEEE calculation of the final result the
last digit was a 4 rather than a 5 and this 'corrected' the result.  Eli, my
apologies, this was indeed a rounding/truncation problem caused by saving the
intermediate result to a double.

-- Art S. Kagel, kagel AT quasar DOT bloomberg DOT com

A proverb is no proverb to you 'till life has illustrated it.  -- John Keats