The “malicious” rounding that you are attempting to protect against cannot actually occur in IEEE-754 floating-point. You are jumping at shadows, making lots of results less accurate to protect against something illusory.
Why is this the case?
In order for ceil(L/M) to produce an incorrect result, we would need to have the mathematically exact result be of the form N + epsilon, with 0 < epsilon <= ulp(N)/2, and N an exact integer. But, can that ever actually happen? No. If it did, we would have a floating-point number L = M(N+epsilon) exactly, or L = MN + delta, with 0 < delta < ulp(MN)[1], which can never happen by the definition of ulp.
This illustrates a difficulty naive backwards-error analysis as commonly taught in numerical analysis courses. While it provides reasonable conservative bounds, which is often sufficient for engineering needs, it misses much of the subtlety of floating-point arithmetic, and leads people to take unnecessary steps to guard against “errors” that aren’t really there.
[1] Lemma: For any binary floating-point type, ulp(xy) > x ulp(y)/2.
Proof: First note that rounding of the product xy will only ever make ulp(xy) larger, so we can safely ignore it. Use the formula ulp(x) = 2^(floor(log2(x)) - P) where P is the precision in bits and expand:
I suspect a bug in your proof right at the beginning, where you say "exact result be of the form N + epsilon, with 0 < epsilon <= ulp(N)/2, and N an exact integer". Assuming you're referring to the result of the entire expression including ceil, can we still say that, considering that ceil() is far from being a continuous function?
Anyway, my example was a bit too theoretical. A better version is one where L and M are integers which are not necessarily representable. Here: http://ideone.com/1Xkhjp
Yes, of course it’s possible to exhibit errors if you change the problem (though it should be noted that the source of error here is entirely in the fact that l is not representable). As soon as you allow representation error, of course, your tolerance-based approach goes out the window; what if the representation error in l and m is larger than your tolerance accounts for?
If I wanted a correct bound for your revised problem, I would require that l be an upper bound on the true value (by forcing it to be rounded up earlier computations, for example) and that m be a lower bound on the true value. Then the simple computation delivers a correct result regardless of the magnitude of the error in l or m.
It is the technique of backwards-error analysis that is naive, not the resulting bound. A “naive” method is not “simple” or “dumb” (or “aggressive” or whatever you think the opposite of “conservative” might be). It is simply a method that doesn’t use some a priori knowledge that’s available for the specific problem at hand. Backwards-error analysis is a naive technique in that it doesn’t benefit from specialized theorems of floating-point that allow one to establish tighter bounds, like the one I proved.
Comments
The “malicious” rounding that you are attempting to protect against cannot actually occur in IEEE-754 floating-point. You are jumping at shadows, making lots of results less accurate to protect against something illusory.
Why is this the case?
In order for ceil(L/M) to produce an incorrect result, we would need to have the mathematically exact result be of the form N + epsilon, with 0 < epsilon <= ulp(N)/2, and N an exact integer. But, can that ever actually happen? No. If it did, we would have a floating-point number L = M(N+epsilon) exactly, or L = MN + delta, with 0 < delta < ulp(MN)[1], which can never happen by the definition of ulp.
This illustrates a difficulty naive backwards-error analysis as commonly taught in numerical analysis courses. While it provides reasonable conservative bounds, which is often sufficient for engineering needs, it misses much of the subtlety of floating-point arithmetic, and leads people to take unnecessary steps to guard against “errors” that aren’t really there.
[1] Lemma: For any binary floating-point type, ulp(xy) > x ulp(y)/2.
Proof: First note that rounding of the product xy will only ever make ulp(xy) larger, so we can safely ignore it. Use the formula ulp(x) = 2^(floor(log2(x)) - P) where P is the precision in bits and expand:
I suspect a bug in your proof right at the beginning, where you say "exact result be of the form N + epsilon, with 0 < epsilon <= ulp(N)/2, and N an exact integer". Assuming you're referring to the result of the entire expression including ceil, can we still say that, considering that ceil() is far from being a continuous function?
I’m referring to the mathematically exact result L/M. Sorry for being unclear.
Anyway, my example was a bit too theoretical. A better version is one where L and M are integers which are not necessarily representable. Here: http://ideone.com/1Xkhjp
Yes, of course it’s possible to exhibit errors if you change the problem (though it should be noted that the source of error here is entirely in the fact that l is not representable). As soon as you allow representation error, of course, your tolerance-based approach goes out the window; what if the representation error in l and m is larger than your tolerance accounts for?
If I wanted a correct bound for your revised problem, I would require that l be an upper bound on the true value (by forcing it to be rounded up earlier computations, for example) and that m be a lower bound on the true value. Then the simple computation delivers a correct result regardless of the magnitude of the error in l or m.
"This illustrates a difficulty naive backwards-error analysis"
It's not naive. It's conservative, the exact opposite.
I pointed part of this out:
It is the technique of backwards-error analysis that is naive, not the resulting bound. A “naive” method is not “simple” or “dumb” (or “aggressive” or whatever you think the opposite of “conservative” might be). It is simply a method that doesn’t use some a priori knowledge that’s available for the specific problem at hand. Backwards-error analysis is a naive technique in that it doesn’t benefit from specialized theorems of floating-point that allow one to establish tighter bounds, like the one I proved.