baffling... rounding error of not so small numbers

The following results tooked me by surprise.

>> floor(1/1e-3)

ans =

        1000

>> floor(1/1e-4)

ans =

       10000

>> floor(1/1e-5)

ans =

       99999

>> floor(1/1e-6)

ans =

     1000000

Impossible!

Answer 1

[Walter Roberson]

Not impossible. The pattern probably repeats with a period of 3.

http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F

Answer 2

[Geoff]

That's not unexpected at all. By taking the floor, you are pretending that floating point precision limits do not exist. In computing, this can be dangerous.

If you are expecting an answer to be 'very near' a value, you can anticipate it to be out by a very little bit. This value is commonly called 'epsilon'.

You need to allow for floating point limits, and it is often best to overestimate the error. If an error of 1e-6 makes absolutely no difference to your algorithm, then allowing it to fail on an error of 1e-11 is ludicrous.

So let's say you set your epsilon to a generous 1e-8.

Now:

floor(1/1e-5 + epsilon)

will return the correct answer.

There is a function called eps in MatLab, which returns the epsilon for different sized numbers (read help eps). Try this:

>> 1e5 - 1/1e-5 - eps(1e5)

ans =

       0

Here is the limit for epsilon. The precision error for 1e5 just happens to be precisely the error you found in that division.

So, while this gives your expected result...

floor(1/1e-5 + eps(1e5))

it is pushing your error tolerance right to the limit, and

floor(1/1e-5 + eps(1e4))

will give your original result.

Try to use epsilons of one order of magnitude higher than any of your intermediate results. Note also that when performing many operations your error will accumulate, and so you may need an even larger epsilon.

Hope that helps, and I didn't ramble too much.

-g-

[Andrew Newell]

Better yet, use round instead of floor.

[Geoff]

Hehe, well yes but that wasn't the question. =)

There are times when you want to take the floor for use as an index. Let's say you're dividing a UNIX epoch time by 86400 (seconds) to give a day index into some array. And imagine that the zeroth second of a new day sometimes has a precision error that puts it in the previous day.

In this case, you can't use _round_, or everything after midday will get put in the next day. Instead you need some kind of epsilon, which would be maybe 0.5 / 86400.

Maybe this example is contrived, but it's the quickest thing I can think of to illustrate a point. This would apply to enumerating any cyclic data. So... Don't say "use datenum()" =P Actually, for the record, I did run into this very problem doing timezone conversions on MatLab dates using UNIX times.

[Lake-Ee]

I now felt silly asking the question, since the problem is ubiquitous. Geoff gave a good remedy. Nonetheless, it is scary that an average user is probably not aware of this.

Actually, this problem came about from wanting to round a number to the 5th decimal. This (round(X*1e5)/1e5) worked (X=3,6,7,11,...), but not this (round(X/1e-5)*1e-5). I switched to floor just to make the problem more apparent.

[Geoff]

I don't think it's a silly question. You are correct in that many people do not suspect or understand precision limits in computing. At some point though, most people who perform calculations on computers are going to come across this problem, and only some will try to find out the reason for it.