floating point arithmetics on ix86

[Matthias Klose <doko at cs dot tu-berlin dot de>]

I am looking for some pointers / FAQ's how to use fp arithmetics on
the ix86 architecture. The background are "bug" reports as found on
http://www.debian.org/Bugs/db/36/36876.html and
http://www.debian.org/Bugs/db/40/40189.html (or appended here).
AFAIK these are issues with the machine precision, which could be
described in a FAQ added to the binary packages.

First Report:
=============

The following program inconsistently chooses minimal
array elements in case of equal elements being present.

Here is the code (ttr.c):
------------------------
#include <stdio.h>
#include <math.h>


int
main() {

  double
a[4]={0.006776,2.5,3.0,0.006776},w[4]={0.000030,3.7,3.0,0.000030},t,amin;
  int ir,i,nr=4;
  double tol=9.31322574615479e-10*3.0;

  for (i=0,ir=-1,amin=DBL_MAX;i<nr;i++) 
    if (a[i]>tol) {
      t=w[i]/a[i];
      if (t<amin || ir==-1) {
        amin=t;
        ir=i;
      }
    }
  
  printf("ir=%d\n",ir);

  for (i=0,ir=-1,amin=DBL_MAX;i<nr;i++) 
    if (a[i]>tol && ((t=w[i]/a[i])<amin || ir==-1)) {
        amin=t;
        ir=i;
      }
  
  printf("ir=%d\n",ir);

  exit(0);

}
---------------  
and here is the execution:
---------------
intech19:~$ egcc -g -Wall ttr.c -o ttr -lm
intech19:~$ ./ttr
ir=0
ir=3
-----------------

Second report:
==============

The following program gives erratic/correct answers under different
conditions. The program is self explanatory.

The expected behaviour is to get the value of the variable 'q' to
equal the value of the variable 'shift' . Erratic behaviour refers to
deviation from this equality. The deviation can be controlled by
increasing/decreasing the value of 'factor' . The variable take
GENERIC values.

The expected behaviour is seen on SGI and SPARC machines with their
native compilers and also with gcc (2.7x), with or without
optimization.

On pentium with gcc/egcs the following happens:

Without optimization:   only those x which are powers of 2 give q = shift.
With optimization:      all give q = shift.

The deviations are NOT an artifact of printf formats. If subsequents
iterations are done with further code, these deviations actually get
amplified, and thus are genuine. In a non trivial code, this
uncontrolled error accumulation is a disaster.

While it is true that beyond 15 decimal places, one has garbage bits
in the internal representation, a suitable round-off/truncation etc
should discard these bits correctly and consistantly for generic
values. (Other platforms mentioned above seem to take care of this so
also optimization on pentium platform.)

In another context, with a somewhat different program, the expected
behaviour resulted WITH-OUT optimization but WITH optimization, the
behaviour was erratic.

/*---------------------------------------------------------------------*/
# include <stdio.h>

main(argc, argv)
int argc;
char **argv ;

{
        double q, q0, x, factor, shift;

        factor = 100000000.000001L ;    /* Choose any typical fraction */
        shift = 0.0L ;

        for(x = 0.0L; x < 257.0L; x = x + 1.0L )        {
                q0 = (factor)*x ;
                q =  q0 - factor*x + shift ;
                if( q == shift )
                        printf(" TRUE:\t%20.8E%20.8E%20.8E\n", x, q0, q);
                else            
                        printf("FALSE:\t%20.8E%20.8E%20.8E\n", x, q0, q);
        }
}

/*---------------------------------------------------------------------*/