Re: FWD: FLOATING-POINT CONSISTENCY, -FFLOAT-STORE, AND X86

[Toon Moene <toon at moene dot indiv dot nluug dot nl>]

Edward, Joe and Craig,

I'm not going to address all the posts separately - lets just
concentrate on Edward's:

> Also, you can use genuinely _better_ algorithms when you can rely on 
> something very close to IEEE, and that is currently pretty hard on 
> x86 with gcc.  And a touch of extended precision can really lead to 
> algorithms that give huge performance improvements (factors of 20-40
> for normal eqns. v. QR for least squares), although those examples
> are beyond the current discussion.

I'll grant you that one - I sure never researched the boundaries of IEEE
754 arithmetic.  However, what I'm challenging is that we should burden
the compiler with these considerations *by default* (I don't mind if
it's hidden behind a compile time option like -pedantic-numerics).

> And Mr. Buck's example does happen in real code.

Yes, but that doesn't make it correct, even on a strict, one size fits
all IEEE 754 machine.  The point is that the following code:

      REAL FUNCTION FINDROOT(FIRSTGUESS)
  10  FINDROOT=<expression involving FIRSTGUESS>
      IF (FINDROOT .EQ. FIRSTGUESS) RETURN
      FIRSTGUESS = FINDROOT
      GOTO 10
      END

simply is not guaranteed to work (I discussed this on comp.compilers
some months ago).  For an arbitrary choice of FIRSTGUESS and <expression
involving FIRSTGUESS> one _cannot_ prove that this won't eternally
oscillate between two numbers just one bit apart *in any precision*.

So this is the wrong way to solve such a problem.

The correct termination comparison is:

      IF (ABS(FINDROOT - FIRSTGUESS) .LT.
     ,    TOLERANCE * FIRSTGUESS) RETURN

with a suitable value of TOLERANCE (dependent on whether computations
are with 32, 64, or 80 bit REALS (Fortran 90 offers intrinsics to
parametrise this).

>  - My main concern is that there is a grid spacing that will render the
>  - basic equation of geostrophy badly approximated in 32-bit arithmetic:
> 
> And at the moment, that entirely depends on which variables happen
> to be spilled and which don't.

[ Sorry, I meant: My *only* concern is ... ]

No, because _we_ know what we're doing because we estimated the error
propagation in an independent way.

What I tried to get across is that it is not *reasonable* to punish
32-bit applications with the burden of either 1) unaligned 80 bit spills
or 2) aligned 80 bit spills that are 4 times as large as necesssary.

Yes, I know that you believe that floating point register spills are
scarce.  You probably also believe that "Real Programmers are not afraid
of 5 page long DO loops" was meant as a joke.

Cheers,

-- 
Toon Moene (toon@moene.indiv.nluug.nl)
Saturnushof 14, 3738 XG  Maartensdijk, The Netherlands
Phone: +31 346 214290; Fax: +31 346 214286
g77 Support: fortran@gnu.org; egcs: egcs-bugs@cygnus.com