## Why does Mathematica results differ from C++ results within machine precision?

30

4

I am trying to wrap my head around floating-point precision. Performing the following numeric operation:

275./6.*1.03692775514337


in Mathematica and C++ and taking the difference of the two results in ~0.7e-14. I expected the difference to be zero within my \$MachinePrecision of ~15.96. C++ uses double as variable type for each number. In addition C++ and Mathematica follow the IEEE 754, which should make division and multiplication exactly rounded operations.

In general I need to know why Mathematica is rounding multiplication and division differently than my C++ program, while both should yield the same result?

For anybody interested in the C++ code:

#include <iostream>

int main() {
std::cout << std::setprecision(17);
std::cout << 275./6.*1.03692775514337 - 47.52585544407112

return(0);
}


1Can you add the C++ code you used to the question and ideally the bit strings of the two resulting double precision numbers? This would make tracking down the differences (e.g. different order of operations) easier. – Thies Heidecke – 2017-12-10T20:45:16.737

1Mathematica rounds 275./6. towards zero (on my machine) even though rounding away from zero introduces less absolute error. But what's you source of "Mathematica follow the IEEE 754"? – Coolwater – 2017-12-10T20:48:56.840

@Coolwater know that you ask I am not even sure, but in its documentation (e.g. MachinePrecision in the Background & Context part) it mentions it several times so I assumed that they use it.

– Knowledge – 2017-12-11T07:34:04.287

53

Something important to keep in mind is that Mathematica parses x / y as

Times[x, Power[y, -1]]


For actual floating point division, use Divide:

Divide[275., 6.]*1.03692775514337 // InputForm

(* 47.52585544407113 *)


which should agree with the C++ result.

8Why does it parse it like that? – Rivasa – 2017-12-11T15:30:01.893

2

It was a design decision taken a long time (at least 29 years) ago... if I had to guess, I'd probably say for ease of symbolic manipulation (e.g. Times is Flat and Orderless). There are some interesting comments in an earlier Q & A, (39200).

– ilian – 2017-12-11T23:01:24.763

2

@Annabelle Probably because multiplication is performed faster than division on most hardware, and modern C/C++ compilers do this optimization as well. Another reason I can think of is because there have been known bugs in the hardware implementation of division which could render the results unreliable on some CPUs.

– undercat applauds Monica – 2017-12-12T07:29:00.340

Wow. The documentation of Divide claims to also be equivalent to x y^-1 yet it's not... and even more confusingly, try evaluating this: {Inactivate[x/y], Inactivate[Divide[x, y]]} it prints out {x * y^(-1), x / y} which makes no sense!! – user541686 – 2017-12-12T10:29:14.313

3@undercat That is a possibility, although 1989 predates the FDIV bug by some years and I am not sure what compilers did back then. I'll also note that the original behavior of Divide was to multiply by the reciprocal prior to version 3.0 (1997). – ilian – 2017-12-12T15:23:09.643

@ilian Division has always been a performance hog, see for example the latency table for Intel 8087 (circa 1980) for fmul and fdiv.

– undercat applauds Monica – 2017-12-13T02:21:10.937

16

Without code and your actual results, this question cannot be answered. Here is one thing that might help: We have a compiler that can compile to C it can show you the code it creates. So why don't you try this?

a = 275.;
b = 6.;
c = 1.03692775514337

fC = Compile[{{a, _Real}, {b, _Real}, {c, _Real}},
a/b*c,
CompilationTarget -> "C"
]


Now we can compare the two results:

fC[a, b, c] - (a/b*c)


This gives not difference on my machine. Let's look at the important part of the created C code:

<< CCodeGenerator
CCodeStringGenerate[fC, "fun"]


With this, we get the core calculation:

mreal R0_0;
mreal R0_1;
mreal R0_2;
mreal R0_3;
mreal R0_4;
R0_0 = A1;
R0_1 = A2;
R0_2 = A3;
R0_3 = 1 / R0_1;
R0_4 = R0_0 * R0_3;
R0_4 = R0_4 * R0_2;
*Res = R0_4;


So the question is, how did you calculate the result?

+1 from me, this answer gives additional context that is useful. By the way, 275./6.*1.03692775514337` is a valid C expression so I simply assumed that's what OP computed. But I understand the point of view that it could have been any other code too. – ilian – 2017-12-12T17:16:47.577

Comments are not for extended discussion; this conversation has been moved to chat.

– Kuba – 2017-12-12T20:14:43.810