Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!uunet!munnari.oz.au!bruce!monu0.cc.monash.edu.au!monu6!minyos.xx.rmit.oz.au!goanna!ok
From: ok@goanna.cs.rmit.oz.au (Richard A. O'Keefe)
Newsgroups: comp.arch
Subject: Re: IEEE arithmetic (Goldberg paper)
Message-ID: <6260@goanna.cs.rmit.oz.au>
Date: 13 Jun 91 12:18:41 GMT
References: <9106070109.AA02137@ucbvax.Berkeley.EDU> <1991Jun10.234313.2157@ingres.Ingres.COM>
Organization: Comp Sci, RMIT, Melbourne, Australia
Lines: 84

In article <1991Jun10.234313.2157@ingres.Ingres.COM>, jpk@ingres.com (Jon Krueger) writes:
> jbs@WATSON.IBM.COM (James B. Shearer) writes:

> > One of the problems with IEEE's infs and
> > NaNs is that they do not fit well with fortran.

> Is this a problem with IEEE floats or with FORTRAN?

With all this noise about IEEE 754/854 and Fortran and 0.0*X and so on,
I thought it was time to hear a voice from the past.

    We avoided algebraic constructions leading to known numerical
    problems.  For example, when floating-point arithmetic is truncated
    or rounding occurs before renormalization, the expression 2.0*x
    often contains error.  We used the expression x+x instead, because
    it is always error-free on binary machines and never worse than
    2.0*x on other machines.  Similarly, the expression 1.0-x was
    replaced, on machines that lacked guard digits, by the expression
    (0.5-x)+0.5, to protect precision for 0.0 < x < 1.0.
	...	
    Perhaps the most elegant example of error pinpointing occurred in
    initial tests at Wisconsin for the complete elliptic integral K(m).
    For m close to 1.0, function accuracy depends critically upon
    accurate evaluation of the expression (1.0 - m).  As described
    earlier, this expression was replaced by the more stable expression
    (0.5 - m) + 0.5.  Despite this tactic, the tests detected large
    errors for m close to 1.0.  ...  Univac\s optimizing compiler had
    ``helpfully'' restored the earlier expression  ... [we] designed a
    simple counter-measure.


	FUNPACK -- a Package of Special Function Routines.
	W. J. Cody.

MORAL:  _any_ compiler which applies "algebraic" identities that apply
to mathematical reals but not to the target machine's floats is going to
make it rather harder for people to write accurate

In another chapter of the same book (Sources and Development of
Mathematical Software), Cody observes

	High-quality software is supposed to be fail-safe and transportable;
	it must work properly regardless of quirks in the host arithmetic
	system.  Software production is seriously hampered when computer
	arithmetic violates simple arithmetic properties such as
		1.0 * X = X
		X * Y = Y * X
	and	X + X = 2.0 * X
	There exist mainframes of recent [1980] design in which each of
	these properties fails for appropriate floating-point X and Y.
	[He's talking about *finite* values, folks.  Not inf, nan, or eps.]
	Worse yet, on some machines there exist floating-point X > 0.0 such
	that
		1.0 * X = 0.0
		X + X = 0.0
	or	SQRT(X)**2 = overflow or underflow.
	ALL THESE ANOMALIES ARE TRACEABLE TO ENGINEERING ECONOMIES.
	Computer architects repeatedly ignore complaints about such
	mathematical atrocities, and new anomalies seem to appear
	with each new machine.
	...
	Not only is [IEEE 754] fee of anomalies, but it also contains
	new features SPECIFICALLY REQUESTED AND DESIGNED BY NUMERICAL
	ANALYSTS WITH SOFTWARE EXPERIENCE.

Elsewhere in the same book (a collection of papers edited by W.R.Cowell
and published in 1984) someone else complains that on the Cray-1 there
is a _finite_ value X such that 1.0*X overflows.  What price the
optimisation 1.0*X => X on that machine?

I suggest that if anyone wants to partially evaluate something like
	sum(i*c(i+1), i = 0,n)
and have it turn first into 0*c(1) + 1*c(2) + 2*c(3) and then into
c(2) + 2*c(3), the appropriate thing is for the partial evaluator to do
that, not the compiler, and to do it only when licensed to by an
explicit declaration "USE-REAL-AXIOMS(c(*))" (or when it knows that the
target machine obeys the laws 0*x = 0, 1*x = x, 0+x = x). It took me a
couple of years to get over the "optimisation" I was taught to use of
replacing X/Y by X*YINV, with YINV precalculated; the last thing I want
is an ``optimising'' compiler putting the less accurate expression back.

-- 
Q:  What should I know about quicksort?   A:  That it is *slow*.
Q:  When should I use it?  A:  When you have only 256 words of main storage.