Why do I get a wrong output when using REM in MATLAB for very big numbers?

I get an erroneous output when I execute the following:

rem(6.3055e16, 7)

The command above returned the value 8 which is larger than the divisor.

Answer 1

[MathWorks Support Team]

This output might be attributed to the machine precision error. Machine precision is defined as the smallest difference between two numbers that the computer recognizes. In other words, on any computer, there is a small gap between each double-precision number and the next larger double-precision number. We can determine the size of this gap, which limits the precision of results, using the EPS function. For example, to find the distance between 6 and the next larger double-precision number, we enter:

format long

eps(6)

This tells us that there are no double-precision numbers between 6 and 6 + eps(6).

Any number larger than 4.5036e+015 or of similar order has an EPS greater than 1. If we have a number whose EPS is greater than 1, we cannot guarantee the output of functions like MOD and REM. Although it might seem like it is a bug, it is not a bug. This is an expected behavior.

We have two workarounds for this problem:

First, if we know that the first argument can be written in the form "a^b", where both "a" and "b" are smaller than something of the order 4.5036e+015, then there is code at MATLAB Central that uses an RSA-type algorithm to get the correct answer. The same can be found at:

 <http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=7908>

Second, we can use a longhand approach, for example,

A = 6.3055e16;

B = A - 7*9e15; % This number is a multiple of 7, so it doesn't affect the REM or MOD functions with base 7

rem(B,7) % REM works properly now

A better explanation to the stated problem is given as follows:

Referring to the help document for REM, the formula is given as

x - n.*y where n = fix(x./y).

Now we do the calculation at the MATLAB prompt directly to see what happens:

x = 6.3055e16; y = 7;

n = fix(x/y)

dif = x - n*y

we receive the following output:

%%%BEGIN PRE%%

n =

9.0079e+015

dif =

     8

%%%END PRE%%

This is the same answer as REM got. Now we note that

>> eps(n)

ans =

     2

>> eps(x)

ans =

     8

>> eps(n*y)

ans =

     8

So we are subtracting two floating point numbers, and the epsilon for each operand is 8. We cannot hope to get a result less than 8, except possibly 0. In fact, since eps(x/y) = 2, the call to fix did not do anything

>> x/y == fix(x/y)

ans =

     1

and we are really just evaluating

   >> x - (x/y)*y

ans =

     8