bug in multiplication operations

this code returns wrong result

disp(string(1899243482440473*5+1))

9496217412202364

Why? It's a big problem. Correct result is 9496217412202366

Answer 1

[Steven Lord]

Not all integer values that are in range of double precision are exactly representable as doubles. If we look at your constant:

x = 1899243482440473*5+1

x = 9.4962e+15

the spacing between representable double values for numbers of that magnitude is:

distanceToNext = eps(x)

distanceToNext = 2

This is not a bug in MATLAB. This is a consequence of how numbers are stored in IEEE double precision. See section 1.7 in the Introduction to MATLAB chapter of Cleve Moler's Numerical Computing with MATLAB for more information about floating point.

If you're using double precision numbers that are larger than flintmax you're going to encounter this behavior.

format longg
flintmax

ans = 
      9.00719925474099e+15

If your numbers are larger than this but smaller than the maximum integer that a 64-bit integer can represent, using int64 or uint64 as @Vilém Frynta suggested is a possibility. If they're larger than intmax for the 64-bit integer types, you're going to need to use Symbolic Math Toolbox or another arbitrary-precision package.

rangeOfUint64 = [intmin('uint64') intmax('uint64')]

rangeOfUint64 = 1×2

0 18446744073709551615

rangeOfInt64 = [intmin('int64') intmax('int64')]

rangeOfInt64 = 1×2

-9223372036854775808 9223372036854775807

If you do use 64-bit integers, be careful not to perform computations in double then convert if the results of the computations could be greater than flintmax. @Vilém Frynta was careful to do this in their suggestion.

Answer 2

[Vilém Frynta]

There's probably a small error due to number being represented by bits. Number is too long and detailed, and Matlab is having a problem to accurately describe it.

Use int64.

disp(string(int64(1899243482440473)*5+1))

9496217412202366

[Taras]

I can't use Int64, because I'm working with large numbers.

Are you going to fix this problem?

[Vilém Frynta]

Other thing you can do is use symbolic toolbox or use File Exchange functions.

syms x
y = x*5 + 1;
result = vpa(subs(y, x, 1899243482440473));
disp(result);

9496217412202366.0

[Taras]

Thank you @Vilém Frynta for your help.

But it's only a workaround, not a fix. Such an approach imposes restrictions on me.

[Vilém Frynta]

yes, it's only a workaround. i'm afraid i cannot help futher, i'll leave the space to someone else who could assist you better. perhaps a Mathworks employee will appear.

if you found my answer atleast a bit helpful, i'd be happy if you could accept it.

[Rik]

What you want is not possible with the native datatype of Matlab (the IEEE-defined double precision floating point value):

flintmax<1899243482440473

ans = logical
   0

The number you want to use is larger than the largest number where x+1 can be represented. If I give you 5 characters, the largest number you can uniquely represent is 99999. If I tell you to store 158021, you will have to round to 1.6e4.

The only work-around is to use a different data type. @Taras already suggested vpa. You might be interested in reading this page if you want to use quad precision.

[Walter Roberson]

"Are you going to fix this problem?"

Frankly:

No, Mathworks is not going to fix this problem. This problem is inherent in any fixed-width floating point number system, no matter what the number base (base 2, base 10, base 60...). If Mathworks were to switch to 96 bit or 128 bit floating point numbers, the same problem would still exist, just at different ranges.

To "fix" the problem, Mathworks would have to switch to variable-precision numbers. There has not been hardware support for variable precision numbers since roughly the days of the IBM 1610, so calculations would all have to be done in software... pretty much the Symbolic Toolbox. Much much slower.

[Walter Roberson]

There are other inherent consequences to fixed-width mathematics that cannot be "fixed" except by not using fixed-width mathematics.

Suppose you are working in 5 digit decimal, and you calculate 1/3 = 0.33333 . Add three of those together and you get 0.99999 . But 3 * (1/3) = 1 not 0.99999 . Should 5-digit decimal have represented 1/3 as 0.33334 ? Add three of those together and you get 1.00002 .

Would working in more decimal digits solve the problem? Suppose you work in 20 decimal digits, then 1/3 -> 0.33333333333333333333 and add three of those you get 0.9999999999999999999999999 not 1 exactly. If you worked with 1000 decimal digits then adding three of the representation of 1/3 would get you 1000 nines, not 1 exactly. Decimal has no finite representation of 1/3 . Switching to binary does not solve the problem: binary has no finite representation of 1/3 either. If you were to switch to base 60, then yes base 60 has an exact representation of 1/3 (as 20 / 60)... but finite base 60 representation cannot exactly represent 1/7, 1/11, 1/13, 1/17...

No matter what integer numeric base you choose, base N, base N has no finite representation for the number 1/(N+1)

The "bug" is in finite mathematics. And the only two ways to fix it are:

1. Change the universe so that mathematics works differently; or
2. Do not use finite mathematics.

I came across a paper a number of years ago that showed that there is a family of functions that if you calculate to any specific finite precision, that the error in the answer relative to the actual answer can be arbitrarily large. And that is a problem that affects variable-precision mathematics, where the user is permitted to specify the precision of calculations. Besides that being slow, it can give arbitrarily bad results (relative to the actual results, not in an absolute sense.) So variable-precision calculations will have "bugs" too -- bugs that cannot be fixed in any possible implementation.