I am getting different floating-point rounding under different build/execute scenarios. Notice the `2498`

in the second run below...

#include <iostream> #include <cassert> #include <typeinfo> using std::cerr; void domath( int n, double c, double & q1, double & q2 ) { q1=n*c; q2=int(n*c); } int main() { int n=2550; double c=0.98, q1, q2; domath( n, c, q1, q2 ); cerr<<"sizeof(int)="<<sizeof(int)<<", sizeof(double)="<<sizeof(double)<<", sizeof(n*c)="<<sizeof(n*c)<<"\n"; cerr<<"n="<<n<<", int(q1)="<<int(q1)<<", int(q2)="<<int(q2)<<"\n"; assert( typeid(q1) == typeid(n*c) ); }

Running as a 64-bit executable...

$ g++ -m64 -Wall rounding_test.cpp -o rounding_test && ./rounding_test sizeof(int)=4, sizeof(double)=8, sizeof(n*c)=8 n=2550, int(q1)=2499, int(q2)=2499

Running as a 32-bit executable...

$ g++ -m32 -Wall rounding_test.cpp -o rounding_test && ./rounding_test sizeof(int)=4, sizeof(double)=8, sizeof(n*c)=8 n=2550, int(q1)=2499, int(q2)=2498

Running as a 32-bit executable under valgrind...

$ g++ -m32 -Wall rounding_test.cpp -o rounding_test && valgrind --quiet ./rounding_test sizeof(int)=4, sizeof(double)=8, sizeof(n*c)=8 n=2550, int(q1)=2499, int(q2)=2499

Why am I seeing different results when compiling with `-m32`

, and why are the results different again when running valgrind?

My system is Ubuntu 14.04.1 LTS x86_64, and my gcc is version 4.8.2.

EDIT:

In response to the request for disassembly, I have refactored the code a bit so that I could isolate the relevant portion. The approach taken between `-m64`

and `-m32`

is clearly much different. I'm not too concerned about why these give a different rounding result since I can fix that by applying the `round()`

function. The most interesting question is: why does valgrind change the result?

rounding_test: file format elf64-x86-64 < 000000000040090d <_Z6domathidRdS_>: < 40090d: 55 push %rbp < 40090e: 48 89 e5 mov %rsp,%rbp < 400911: 89 7d fc mov %edi,-0x4(%rbp < 400914: f2 0f 11 45 f0 movsd %xmm0,-0x10(%r < 400919: 48 89 75 e8 mov %rsi,-0x18(%rb < 40091d: 48 89 55 e0 mov %rdx,-0x20(%rb < 400921: f2 0f 2a 45 fc cvtsi2sdl -0x4(%rbp), < 400926: f2 0f 59 45 f0 mulsd -0x10(%rbp),%x < 40092b: 48 8b 45 e8 mov -0x18(%rbp),%r < 40092f: f2 0f 11 00 movsd %xmm0,(%rax) < 400933: f2 0f 2a 45 fc cvtsi2sdl -0x4(%rbp), < 400938: f2 0f 59 45 f0 mulsd -0x10(%rbp),%x < 40093d: f2 0f 2c c0 cvttsd2si %xmm0,%eax < 400941: f2 0f 2a c0 cvtsi2sd %eax,%xmm0 < 400945: 48 8b 45 e0 mov -0x20(%rbp),%r < 400949: f2 0f 11 00 movsd %xmm0,(%rax) < 40094d: 5d pop %rbp < 40094e: c3 retq < | rounding_test: file format elf32-i386 > 0804871d <_Z6domathidRdS_>: > 804871d: 55 push %ebp > 804871e: 89 e5 mov %esp,%ebp > 8048720: 83 ec 10 sub $0x10,%esp > 8048723: 8b 45 0c mov 0xc(%ebp),%eax > 8048726: 89 45 f8 mov %eax,-0x8(%ebp > 8048729: 8b 45 10 mov 0x10(%ebp),%ea > 804872c: 89 45 fc mov %eax,-0x4(%ebp > 804872f: db 45 08 fildl 0x8(%ebp) > 8048732: dc 4d f8 fmull -0x8(%ebp) > 8048735: 8b 45 14 mov 0x14(%ebp),%ea > 8048738: dd 18 fstpl (%eax) > 804873a: db 45 08 fildl 0x8(%ebp) > 804873d: dc 4d f8 fmull -0x8(%ebp) > 8048740: d9 7d f6 fnstcw -0xa(%ebp) > 8048743: 0f b7 45 f6 movzwl -0xa(%ebp),%ea > 8048747: b4 0c mov $0xc,%ah > 8048749: 66 89 45 f4 mov %ax,-0xc(%ebp) > 804874d: d9 6d f4 fldcw -0xc(%ebp) > 8048750: db 5d f0 fistpl -0x10(%ebp) > 8048753: d9 6d f6 fldcw -0xa(%ebp) > 8048756: 8b 45 f0 mov -0x10(%ebp),%e > 8048759: 89 45 f0 mov %eax,-0x10(%eb > 804875c: db 45 f0 fildl -0x10(%ebp) > 804875f: 8b 45 18 mov 0x18(%ebp),%ea > 8048762: dd 18 fstpl (%eax) > 8048764: c9 leave > 8048765: c3 ret

Edit:It would seem that, at least a long time back, valgrind's floating point calculations wheren't quite as accurate as the "real" calculations. In other words, this MAY explain why you get different results. See this question and answer on the valgrind mailing list.

Edit2:And the current valgrind.org documentation has it in it's "core limitations" section here - so I would expect that it is indeed "still valid". In other words the documentation for valgrind says to expect differences between valgrind and x87 FPU calculations. "You have been warned!" (And as we can see, using sse instructions to do the same math produces the same result as valgrind, confirming that it's a "rounding from 80 bits to 64 bits" difference)

Floating point calculations WILL differ slightly depending on exactly how the calculation is performed. I'm not sure exactly what you want to have an answer to, so here's a long rambling "answer of a sort".

Valgrind DOES indeed change the exact behaviour of your program in various ways (it emulates certain instructions, rather than actually executing the real instructions - which may include saving the intermediate results of calculations). Also, floating point calculations are well known to "not be precise" - it's just a matter of luck/bad luck if the calculation comes out precise or not. 0.98 is one of many, many numbers that can't be described precisely in floating point format [at least not the common IEEE formats].

By adding:

cerr<<"c="<<std::setprecision(30)<<c <<"\n";

we see that the output is `c=0.979999999999999982236431605997`

(yes, the actual value is 0.979999...99982 or some such, remaining digits is just the residual value, since it's not an "even" binary number, there's always going to be something left.

This is the `n = 2550;`

, `c = 0.98`

and `q = n * c`

part of the code as generated by gcc:

movl $2550, -28(%ebp) ; n fldl .LC0 fstpl -40(%ebp) ; c fildl -28(%ebp) fmull -40(%ebp) fstpl -48(%ebp) ; q - note that this is stored as a rouned 64-bit value.

This is the `int(q)`

and `int(n*c)`

part of the code:

fildl -28(%ebp) ; n fmull -40(%ebp) ; c fnstcw -58(%ebp) ; Save control word movzwl -58(%ebp), %eax movb $12, %ah movw %ax, -60(%ebp) ; Save float control word. fldcw -60(%ebp) fistpl -64(%ebp) ; Store as integer (directly from 80-bit result) fldcw -58(%ebp) ; restore float control word. movl -64(%ebp), %ebx ; result of int(n * c) fldl -48(%ebp) ; q fldcw -60(%ebp) ; Load float control word as saved above. fistpl -64(%ebp) ; Store as integer. fldcw -58(%ebp) ; Restore control word. movl -64(%ebp), %esi ; result of int(q)

Now, if the intermediate result is stored (and thus rounded) from the internal 80-bit precision in the middle of one of those calculations, the result may be subtly different from the result if the calculation happens without saving intermediate values.

I get identical results from both g++ 4.9.2 and clang++ -mno-sse - but if I enable sse in the clang case, it gives the same result as 64-bit build. Using `gcc -msse2 -m32`

gives the 2499 answer everywhere. This indicates that the error is caused by "storing intermediate results" in some way or another.

Likewise, optimising in gcc to -O1 will give the 2499 in all places - but this is a coincidence, not a result of some "clever thinking". If you want correctly rounded integer values of your calculations, you're much better off rounding yourself, because sooner or later `int(someDoubleValue)`

will come up "one short".

Edit3:And finally, using `g++ -mno-sse -m64`

will also produce the same `2498`

answer, where using `valgrind`

on the same binary produces the `2499`

answer.