5th IEEE op

[James Demmel]

There was discussion of this on the net some time ago, and it appears we are 
revisiting the same issues. Below I enclose the message I sent out some time
ago on this subject (so please delete it if you have read it already!). 

   Jim Demmel, UC Berkeley, demmelacs.berkeley.edu

---------------------------------------------------------------------------

There is a proposal going around to incorporate Kulisch/Miranker style
arithmetic as an extension of IEEE floating point arithmetic. I'd like
to summarize this proposal briefly and argue against it.

Kulisch/Miranker propose to add the ability to compute dot products of
arbitrary vectors exactly. In other words, 

   (1,1e38,-1e38).(1,1e38,1e38) = ( 1*1 + 1e38*1e38 ) - 1e38*1e38 = 1

would be computed correctly instead of as 0. This would appear to
require carrying the equivalent of at least about  2*2*38+7=159 decimal digits 
in IEEE single, and far more in double. This has been done in the past in
microcode and cache, as on the IBM 4361. They propose to add a single type,
DOTPRECISION, which one could only use to accumulate dot products, add and
subtract, and round (up, down or to nearest) to the other IEEE precisions.

Their motivation is to build algorithms using interval arithmetic with
standard precision for most operations and DOTPRECISION for a few
critical ones. Their algorithms, as implemented for example in the IBM
ACRITH package, claim to solve linear systems and get minimally wide
output intervals guaranteed to contain each solution component; i.e. the
interval bounds are adjacent floating point numbers. They have similar
codes for some eigenproblems, zeros of polynomials, linear programming,
and other problems.

A conference proceedings on this material may be found in
"A new approach to scientific computation", U. Kulisch and W. Miranker, eds.,
Academic Press, 1983.  Some of their claims are exaggerated, as described in
"Anomalies in the IBM ACRITH Package," W. Kahan and E. Leblanc, Proc.
7th IEEE Computer Arithmetic Symposium, 1985.

Here are the main arguments against adopting Kulisch's scheme:

1) None of their algorithms that I know of (I should go read their most
recent stuff to reconfirm this, though) depend for correctness on an
exact inner product. And their accuracy also depends only weakly on
how much extra precision you have: As the precision grows, you can
solve more problems to "full accuracy", but the number of these problems
is very small and gets smaller as the precision increases. You have
to work hard to concoct examples which need anywhere near the precision
required for exact inner products in order to be solved fully accurately.
So the benefit from all those extra bits is small.

2) Many of their iterative schemes (such as for linear equation solving)
depend on working precision Gaussian elimination being accurate
enough to get a convergent iteration going. As long as there is a little
bit of convergence, they compute residuals accurately to get lots of
accuracy in the answer. But if the working precision Gaussian elimination
is not good enough, their extra precision does not help them. You would
need to do the Gaussian elimination itself in higher precision, which they
cannot do with their proposed data types and instruction set.

3) They do not address exponent range limitations, which are at least
as problematic in some software development as accuracy limitations.
This is true in polynomial evaluation, for example.

4) There is really a statistical ("RISC vs CISC") question underlying this: 
Of the following reasons for failure to get a good enough answer:

   (1) not enough precision in the residuals for iterative refinement
       to work
   (2) not enough precision in the initial Gaussian elimination to get
       convergence (however slow)
   (3) exponent range limitation

which is most likely to cause failure? Kulisch only addresses (1).
Quad addresses all three, but (1) to a somewhat weaker extent
than Kulisch. I do not have statistics (I doubt anyone does) but I
feel strongly that the extra problems one can solve with quad because
it deals with (2) and (3) far outweigh the problems Kulisch can
solve but quad can't because of (1).

5) There are faster and simpler algorithms using (I believe) the same
hardware resources which they cannot use, because they limit themselves
to dot products. Again, consider polynomial evaluation: if one could
multiply a long accumulator by a single precision variable, using a
loop like:

   real s,tmp
   DOTPRECISION d1, d1
   /* compute d2=s*d1 */
   d2=0
   repeat
     tmp=d1
     d1=d1-tmp
     d2=d2+s*tmp
   until d1=0
    
then one could just use the standard algorithm to evaluate polynomials
(Horner's rule) and it would work. Instead, they need a complicated 
library routine. In other words, their proposed instruction set makes
programmers have to completely rethink things unnecessarily.

So from a software point of view, it is much easier to develop algorithms
using quad than shoehorning everything into dot products. Kahan has more
to say on this in his article cited above.

6) Trying to get tight guaranteed error bounds is a laudable goal, and
interval arithmetic is a valuable tool, but it will not be possible to do
this reliably in systems without variable precision. The reason is that if one
tries to build larger applications out of interval building blocks
(like ACRITH routines), then even if the input data is exact 
(0-width intervals), the output of the first interval routine will be 
intervals. These intervals are then fed in to a second routine. How wide 
must the output intervals of the second routine be? It is easy to show that 
in general they must be COND2 times wider than the input intervals, 
where COND2 is the condition number of the problem being solved by the 
second routine. If these intervals are then used as input to a third routine, 
their width will be magnified by another condition number factor COND3. 
If we are computing with a fixed maximum precision, eventually these 
products of COND factors will make the intermediate intervals so large 
as to be useless. The only way around this is to completely reformulate 
the problem to avoid these condition numbers, which may be very difficult 
and time-consuming (and defeats the purpose of using these libraries, which is 
getting the answer faster as well as more accurately), or else to use 
variable precision.

Thus, I do not foresee wide use of interval arithmetic until variable 
precision is adequately supported, and this cannot be done without a 
more major architectural revision than a small change in the IEEE floating
point standard.