Re: Relying on precise integer calculation with double

[Andrew Haley <aph at redhat dot com>]

Daniel Bratell writes:

 > This is likely to be the one FAQ you've all learned to answer but
 > as I was bitten by it, I just wanted to make sure that what I saw
 > was expected.
 > 
 > I've used -ffast-math for a slight speedup of floating point
 > arithmetics, but I've also used doubles to store integers. Since
 > integers (up to a certain size) can be stored without loss of
 > precision in a double it's been no problems until I started using
 > gcc 4.0 (this is with gcc 4.0.1 on a fedora core 3, x86 system).
 > 
 > Now I see my "exact" integer arithmetic converted to inexact
 > fractional double arithmetic. More specifically a division by 7
 > (can be exactly represented in a double) is converted to a
 > multiplication by 0.14285...  (about 1/7) and then rounding errors
 > can make my "integers" get another value.

Yes.  That's one of the things -ffast-math does.

 > Is this to be expected? Should I stop using
 > -funsafe-math-optimizations? Somehow I've been tricked into
 > thinking that it would be safe to use that flag unless I relied on
 > extreme precision or behaviour related to NaN or +-Inf, and that
 > normal simple arithmetics should still give the same result.
 > 
 > So there you have it. Now you can kick me. ;-)
 > 
 > Test program below. Expected 1, get 3. Compiled with
 > gcc -O2 -ffast-math
 > 
 > ---------------------------------
 > 
 > #include <stdio.h>
 > 
 > int foo(double diff_days)
 > {
 >    if(diff_days/7.0 != (int)(diff_days/7.0))
 >      return 3;
 >    return 1;
 > }
 > 
 > int main(int argc, char** argv)
 > {
 >    printf("Return value (7.0/7.0) = %d\n", foo(7.0));
 > }
 > 
 > --------------------------------

This is due to the most popular gcc bug of all, PR323.  See
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323.


gcc implements division by 7 by multiplying by
0.14285714285714284921269268124888185.

The result of this is 0.99999999999999994448884876874217298 (raw 0x3ffefffffffffffffc00).

If you stored this in a variable and used -ffloat-store, the problem
you see wouldn't happen.

  //   if(diff_days/7.0 != (int)(diff_days/7.0))
   if(x != (int)(diff_days/7.0))

 $ gcc -O2 -ffast-math p.c -save-temps -g -ffloat-store
 $ ./a.out 
                             1

Andrew.