mod gives incorrect result

Here is a mathematical error with a simple expression:

>> mod(-sin(pi),512)

ans =

   512

Of course, pi isn't exact, so sin(pi) gives a small negative number. That isn't good, but it is perhaps excusable. However, that mod(X, 512) ever returns 512 is totally unacceptable.

[Steve]

I don't want to use rem because I want mod(-1,512)=511. I really do want mod, but I want 0 <= mod(x,512) <5 12. I think anyone using the mod operator would want that.

I don't know if there is a better fix, but I am doing the following workaround (with floor to convert to an integer index which must be in the range 0 to N-1.):

i = floor(mod(x, N)); if (i == N) i = 0; end;

[the cyclist]

Another workaround would be to apply the mod function twice. [This is another indication that mod()'s output is a bit mathematically silly in this case.]

Answer 1

[Walter Roberson]

In order to understand this, you need to look at the fact that -sin(pi) is a negative value that has a very small magnitude (due to standard round-off issues), so mod(-sin(pi), 512) is 512-sin(pi) . But sin(pi) is less than eps(512), so 512-sin(pi) is most closely representable as 512 itself.

This is, I agree, not what I would have expected, but I had never really thought about the matter before.

Would it be better if mod(x, P) with positive P and x in (-eps(P), 0) (exclusive) be 0? Arguable, but it is not obvious to me that that possible outcome would be better than the current outcome.

[the cyclist]

In essence, what is happening here is that the tiny floating point error is going through an "amplifier". It is analogous to expecting

(1 - 2/3 - 1/3)*1e17

to be 0, but it turns out to be about 5.5.

So, there is the usual aspect of "let the buyer beware" if floating point error can give a big difference in your results.

In this particular case, though, I would argue that 512 as the "most closely representable" is less compelling, because 512 is not in the theoretical range of the mod() function, mathematically. I know it sounds odd, but I think that 0 is the most closely representable number to 512-eps in this case. Mathematically,

Limit   mod(512-eps,512)
{eps->0}

is zero, not 512.

Answer 2

[Wayne King]

Are you sure you're not confusing mod() with rem()?

If you read the help for mod(), it says that

 mod(x,y) returns  x-floor(x/y)*y

 rem(x,y) returns x-fix(x/y)*y

   -sin(pi)-floor(-sin(pi)/512)*512

is 512 because

   floor(-sin(pi)/512)

is -1.

so essentially 0+(1)*512 = 512

on the other hand

   -sin(pi)-fix(-sin(pi)/512)*512

yields essentially

 0-(0)*512=0

[the cyclist]

It's admirable that MATLAB is doing what it's documentation states, but mathematically

mod(x,N)

has a range of [0,N), strictly less than N.

[Steve]

Thanks for your reply.

Yes, if -sin(pi) isn't zero, then floor(-sin(pi)/512) is -1.

Although you write "essentially 0", it seems unwarranted to treat -sin(pi) as unequal to zero in "floor(x/y)" and equal to zero as "x". Rather than "0+(1)*512", as your write, isn't it -sin(x)-(-1)*512, which should be slightly less than 512.

[Wayne King]

I think "essentially 0" is warranted here since -sin(pi) is on the order of 10^(-16). For example:

   rng default;
   x = randn(10,1);
   max(abs(ifft(fft(x))-x))

gives a value of 8x10^(-16) but nobody would say that the example above does NOT demonstrate perfect inversion of x.

Answer 3

[Wayne King]

If anybody's interested in this, google:

"division and modulus for computer scientists"

mod() is implementing what Knuth termed "floored division". It's also described here:

http://en.wikipedia.org/wiki/Modulo_operation

Just FYI, Python implements the same operation for % as MATLAB mod(), but R appears to be different. R's %% operator appears to mimic rem() in MATLAB.

So in Python:

>>>import Numpy as np
>>>-np.sin(np.pi)%512

returns 512.

In R

> -sin(pi)%%512

returns zero like MATLAB

rem(-sin(pi),512)