fused multiply add changes exceptions too

[David Hough]

Bob Alverson pointed out something that I had overlooked previously:
fused multiply-add will affect the results by more than a bit when 
intermediate exceptions disappear that would otherwise intervene.
An example:

#include <math.h>

double addmul(x,y,z)
double x,y,z;
{
return y*z-x;
}

main()
{
double huge, big;
huge=exp(1000.0);
big=1.0e200;
printf(" huge %g \n",huge);
printf(" big %g \n",big);
printf(" addmul %g \n",addmul(huge,big,big));
}

On an RS/6000:

 cc -O addmul.c -lm -qnomaf ; a.out
 huge INF 
 big 1e+200 
 addmul NaNQ 


This is a conventional IEEE result.  But

 cc -O addmul.c -lm ; a.out
 huge INF 
 big 1e+200 
 addmul -INF 

exploits the "infinite precision" intermediate result to avoid an 
exception.  This is different, but better.   Most other fused multiply-add
implementations aren't so careful and are neither more accurate nor less
prone to exceptions, just faster while
gratuitously different.

Of course this doesn't help in a dot product accumulation, of course;
to avoid gratuitous exceptions in the final result, the entire accumulation
has to occur with a wider exponent; that suggests operators like Fortran
DPROD which are defined in paper hardware at least for SPARC V8.