CODEWITHSUNDEEP
matlabfloating-pointprecision
nbro
nbro

Reputation: 15837

Subtracting two equal floating-point numbers in MATLAB is not equal 0

Suppose we want to calculate (a + b)2 in two different ways, that is

  1. (a + b) * (a + b)

  2. a2 + 2 a b + b2

Now, suppose a = 1.4 and b = -2.7. If we plug those two numbers in the formulas with format long we obtain in both cases 1.690000000000001, that is, if I run the following script:

a = 1.4; 
b = -2.7;

format long

r = (a + b) * (a + b)

r2 = a^2 + 2*a*b + b^2

abs_diff = abs(r - r2)

I obtain 

r = 1.690000000000001

r2 = 1.690000000000001

abs_diff = 6.661338147750939e-16

What's going on here? I could preview different results for r or r2 (because Matlab would be executing different floating-point operations), but not for the absolute value of their difference.

I also noticed that the relative error of r and r2 are different, that is, if I do

rel_err1 = abs(1.69 - r) / 1.69
rel_err2 = abs(1.69 - r2) / 1.69

I obtain

rel_err1 = 3.941620205769786e-16

rel_err2 = 7.883240411539573e-16

This only makes me think that r are not actually the same r2. Is there a way to see them completely then, if they are really different? If not, what's happening?

Also, both relative errors are not less than eps / 2, does this mean that an overflow has happened? If yes, where?

Note: This is a specific case. I understood that we're dealing with floating-point numbers and rounding errors. But I would like a to understand better them by going through this example.

Upvotes: 0

Views: 419

Answers (2)

Amro
Amro

Reputation: 124543

Don't rely on the output of format long to conclude that two numbers are equal...

a = 1.4;
b = -2.7
r1 = (a + b) * (a + b);
r2 = a^2 + 2*a*b + b^2;
r3 = (a+b)^2;

Instead you can check their hex representation using:

>> num2hex([r1 r2 r3])
ans =
3ffb0a3d70a3d70d
3ffb0a3d70a3d710
3ffb0a3d70a3d70d

or the printf family of functions:

>> fprintf('%bx\n', r1, r2, r3)
3ffb0a3d70a3d70d
3ffb0a3d70a3d710
3ffb0a3d70a3d70d

or even:

>> format hex
>> disp([r1; r2; r3])
   3ffb0a3d70a3d70d
   3ffb0a3d70a3d710
   3ffb0a3d70a3d70d

Upvotes: 3

Ander Biguri
Ander Biguri

Reputation: 35525

Floating point arithmetic is not associative.

While mathematically these two are equal, they are not in floating point maths.

r = (a + b) * (a + b)

r2 = a^2 + 2*a*b + b^2

The order operations are executed in floating point maths is very relevant. That is why when you do floating point maths you need to be very careful of the order of you r multiplications/divisions, specially when working with very big numbers together with really small numbers.

Upvotes: 1

Related Questions