Some problem about ceil()

The result of ceil(215.08/0.152) is 1416. However, it should be 1415 in practice. This makes a further mistake for the latter calculation. How can I avoid this problem?

[Star Strider]

Your image did not appear.

Using round instead of ceil returns 1415.

[Edmond Dantes]

Thanks for your comment. But I really need to use ceil(), also to deal with other cases. The meaning here I need is just the nearest bigger integer, not round.

Answer 1

[James Tursa]

Yeah ... floating point arithmetic bites again:

>> num2strexact(215.08/0.152)
ans =
1.415000000000000227373675443232059478759765625e3
>> num2strexact(215080/152)
ans =
1.415e3

The trailing bits of the floating point calculation will make a difference in the result. If you need the 2nd result in your calculations, your code is not robust against these floating point differences and you will need to rewrite your code. You can see that trailing bit in the hex representation:

>> num2hex(215.08/0.152)
ans =
40961c0000000001

[Image Analyst]

Why does this:

fprintf('%.30f\n', 215.08/0.152);

give a different answer:

1415.000000000000200000000000000000

[James Tursa]

Because I am guessing you are on a Windows machine and not a Unix machine. For some reason, the algorithm behind the "printf" family of functions on a PC is different from the one TMW uses on Unix machines. The PC version only prints out the "significant" digits of the number regardless of how many you request be printed out. In fact, that difference is what inspired me to write the num2strexact mex function in the first place.

[Image Analyst]

Thanks. It could be helpful/speedy/convenient to people if you provided a link to your num2strexact function.

[James Tursa]

http://www.mathworks.com/matlabcentral/fileexchange/22239-num2strexact--exact-version-of-num2str-

[Edmond Dantes]

I can get your meaning, but I still have no idea about how to write my code to avoid this problem. I have some other calculations, such as 300/0.139. I need the result of 2159. So I have to use ceil() to get the results.

[Edmond Dantes]

I have tried some other calculations, ceil(3.002/0.158)=19, ceil(1.807/0.139)=13, ceil(20.216/0.152)=133. They are all OK. Why only ceil(215.08/0.152) failed?

[Star Strider]

I suspect you’re seeing the result of floating-point approximation error.

Example —

a = sym(215.08);
b = sym(0.152);
c = a/b
d = ceil(215.08/0.152)

c =
1415

d =
     1416

The symbolic calculation produces the result you expect, the double-precision calculation does not.

[Walter Roberson]

The error here is in thinking that any floating point number such as .152 has an exact binary representation. 1/10 is an infinite repeating number in binary, just the same way that 1/3 or 1/7 is an infinite repeating number in decimal. If you approximate 1/3 as 0.333 decimal and add that together three times you would get 0.999 not the 1 you would expect from algebra. You are going to suffer from round-off error for any given finite precision when you calculate a fraction whose denominator is not a power of the number base.

So when you write 215.08/0.152 you are not dividing (21508/100) by (152/1000): each of the numerator and denominator are approximated in base 2, giving you values that are not exactly the same as the base 10 values you see displayed. In each case for a literal number, MATLAB picks the representable floating point number which is nearest to the decimal number, which might be less than the decimal number or might be greater than the decimal number. The representation error will propagate after that.

I was going to suggest that you switch to the symbolic toolbox and use, for example, sym('215.08') as I thought that the symbolic toolbox represented such numbers in decimal internally. However, I did some tests this evening and discovered that I was mistaken, that the actual internal representation is binary. For example,

>> sym('.3')*sym(100) - sym(30)
ans =
5.8774717541114375398436826861112e-39

Perhaps https://www.mathworks.com/matlabcentral/fileexchange/6446-multiple-precision-toolbox-for-matlab would be suitable for your work.

[Edmond Dantes]

I think it is different between sym('0.3') and sym(0.3). sym('.3') returns 0.3, but sym(0.3) returns 3/10.

I have tested both in Matlab and C++. If only write 215.08, it is considered as a double-precision number. And as you said, it is an infinite repeating number in binary and will cause round-off error.

If I write ceil(sym(215.08)/sym(0.151)), it will return a right result. However, I seldom saw others using sym in Matlab programs, such as y = sym(a)/sym(b). Maybe it is suitable for this problem, but no need for a common calculation?

Nevertheless, how can I deal with this problem in C++? I have not found a same method as sym in C++.

[Walter Roberson]

sym() is slower for calculations, and a lot of the time the extra precision is not required.

Floating point operations are not transitive or distributive.

Programs that are fragile to single bit round-off are usually not designed with proper numeric error analysis and so tend to break for other reasons. Like failing to recognize that a calculation will overflow or underflow under circumstances that were thought to be handled. (If you use the gamma function, or one of the Bessel-related functions, or the ratio of factorials, or if you use a polynomial of degree higher than seven over a range outside [-1, +1], then chances are your code breaks in ways you did not plan for.)

Writing robust floating point calculations is not easy.

[Edmond Dantes]

Thanks a lot for you all. I will use sym to solve this problem. y = ceil(sym(a)/sym(b)).

Answer 2

[Guillaume]

As others have commented, this is normal behaviour for computers using ieee floating points. Switching to some non-binary storage system (e.g. .Net System.Decimal) is an option. Alternatively, you can round your result to something with less decimal, then use ceil on that:

ceil(round(215.08/0.152, 8))

While the above works for this particular case, I'm sure that some counterexamples could be found.