How does the mod function compensate for floating point round off?

The documentation for the mod function states, "mod attempts to compensate for floating-point round-off effects to produce exact integer results when possible." What exactly does this mean?

Answer 1

[Joel Handy]

Answering my own question. I wrapped the mod function and generated c code. if the quotient is within eps * quotient of bieing an integer value, the mod function returns zero.

[John D'Errico]

x = 1 + eps;
mod(x,1)

ans = 2.2204e-16

Does that suggest your answer may have a subtle flaw? There is no number repersentable in MATLAB between 1 and (1 + eps).

[Bobby Cheng]

Two things:

1) The compensation only happens when the second input is not of integer values. That is why mod(1+eps,1) returns eps. Otherwise it is what Joel said.

2) I can see there is room for improvement for MATLAB documentation to better clarify this.

Answer 2

[Matt J]

I suspect it is just meant to tell you that it doesn't do a naive computation like in modNaive below.

m=1000*pi; n=pi;

mod(m,n)

ans = 0

modNaive(m,n)

ans = 3.1416

It isn't something you should try to rely on, though. There are definitely cases where the result won't be an exact integer, even when it's clearly what would be ideal, e.g.,

mod(117*sqrt(1001)+1,sqrt(1001))-1

ans = -3.5527e-14

function out=modNaive(m,n)

  out=m-floor(m/n)*n;

end

[Joel Handy]

Thanks for the reply, but I should give a little more background information. I understand the implications to the results of the mod function, but not the specific conditions under which the mod function will coerce the result.

I've used matlab coder and embedded coder to generate code for an algorithm. There is a section of code that, for various reasons, I cannot generate code for, and needs to be hand coded. The section in question contains the mod function.

I'm attempting to test the C++ version against the matlab version. My hand coded mod function produces results along the lines of modNaive which unfortunately has the effect of missaligning all the output data as well producing slightly different results.

in the end, I'm looking to code a version of mod that mataches matlab.

plus I think its important to fully understand any quirks of the functions I use and the hand wavey explanation in the mod documentation doesnt sit well with me.

[Matt J]

in the end, I'm looking to code a version of mod that mataches matlab.

Matching Matlab will serve no purpose. Even if you were to match its implementation, mod is inherently a discontinuous function, which means that if you run your code on different CPU architectures, or with differences in multithreading, you will always run the risk of getting different results in the neighborhood of the discontinuity.

plus I think its important to fully understand any quirks of the functions I use and the hand wavey explanation in the mod documentation doesnt sit well with me.

The hand-wavey explanation likely means that the precise implementation of mod is not something MathWorks wants to commit to. For all we know, it has been different in every Matlab release.

Again, though, it shouldn't matter. There is no implementation of mod that can be counted on to give the same results for all inputs on all machines. The only circumstance in which you can hope for stable results is if you happen to know the inputs will always be integers (in which case Bruno's implementation would probably work well).

[Bruno Luong]

"mod is inherently a discontinuous function, which means that if you run your code on different CPU architectures"

I 'm not sure I can agree, mod must be identical sequence of instructions based on IEEE 754 division and remainder. Any CPU architecture should give identical result.

As for release dependent it seems mod code ghas changed since 2016 if I believeto the H1 line.

[Joel Handy]

Matching Matlab will serve no purpose. Even if you were to match its implementation, mod is inherently a discontinuous function, which means that if you run your code on different CPU architectures, you will always run the risk of getting different results in the neighborhood of the discontinuity.

I agree that mod is inheritantly discontinuous and expect different architectures to give slilghtly different results. I disagree that matching this feature serves no purpose. This isnt exclusively a machine precision issue. This is additional logic. It is intentional manipulation of your results.

[Matt J]

@Joel Handy It is additional logic, and we don't know authoritatively what it is. However, the only thing it could rationally be doing is assessing a tolerance on the error,

err=abs(m/n-round(m/n))

If err<tolerance, it is assumed that m/n is intended to be an integer and the output of mod defaults to zero. Otherwise, mod is computed as in modNaive.

[Matt J]

I 'm not sure I can agree, mod must be identical sequence of instructions based on IEEE 754 division and remainder. Any CPU architecture should give identical result.

That may apply to the operations executed within mod itself. However, the preceding computations that generated the inputs to mod may have different floating point errors affecting them. That would be enough to generate different results in the neighborhood of mod's discontinuities.

Answer 3

[Bruno Luong]

I try to replicate MATLAB mod ith this function, it seems working well for 2 examples, no warranty beside that.

x=(rand(1,1e6)-0.5)*(10*pi);
all(rmod(x,pi) == mod(x,pi))

ans = logical
   1

x = (-1000:1000)*pi;
all(rmod(x,pi) == mod(x,pi))

ans = logical
   1

function r = rmod(x, a)
k = round(x/a);
r = x - a*k;
r(abs(r) < eps(x)) = 0; % EDIT test probably non effective
r(r < 0) = r(r < 0) + a;
end

Answer 4

[Bruno Luong]

Second version, correction under stricter condition than in the first version

function r = rmod2(x, a)
k = round(x/a);
r = x - a*k;
r(abs(r) < eps(a)) = 0;
r(r < 0) = r(r < 0) + a;
end

[Bruno Luong]

I think my test with espsilon is actually useless, this seem to match mod, unless someone can come with a counter example

function r = rmod(x, a)
k = round(x/a);
r = x - a*k;
r(r < 0) = r(r < 0) + a;
end