Incorrect answer in multiplication of doubles

I have a vector of 14 doubles: A = [32 27 25 7 11 13 17 19 23 29 31 37 41 43]. The Matlab function prod gives me the following result: >>prod(A) ans = 9419588158802421760 Let N = prod(A) I get from Matlab: >>mod(N,27) ans = 25 How is it possible? N is the product of 14 factors, one of which is 27 but mod(N,27) is 25? Also the following logical comparison returns true: >>N == (N+200) ans = 1 I am really puzzled by this. Note that N is a double and its value is far less than realmax so I expect it should compute the correct value.

[Stephen23]

"its value is far less than realmax"

Yes, but realmax is irrelevant. double can only represent every consecutive integer up to flintmax.

Answer 1

[John D'Errico]

However, realmax is NOT the limit of a double such that it will represent an exact integer. That limit is 2^53. Your product exceeds that limit.

For example, we see that when a double exceeds 2^53 (called flintmax), the double is not represented down to the least significant bit. For example, we can add 1 to 2^54, and get the same number.

2^54 == 2^54 + 1
ans =
  logical
   1

You could have used a big integer form, for exampel, my VPI toolbox.

A = vpi([32 27 25 7 11 13 17 19 23 29 31 37 41 43]);
mod(prod(A),27)
ans =
    0

Or use sym.

A = sym([32 27 25 7 11 13 17 19 23 29 31 37 41 43]);
mod(prod(A),27)
ans =
0

In either case, you needed a tool that could handle the large integers involved.

Since the product of your numbers is less than the limit of uint64, one might think you could use uint64 to do this.

A = uint64([32 27 25 7 11 13 17 19 23 29 31 37 41 43]);
class(prod(A))
ans =
double

So it looks like the prod of a uint64 array is still returned as a double.

One other thing you could have done, which works since many of the elements in A are at least 27, or larger:

A = [32 27 25 7 11 13 17 19 23 29 31 37 41 43]
mod(prod(mod(A,27)),27)
ans =
     0

In fact, since one element is exactly 27, this too will work:

prod(mod(A,27))
ans =
     0

[Stephen23]

@John D'Errico: note that many arithmetic operations, such prod, sum, etc, also allow a 'native' option, which returns an array of the same class as the input (and presumably does the calculation with that class too):

>> A = uint64([32 27 25 7 11 13 17 19 23 29 31 37 41 43]);
>> prod(A,'native')
ans = 9419588158802421600

[John D'Errico]

Thanks Stephen. There is always something new that we can learn.