mod(a,b) floating-point rounding error

Hello, all. I want to generate three vectors with the same range but different "resolutions." That is, the first vector will contain only integers, the second vector will contain integers and 1/10 multiples, the third will include 1/100 multiples. Here is my code:

delta = [1, 0.1, 0.01];
for i = length(delta):-1:1
  if i == 3
      x(i,:) = -1:delta(i):6;
  else
      x(i,:) = x(3,:).*(mod(x(3,:),delta(i))<1e-15);
  end
end

As you can see, I avoid floating-point rounding errors by stipulating the mod() output be less than 1e-15. However, this particular code yields a rather egregious error when building the 1/10 vector. In particular, it gives

mod( x(3,121), delta(2) ) = 0.1

When I'm debugging, the variables are shown to hold the following values:

x(3,121) = 0.2

delta(2) = 0.1

Unless I'm missing something, this is not right.

I'll note that the rest of the vector comes out correct! It's only when x(3,k) = 0.2 that the code outputs the wrong value.

Any help is much appreciated.

Thanks, Jeff

[Geoff Hayes]

Jeff - x(3,121) might not be 0.2 exactly but something close to it. If I run your above code and then call

 fprintf('%.30f\n',x(3,121))

I get an answer of

 0.199999999999999955591079014994

and so

 mod(x(3,121), 0.1) appears to be 0.1

which is what you have.

If I assign 0.2 to x(3,121) then

 x(3,121) = 0.2;
 fprintf('%.30f\n',x(3,121))

is

 0.200000000000000011102230246252

with

 mod(x(3,121), 0.1) appears to be 0

which is what you are expecting or hoping for. See floating point numbers for more information on this topic and why you don't always get the answer you are expecting when using floating point numbers...

[Stephen23]

"mod(a,b) yields incorrect answer (not a floating-point rounding error)"

Actually mod(a,b) yields the correct answer (due to floating-point rounding error).

Because of floating point error, x(3,121) is slightly less than its nominal value, and so mod returns the correct answer. It is not clear how or why you expect mod(...)<1e-15 to make any difference to this, as by the time you make that comparison it is already too late to make any difference. Perhaps you assume that floating point error is always positive, giving values larger than the nominal value, but as this example shows, this is not the case.

One simple solution is to not take mod of floating point values: simply perform your operations on integers, and then afterwards scale down to the magnitude you need.

Answer 1

[OCDER]

Nope, still a tricky rounding error.

isequal(x(3, 121), 0.2) 
ans =  
       0

So x(3, 121) is NOT 0.200. Instead, it could be 0.19999999999999999999999999999999999999999 .

When you do

mod(x(3, 121), 0.2))

you'll get

0.09999999999999999999999999999999999999999

which is why

mod(x(3, 121), 0.2) = 0.1

And so the following statement will return 0:

x(2, 121) = x(3, 121).*(mod(x(3, 121), 0.2) < 1e-15)   % = 0

To fix all this, use ismembertol like this:

delta = [1, 0.1, 0.01];
for i = length(delta):-1:1
  if i == 3
      x(i,:) = -1:delta(i):6;
  else
      GetLoc = ismembertol(x(3, :), -1:delta(i):6, 1e-15);
      x(i,GetLoc) = x(3,GetLoc);
  end
end

Answer 2

[Aquatris]

x(3,121) is not exactly 0.2. It is slightly less than 0.2. That is why when Matlabs mod() function does its algorithm, which is;

mod(x,y) -> x - floor(x./y).*y

The floor(x(3,121)/delta(2)) returns 1 instead of 2. This might be because elements of x3 did not initialized with exactly 0.01 apart (once again numeric issue).

One solution can be to use symbolic toolbox and replace your code with (I do not know what you are trying to do with 1e-15 section but you cannot use it with symbolic) ;

delta = [1, 0.1, 0.01];
for i = length(delta):-1:1
  if i == 3
      x(i,:) = sym(-1:delta(i):6);
  else
      x(i,:) = sym(x(3,:).*(mod(x(3,:),delta(i))));
  end
end

Answer 3

[David Goodmanson]

Hi Jeff,

This is definitely a floating point issue. After running the code,

x(3,121) -.2
ans = -5.551115123125783e-17

In this situation mod (... , .1) has the effect of rounding down to .1 so the test fails and you get 0. You could come up with a modified test, but it makes a lot more sense to avoid working in floating point altogether. For example

x = -1:.01:6;              % as before
n = length(x);
x = [zeros(2,n);x];
x(2,1:10:n) = x(3,1:10:n);
x(1,1:100:n) = x(3,1:100:n);