pentium bug found by Kahan's number-theory methods if you wait long enough

[David Hough]

Dr. K-C Ng has coded one of Kahan's number-theoretic tests for correctly
rounded arithmetic and I have discovered that it could have discovered the
Pentium bug if allowed to run long enough (a couple of hours). 
There are also corresponding tests for multiplication and sqrt, the IEEE
operations most likely to be shortcut in sloppy hardware.
The test are designed to generate difficult cases for round-to-nearest.

The reduced test program follows.   It can be compiled with Sun's compilers
for Solaris 2 for Intel, or with suitable changes presumably on other
systems with other compilers.   The ieee_flags calls do the obvious thing.

Correct output on SPARC system:

 correct ratio     1.6458333134650994603787680105540058975268e+00 
 long double x 20709356 y 12582900 delta 0 
      double x 20709356 y 12582900 delta 7.40329e-17 
       float x 20709356 y 12582900 delta -5.96046e-08 

Correct output on 486 system:

 correct ratio     1.6458333134650994603606979915788599555526e+00 
 long double x 20709356 y 12582900 delta 0 
      double x 20709356 y 12582900 delta 7.4051e-17 
       float x 20709356 y 12582900 delta -5.96046e-08 

Incorrect output on Pentium system:

 correct ratio     1.6458333134650994603606979915788599555526e+00 
 long double x 20709356 y 12582900 delta -1.27157e-06 
      double x 20709356 y 12582900 delta -1.27157e-06 
       float x 20709356 y 12582900 delta -1.2517e-06 

Source code:

/* compile -fsingle and link with -lsunmath or -lm depending on compiler
   release */

main()
{
volatile float fx, fy, fz;
volatile double dx, dy, dz;
volatile long double lx, ly, lz, ratio;
char result[80];

ratio=1.6458333134650994603787680105540058975268l;
printf(" correct ratio %50.40Le \n",ratio);

ieee_flags("set","precision","extended",result);
lx=20709356;
ly=12582900;
lz=lx/ly - ratio;
printf(" long double x %d y %d delta %Lg \n",(int)lx,(int)ly,(long double)lz);

ieee_flags("set","precision","double",result);
dx=20709356;
dy=12582900;
dz=dx/dy - ratio;
printf("      double x %d y %d delta %g \n",(int)dx,(int)dy,(double)dz);

ieee_flags("set","precision","single",result);
fx=20709356;
fy=12582900;
fz=fx/fy - ratio;
printf("       float x %d y %d delta %g \n",(int)fx,(int)fy,(double)fz);

ieee_flags("set","precision","extended",result);
return 0;
}