Modulus Problem When Using Floating-Numbers

Hello everyone,

I would like to compute the following formula using MATLAB's mod function.

My first five x0 values are 0,360000000000000 0,921600000000000 0,289013760000000 0,821939226122650 0,585420538734197

M value is fixed and is equal to 256.

When i use mod function directly, i get wrong results such as [255,992187500000 0 0,00781250000000000 40,9843750000000 11,7421875000000] which are not even integers.

Here x_0x10^14 is less than flintmax which is 2^53 so probably it is a precision problem.

Can you help me to solve my problem?

Answer 1

[John D'Errico]

This is only precision problem, if you accept that means a problem of understanding a floating point arithmetic environment.

In MATLAB (as is true in many languages), the number 0.36 is not represented exactly as the fraction 36/100. Instead, MATLAB uses a binary form to represent the number.

x = 0.36

x = 0.3600

It LOOKs like 0.36. But is it? Just like the number 1/3 is NOT representable in a finite number of digits in a decimal form, the fraction 36/100 is not representable in a finite number of binary bits. MATLAB is forced to represent that number as effectively

sum(2.^[-2 -4 -5 -6 -11 -13 -17 -18 -19 -20 -22 -24 -25 -26 -31 -33 -37 -38 -39 -40 -42 -44 -45 -46 -51 -53])

ans = 0.3600

which still looks vaguely like the number you wanted. But is it that? In fact, it should be an infinitely repeating number, with extra binary bits down at the botttom, just like 1/3 or 2/3 cannot be represented in decimal form.

vpa(sum(sym(2).^[-2 -4 -5 -6 -11 -13 -17 -18 -19 -20 -22 -24 -25 -26 -31 -33 -37 -38 -39 -40 -42 -44 -45 -46 -51 -53]) - sym('0.36'))

ans = 

Now, what happens when you multiply 0.36 with 1e14? It is not an integer. That would only be true if 0.36 was exactly representable. And therefore mod gets in trouble. It does what it can. But you passed it a non-integer, so what do you expect? The same is true of each of your other cases.

If you want a valid result, then you could use tools like syms (or my own HPF toolbox) to do all of your work. But expect things to be really incredibly slow.

mod(sym('0.36')*1e14,256)

ans = 

0.0

My HPF will work, because it stores the numbers in a decimal form, and therefore 0.36 is not stored as an approximation.

mod(hpf('0.36')*1e14,256)
ans =
    0

In the end, you will be far better off understanding floating point arithmetic than to use a powerful tool like that though.