Last update: October 2006
Optimized only for Firefox 3 on CRT!
Contents:
Arithmetic with error bounds is usually implemented with a major weakness. We can illustrate this problem with a simple but artificial example. Suppose we want to compute the expression: "x / x", given the value of x.
Since the numerator and denominator are the same variable, whatever value one of them has, the other has the same one, and the answer should be something like: |1.0,0.0|, no matter what the error range is.
Unfortunately, the common approach assumes implicitly that all the variables are independent, which is obviously not true if variables occur more than once in any expression. If x = |1.0,0.1| then x/x may be computed to be |1.0,0.2|.
The example above is artificial, but "(x + 1) / (x + 2)" is not, and it suffers from the same problem. Other examples occur all the time in practical computing.
It seems that the correlated term problem is ignored by software developers. Note the problem will appear between any dependent variables and that's pretty hard to control.
Can we offer a general polynomial (rational) function that takes care of the said problem? Is it possible?
We already do it for powers. It's easy when the function is monotonic but polynomials are not monotonic in general. They are monotonic between local extermum points and we can find these points using a quick iterative algorithm.
We can have an approximate result using the first terms of a power series developed in the middle.
We can divide the interval by N points and take the largest and smallest values. Actually we can find approximate median.
A = (a,b,c,...) x = |u,v|
POLY(x, A) = a + b*x + c*x^2 + ...
What happens near a singularity?
Suppose we are computing "1/x", and the error range of "x" includes zero. The error range of the result is the union of two segments: "(-infinity,-a]", and "[b,+infinity)".
The monster above is the complement of the interval "[-a,b]". We could modify our algorithm to deal with such creatures, however, this wouldn't close interval operations, as more complicated error ranges may be produced, e.g. when such ranges meet another singularity.
In this package we are using only interval error ranges and flag an error when a singularity is encountered.
SIN is our primitive. We will reduce COS and TAN to SIN. There is an alternative primitive implementation of TAN, but it is not used. Using TAN as the primitive makes things simpler but may be less desirable numerically.
Some facts about SIN:
The number of extermum points and their types determines the result. We could count the maxima and minima in our interval using the fact that every maximum has the form pi/2 + 2*n*pi and every minimum -pi/2 + 2*n*pi.
x%lower <= 2*m*pi + pi/2 <= x%upper (maxima) x%lower <= 2*n*pi - pi/2 <= x%upper (minima) 1/2 * (x%lower/pi - 1/2) <= m <= 1/2 * (x%upper/pi - 1/2) 1/2 * (x%lower/pi + 1/2) <= n <= 1/2 * (x%upper/pi + 1/2) eps = epsilon(ZERO) max1 = ceiling(1/2*((x%lower - eps)/pi - 1/2)) max2 = floor(1/2*((x%upper + eps)/pi - 1/2)) min1 = ceiling(1/2*((x%lower - eps)/pi + 1/2)) min2 = floor(1/2*((x%upper + eps)/pi + 1/2)) nmax = max2 - max1 + 1 nmin = min2 - min1 + 1
Using m and n we can easily write the SIN implementation. The problem with this method is of course the rounding error but since we are interested in rounding outwards we are compensating with the eps factors (which may be increased if necessary).
To illustrate that an approach dividing the problem into a lot of cases is no better than counting minima/maxima, we will persue such a method.
let us define for convenience:
x = [x%lower, x%upper] len(x) = x%upper - x%lower D(x) = cos(x%lower) * cos(x%upper)
The numerically interesting range is len(x) < pi. If len(x) >= pi the computation has already gone bad.
In this mini range there there can be at most 1 extermum. if D(x) > 0 then sin has no extermum in interval x and is monotonic. the image of x is:
y%lower = min(sin(x%lower),sin(x%upper))
y%upper = max(sin(x%lower),sin(x%upper))
If D = 0 there is 1 extermum in an endpoint but sin is still monotonic and the above relations hold. One side will be either +1 or -1.
If D < 0 there is an extermum. the signs of D's factors will indicate the type.
If cos(x%lower) > 0 we have a maximum:
y%lower = min(sin(x%lower),sin(x%upper))
y%upper = 1
If cos(x%lower) < 0 we have a minimum:
y%lower = -1
y%upper = max(sin(x%lower),sin(x%upper))
If cos(x%lower) = 0 we have to check the sign of cos(x%upper) which can't be zero.
If cos(x%upper) > 0 we have a minimum ......
What about the case len(x) >= pi ? Sin is periodic with period 2*pi, so if the error range is larger we must have both a minimum and maximum inside and the result is [-1,+1].
The interesting range is thus 2*pi > len(x) >= pi. In this range there are 1-2 extermum points. If D > 0 we have two extermum points and the image is [-1,+1]. if D <= 0 there is one extermum and we should check D's factors.
Controlling the rounding error in such a convoluted algorithm is more difficult than straightforward counting.
Interval arithmetic is usually employed when we deal with error ranges instead of pure mathematical real numbers. It's usually assumed that the middle of an interval is near the "actual" value. The error range is then symmetric about this middle point.
It's easy to see that the middle of an interval is usually not its best scalar approximant. We should employ a point that is "median", i.e. there is equal probability for the "actual" value to be below and above it. With such a "median" the interval becomes non-symmetric error range.
A little example will clarify this idea. Suppose we measure some quantity and find that its possible values lay in the range [1,3], with uniform distribution. Suppose further that for some reason we are interested in the square of this quantity. The squared values clearly lay in the range [1,9]. Is the middle point of this range (i.e. 5) a good representative of the result? Clearly not! The median of the first range (i.e. 2) is mapped to 4 in the second. This mapped value is obviously the median of the result and represent it better.
We usually have to assume that our input has a uniform probability inside the error range but this will change with every operation. Triad arithmetic is an extension of interval arithmetic that computes a median in addition to the error range of the result.