Re: strange double comparison results with -O[12] on x86(-32)

[Vincent Lefevre <vincent+gcc at vinc17 dot org>]

On 2005-06-13 19:41:23 -0400, Robert Dewar wrote:
> Dave Korn wrote:
> > >  ... or it's a bug in the libc/crt-startup, which is where the
> >hardware rounding mode is (or should be) set up ...
> 
> Well if you think that the operations should reflect IEEE 64-bit
> semantics (which is the only rationale for mucking with the rounding
> mode in the startup), then this is only an approximate fix, since
> the intermediate values still have excessive range, so you don't get
> infinities when expected.

Yes, and with optimizations, this makes things difficult to control.
And one may also get infinities when they are not expected. With the
program below, one gets the following results on x86:

With value = 1:
[1a] x = 3.595e+308   y = 3.595e+308   z = -3.595e+308

With value = 2:
[1a] x = 3.595e+308   y = 3.595e+308   z = -3.595e+308
[1b] x = 3.595e+308   y = 3.595e+308   z = -3.595e+308

With value = 3 and CR = 0:
[2a] x = 3.595e+308   y = 3.595e+308   z = -3.595e+308

With value = 3 and CR = 1:
[2a] x = 1.798e+308   y = 1.798e+308   z = -inf

With value = 4 and CR = 0:
[1a] x = 3.595e+308   y = 3.595e+308   z = -3.595e+308
[2a] x = 3.595e+308   y = 1.798e+308   z = -inf

With value = 4 and CR = 1:
[1a] x = 3.595e+308   y = 3.595e+308   z = -1.798e+308
[2a] x = 3.595e+308   y = 1.798e+308   z = -1.798e+308

With value = 5:
[1a] x = 3.595e+308   y = 1.798e+308   z = -1.798e+308
[1b] x = 3.595e+308   y = 1.798e+308   z = -1.798e+308
[2a] x = 1.798e+308   y = 1.798e+308   z = -1.798e+308

In particular, it should have been impossible to get infinities
with this example.


#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <fenv.h>

#pragma STDC FENV_ACCESS ON

#define P1 ((MODE) % 3)
#define P2 ((MODE) / 3)

int main (int argc, char *argv[])
{
  double x, y, z;

  if (argc != 2)
    {
      fprintf (stderr, "Usage: overflow <double>\n");
      exit (1);
    }

  if (fesetround (FE_DOWNWARD))
    {
      fprintf (stderr, "Can't set rounding mode to FE_DOWNWARD\n");
      exit (1);
    }

  x = atof (argv[1]);
  x *= DBL_MAX;
  y = x + 0.0;
  z = x * -1.0;
#if P1
  printf ("[1a] x = %-12.4Lg y = %-12.4Lg z = %-12.4Lg\n",
          (long double) x, (long double) y, (long double) z);
#if P1 == 2
  printf ("[1b] x = %-12.4Lg y = %-12.4Lg z = %-12.4Lg\n",
          (long double) x, (long double) y, (long double) z);
#endif
#endif
#if CR
  printf ("\r");
#endif
  x += 0.0;
  y += 0.0;
  z += 0.0;
#if P2
  printf ("[2a] x = %-12.4Lg y = %-12.4Lg z = %-12.4Lg\n",
          (long double) x, (long double) y, (long double) z);
#if P2 == 2
  printf ("[2b] x = %-12.4Lg y = %-12.4Lg z = %-12.4Lg\n",
          (long double) x, (long double) y, (long double) z);
#endif
#endif
  return 0;
}


I wonder if one can make such a program crash with a FPE signal by
using something like:

  #define _GNU_SOURCE
  #include <fenv.h>
[...]
  feclearexcept (FE_OVERFLOW);
  feenableexcept (FE_OVERFLOW);
  x += 0.0;

I didn't manage.

-- 
Vincent Lefèvre <vincent@vinc17.org> - Web: <http://www.vinc17.org/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.org/blog/>
Work: CR INRIA - computer arithmetic / SPACES project at LORIA