Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Apparently identical math expressions with different output

The following code will output different results for variables 'e' and 'f' on a x86 32 bit machine but the same results on a x86 64 bit machine. Why? Theoretically the same expression is being evaluated, but technically it is not.

#include <cstdio>
main()
{
  double a,b,c,d,e,f;
  a=-8988465674311578540726.0;
  b=+8988465674311578540726.0;
  c=1925283223.0;
  d=4294967296.0;
  e=(c/d)*(b-a)+a;
  printf("%.80f\n",e);
  f=c/d;
  f*=(b-a);
  f+=a;
  printf("%.80f\n",f);
}

Note ... 32 bit x86 code can be generated with 'gcc -m32' ,thanks @Peter Cordes https://stackoverflow.com/users/224132/peter-cordes

See also

is boost::random::uniform_real_distribution supposed to be the same across processors?

--- update for user Madivad

64 bit output 
-930037765265417043968.00000...
-930037765265417043968.00000...

32 bit output
-930037765265416519680.00000...
-930037765265417043968.00000...

The "mathematically correct" output can be given by this python code

from fractions import Fraction
a=-8988465674311578540726
b=8988465674311578540726
c=1925283223
d=4294967296
print "%.80f" % float(Fraction(c,d)*(b-a)+a)

-930037765265416519680.000...
like image 807
don bright Avatar asked Nov 19 '15 04:11

don bright


1 Answers

FLT_EVAL_METHOD.

C allows intermediate FP calculations to occur at higher/wider types depending on FLT_EVAL_METHOD. So when wider types are used and code flow differs, though mathematically equal, slightly different results may occur.


Except for assignment and cast (which remove all extra range and precision), the values yielded by operators with floating operands and values subject to the usual arithmetic conversions and of floating constants are evaluated to a format whose range and precision may be greater than required by the type. The use of evaluation formats is characterized by the implementation-defined value of FLT_EVAL_METHOD:

-1. indeterminable;
0. evaluate all operations and constants just to the range and precision of the type;
1. evaluate operations and constants of type float and double to the range and precision of the double type, evaluate long double operations and constants to the range and precision of the long double type;
2. evaluate all operations and constants to the range and precision of the long double type.
C11dr §5.2.4.2.2 9

[Edit]

@Pascal Cuoq has a useful comment on the veracity on FLT_EVAL_METHOD. In any case, FP code, optimized different along various code paths, may present different results. This may occur when FLT_EVAL_METHOD != 0 or compiler is not strictly conforming.

Concerning a detail of the post: the operation X*Y + Z done in 2 operations of * and then + could be contrasted with fma() which "compute (x × y) + z, rounded as one ternary operation: they compute the value (as if) to infinite precision and round once to the result format, according to the current rounding mode." C11 §7.12.13.1 2. Another candidate for the difference in results could be due to the application "fma" to the line e=(c/d)*(b-a)+a;

like image 185
chux - Reinstate Monica Avatar answered Nov 18 '22 02:11

chux - Reinstate Monica