A common and simple way to treat FP mathematically is by bounding the result with the machine epsion u [1]. For example, given FP numbers x and y:
x <fp_mul> y = (1+e)xy ; abs(e) <= u
The u ideally depends only on the FP format and rounding mode. I say ideally because library functions like sqrt may be less precise than the ideal (which is the result being equal to the exact value rounded according to the rounding mode). Note that the abstraction only as long as you don't go outside the range of normal exponents (it breaks when you reach infinity or subnormal).
EXAMPLE. Suppose we have a FP value L representing a length of something, and want to split it into segments no longer than the FP value M. Assume both integers are small enough to be representable in FP. The answer in exact arithmetic is:
N = ceil(L/M).
However if we evaluate this in FP, we may get a result N'<N for certain edge values. Which means one of the segment lengths we output will be greater than M!. So the result of the algorithm is incorrect!
We can fix it by multiplying (L/M) with a number sufficiently greater than one before the ceil, at the cost of sometimes deciding to split into one more segment than mathematically necessary. The fixed FP algorithm may then be:
N' = ceil(1.00005 <fp_mul> (L <fp_div> M))
After we apply the machine epsilon formulas we come to the conclusion that:
N' = ceil( (1+e1)(1+e2)1.00005 (L/M) ) ; abs(e1)<=u, abs(e2)<=u
We then show that:
(1+e1)(1+e2)*1.00005 >= 1.
This proves that what the algorithm gives to the ceil function is not less the exact value L/M, and that N'>=N. And ceil() can't make errors by definition, assuming the integer value is representable. QED
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
A common and simple way to treat FP mathematically is by bounding the result with the machine epsion u [1]. For example, given FP numbers x and y:
x <fp_mul> y = (1+e)xy ; abs(e) <= u
The u ideally depends only on the FP format and rounding mode. I say ideally because library functions like sqrt may be less precise than the ideal (which is the result being equal to the exact value rounded according to the rounding mode). Note that the abstraction only as long as you don't go outside the range of normal exponents (it breaks when you reach infinity or subnormal).
EXAMPLE. Suppose we have a FP value L representing a length of something, and want to split it into segments no longer than the FP value M. Assume both integers are small enough to be representable in FP. The answer in exact arithmetic is:
N = ceil(L/M).
However if we evaluate this in FP, we may get a result N'<N for certain edge values. Which means one of the segment lengths we output will be greater than M!. So the result of the algorithm is incorrect!
We can fix it by multiplying (L/M) with a number sufficiently greater than one before the ceil, at the cost of sometimes deciding to split into one more segment than mathematically necessary. The fixed FP algorithm may then be:
N' = ceil(1.00005 <fp_mul> (L <fp_div> M))
After we apply the machine epsilon formulas we come to the conclusion that:
N' = ceil( (1+e1)(1+e2)1.00005 (L/M) ) ; abs(e1)<=u, abs(e2)<=u
We then show that:
(1+e1)(1+e2)*1.00005 >= 1.
This proves that what the algorithm gives to the ceil function is not less the exact value L/M, and that N'>=N. And ceil() can't make errors by definition, assuming the integer value is representable. QED
[1] http://en.wikipedia.org/wiki/Machine_epsilon
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.