Mail Archives: djgpp/1997/02/07/10:19:07
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
- Raw text -