29
$\begingroup$

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);
}
$\endgroup$
  • 1
    $\begingroup$ Can 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. $\endgroup$ – Thies Heidecke Dec 10 '17 at 20:45
  • 1
    $\begingroup$ Mathematica 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"? $\endgroup$ – Coolwater Dec 10 '17 at 20:48
  • $\begingroup$ @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. $\endgroup$ – Knowledge Dec 11 '17 at 7:34
52
$\begingroup$

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.

$\endgroup$
  • 8
    $\begingroup$ Why does it parse it like that? $\endgroup$ – Rivasa Dec 11 '17 at 15:30
  • 2
    $\begingroup$ 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). $\endgroup$ – ilian Dec 11 '17 at 23:01
  • 2
    $\begingroup$ @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. $\endgroup$ – undercat Dec 12 '17 at 7:29
  • $\begingroup$ 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!! $\endgroup$ – Mehrdad Dec 12 '17 at 10:29
  • 3
    $\begingroup$ @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). $\endgroup$ – ilian Dec 12 '17 at 15:23
16
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ +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. $\endgroup$ – ilian Dec 12 '17 at 17:16
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ – Kuba Dec 12 '17 at 20:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.