What is the most numerically precise method for dividing sums or differences?
Consider (ab)/(cd) operation, where a , b , c and d are floating-point numbers (namely, double type in C++). Both (ab) and (cd) are ( sum - correction ) pairs, as in Kahan summation algorithm. Briefly, the specific of these ( sum - correction ) pairs is that sum contains a large value relatively to what's in correction . More precisely, correction contains what didn't fit in sum during summation due to numerical limitations (53 bits of mantissa in double type).
What is the numerically most precise way to calculate (ab)/(cd) given the above speciality of the numbers?
Bonus question: it would be better to get the result also as ( sum - correction ), as in Kahan summation algorithm. So to find (ef)=(ab)/(cd) , rather than just e=(ab)/(cd) .
The div2 algorithm of Dekker (1971) is a good approach.
It requires a mul12(p,q) algorithm which can exactly computes a pair u+v = p*q . Dekker uses a method known as Veltkamp splitting, but if you have access to an fma function, then a much simpler method is
u = p*q
v = fma(p,q,-u)
the actual division then looks like (I've had to change some of the signs since Dekker uses additive pairs instead of subtractive):
r = a/c
u,v = mul12(r,c)
s = (a - u - v - b + r*d)/c
The the sum r+s is an accurate approximation to (ab)/(cd) .
UPDATE: The subtraction and addition are assumed to be left-associative, ie
s = ((((a-u)-v)-b)+r*d)/c
This works because if we let rr be the error in the computation of r (ie r + rr = a/c exactly), then since u+v = r*c exactly, we have that rr*c = auv exactly, so therefore (auvb)/c gives a fairly good approximation to the correction term of (ab)/c .
The final r*d arises due to the following:
(a-b)/(c-d) = (a-b)/c * c/(c-d) = (a-b)/c *(1 + d/(c-d))
= [a-b + (a-b)/(c-d) * d]/c
Now r is also a fairly good initial approximation to (ab)/(cd) so we substitute that inside the [...] , so we find that (auv-b+r*d)/c is a good approximation to the correction term of (ab)/(cd)
对于微小的更正，可能会想到
(a - b) / (c - d) = a/b (1 - b/a) / (1 - c/d) ~ a/b (1 - b/a + c/d)
下一篇: 什么是最精确的数字分割方法？