Optimizing a+b-a-b etc.

[Thomas M. Breuel]

|These sorts of expressions appear in code which simulates double precision,
|and are an essential part of least one current LAPACK routine, and I plan
|to use it more widely in an IEEE arithmetic version of LAPACK currently
|under development. A typical code sequence is the following, which returns
|A and B representing the exact sum of C and D, where B is precisely the
|roundoff committed while computing A.
|
|  A = C+D
|  if (abs(C)<abs(D)) swap C and D
|  B = (A-C)-D
|
|Parentheses must be respected. It would break the code to replace B by 0.
|For more details, see, for example, the paper by Doug Priest in the 
|Proceedings of the 10th Symposium on Computer Arithmetic (IEEE, 1991). Jim Demmel

There may indeed be a need for such manipulations of floating point
numbers. However, I believe such unoptimizable floating point operations 
should not be the default. For most algorithms, I want the compiler
to optimize expressions to run faster, as long as the optimized expression
gives results that are no less accurate (but possibly more accurate)
than the unoptimized expression; this means that it should be OK
to optimize "(((a+b)-a)-b)" to "0".

For those (few) cases in which you depend on roundoff errors, another
set of functions could be provided that are guaranteed by the compiler
not to be optimized away (in fact, you can achieve this effect using
separate compilation in many C systems). Using such functions,
you could then write something like "fsub(fsub(fadd(a,b),a),b)"
As I mentioned before, an alternative would be to control optimization
of floating point expression via compiler directives. In either
case, a standard would be useful.

But, please don't spoil C for the "rest of us", who don't want to have
to hand-optimize every floating point expression just because some
esoteric uses of floating point have ursurped the standard 
notation for floating point.

					Thomas.

PS: I believe that code like the one you gave above is already
unportable in C. If you are going to make extensive use of such
operations, you may want to use explicit macros or functions for the
individual floating point operations. This will make porting it to
compilers that have a different notion of "optimization" than you
do much easier.