Floating point standards

[James Demmel]

This is a response to comments from Ed Barkmeyer of NIST on floating point
models and programming. I share his displeasure with the lack of standards
with which to write portable floating point code, and would like to briefly
explain why things are not simple enough for standards like IEEE 754 (a good
one) or the LIA (ISO 10967, aka LCAS - a bad one) to solve all problems.

There are three styles of reasoning one uses for floating point programs:

  Algebraic, where we pretend floats are real numbers, with all the properties
             of reals,
  Analytic,  which is just float(a op b) = (a op b)(1+ eps), where we still
             treat everything as real numbers, and have a bound on abs(eps).
             We may also have over/underflow thresholds, outside of which
             results are 0, undefined or exceptional.
  Discrete,  where we use the fact that floats are bit strings with certain
             properties.

These are very different kinds of arguments. Each has its range of use,
and does not work (or is not useful) if used inappropriately. It is the
art of the numerical analyst to know which to use when, and some correctness
arguments may well use all three. Using a "high" level argument (algebra
being highest and discrete lowest) makes the arguments simplest and most
compact, but may ignore possibly significant details, and using a "low" level
argument can lead to too many details to deal with.  Any attempt to 
standardize or formalize this art in such a way as to exclude one of these 
styles of argument will not help people write better software.

Unfortunately, the discrete model varies significantly from machine to machine,
and compiler to compiler, and will do so for the foreseeable future, even on
machines claiming conformance to IEEE 754. I do not think that we can or
perhaps even should try to stop all these innovations. Instead, we should make
sure the details are well documented and available to the programmer at
run-time and/or compiler-time through environmental enquiries. The LIA is an 
attempt to make a model half way between analytic and discrete, with some 
flavor of both.  There are codes that are portable and correct but cannot be 
proven using LIA, because it does not permit certain discrete arguments to be 
used (see the article "Analysis and Refutation of the LCAS" by W. Kahan,
published in the SIGNUM Newsletter and ACM SIGPLAN 25,1 (Jan 1990)).
In particular, there is a code intended for inclusion in LAPACK that works
on all machines except Cray, but cannot be proven using LIA. In fact, it
can be shown to fail on a hypothetical, nonexistent machine conforming to 
LIA, unless modified to slow it down significantly. (Cray has announced they
will conform to IEEE 754 in the future, and so this problem will go away.)

So if the LCAS is not a solution, what about IEEE 754? Unfortunately, it does
not specify how expressions are evaluated if extra precision is available
(IEEE extended registers on Intel microprocessors, and fused multiply-add
on the RS/6000), nor does it say how exception handling appears to the user,
or how much it costs. In some implementations, IEEE infinity arithmetic is
as fast as standard floating point. On the DEC alpha, it appears to require
a trap, and so is at least 30 times slower than unexceptional floating point.
So a programmer writing a code intended to be portable across IEEE 
implementations and be fast might unfortunately decide to be as paranoid about
avoiding overflow as in the past, since the 30x slowdown might be too
high a price to pay, even if rare.

So what should we do? For expert programmers, I would like very detailed
environmental enquiries available at run-time and/or compile-time, so expert
codes can automatically configure themselves to the local system. This means
how each operation rounds, overflows, has exceptions handled, etc.
This also include complex arithmetic, since some compilers lose half the
exponent range on complex division by implementing it so as to require
squaring entries of the denominator.  It would also be very useful to have 
directives that direct the compiler to compile the next few lines of code 
without optimization, rearrangements, changes of precision, etc., so that we 
can straightforwardly use discrete-style arguments for critical sections of 
code.

This all will make experts happy, but not your average programmer uninterested
in all the detail. For these programmers (the vast majority) we should also
have simpler enquiries which return a single machine epsilon, 
overflow threshold, etc. adequate (if conservative) for every operation.

Since machines and compilers are changing, it is difficult to design detailed
enquiries that will cover all future innovations, perhaps even on IEEE 
machines. Designing these enquiries is best left to the people who will use
them, numerical analysts, but working in conjunction with language and compiler
people. Compiler and language people, who have traditionally made these
decisions, are quite understandably not familiar with the all ramifications of
their decisions.

So the burden is on us numerical analysts to design such a detailed standard.
If we don't supply it, we will spend a lot of time explaining why other
proposals like LIA are inadequate.

  Jim Demmel
  CS Division and Math Dept.
  UC Berkeley