Re: R-alpha: R 0.12 alpha: problem with qt()

Douglas Bates (bates@stat.wisc.edu)
Wed, 6 Nov 96 14:21 CST


Message-Id: <m0vLET7-000hhwC@franz.stat.wisc.edu>
Date: Wed, 6 Nov 96 14:21 CST
From: Douglas Bates <bates@stat.wisc.edu>
To: Peter Dalgaard BSA <pd@kubism.ku.dk>
Subject: Re: R-alpha: R 0.12 alpha: problem with qt()
In-Reply-To: <x2ohhbdn8q.fsf@bush.kubism.ku.dk>
 <x2pw1rdo9i.fsf@bush.kubism.ku.dk>

>>>>> "Peter" == Peter Dalgaard BSA <pd@kubism.ku.dk> writes:

  Peter> Just following up with a little tidbit: Extending the test
  Peter> printout to fprintf(stderr, "adj=%30.16g\n", adj);
  Peter> fprintf(stderr, "tx=%30.16g\n", tx); fprintf(stderr,
  Peter> "prev=%30.16g\n", prev); fprintf(stderr, "y=%30.16g\n", y);
  Peter> fprintf(stderr, "sq=%30.16g\n", sq); fprintf(stderr,
  Peter> "g=%30.16g\n", g); ... makes the endless loop disappear!

OK, so the problem is caused by a comparison that is too restrictive.

One difference between optimized and unoptimized floating point code
is whether values in floating point (FP) registers are stored into
memory and reloaded between operations.  

The floating point unit on Intel processors has 80 bit registers with
an expanded mantissa.  A double precision value is stored in memory as
a 64 bit IEEE floating point number.  Thus storing a value from an FP
register into a memory location is actually a rounding operation as
well as a storage operation.

If you do a comparison of two "double" values from memory, they will
be converted from 64 bits to 80 bits when they are loaded into FP
registers then compared.  (In fact, for equality they may actually be
compared as bit patterns.)  Of course they will be equal if they matched
as 64 bit quantities.  If you compare two values stored in registers
directly after their calculation they will be equal only if they match
as 80 bit quantities.

The reason that your debugging print statements changed the behaviour
is because they caused different usage of the floating point registers
in the critical section of the calculation.

There is actually a compiler option to get around this problem

 `-ffloat-store'
      Do not store floating point variables in registers, and inhibit
      other options that might change whether a floating point value is
      taken from a register or memory.

      This option prevents undesirable excess precision on machines such
      as the 68000 where the floating registers (of the 68881) keep more
      precision than a `double' is supposed to have.  For most programs,
      the excess precision does only good, but a few programs rely on the
      precise definition of IEEE floating point.  Use `-ffloat-store' for
      such programs.

The better solution, of course, is to examine the code and change the
nature of the comparisons being made.
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
r-testers mailing list -- To (un)subscribe, send
subscribe	or	unsubscribe
(in the "body", not the subject !)  To: r-testers-request@stat.math.ethz.ch
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-