A summation algorithm is an algorithm that computes the sum of a finite list of numbers ∑ L [ i ] {\textstyle \sum L[i]} . It is especially relevant in floating-point arithmetic where the associative property ( a + b ) + c = a + ( b + c ) {\displaystyle (a+b)+c=a+(b+c)} does not hold like it does in (mathematical) real numbers, rational numbers, fixed-point numbers, and unsigned integers, so that the order of calculation can affect the final result. A closely related problem is the calculation of dot products, both of which have numerous proposed algorithms that differ in terms of simplicity, speed (single-thread and parallelized), and accuracy..
Background Floating-point works similarly to a scientific notation with a limited range on the exponent and a limited number of digits on the significand (left part, mantissa). When numbers of very different magnitudes (absolute values) are added together, the result may be unchanged from the one with the larger magnitude, making any information in the smaller one lost. When numbers of very similar magnitudes but of opposite signs are added together, cancellation results, with all but the most inaccurate digits being 0 in the result, with potentially catastrophic consequences for precision. More concretely, the summation operator "+" of floating-point arithmetic is correctly-rounded, i.e. it has to match the rounding mode currently set. In the case of the two "to nearest" modes, the error is guaranteed to be ≤ 0.5 ulp, and for the other modes ≤ 1 ulp; in all cases, it is deterministic. A ulp is a unit in the last place: the magnitude difference represented by the last digit of the mantissa changing from a 0 to 1, and its own magnitude is determined by the exponent. The error incurred by floating-point rules, relative to the ideal of infinite-precision arithmetic, is called the roundoff error. For the basic case of calculating a + b {\displaystyle a+b} versus its floating-point version a ⊕ b {\displaystyle a\oplus b} , an algorithm called 2Sum gives both s := a ⊕ b {\displaystyle s:=a\oplus b} and the exact error t = a + b − ( a ⊕ b ) {\displaystyle t=a+b-(a\oplus b)} . 2Sum forms the basis of many accurate summation algorithms. The intrinsic sensitivity of asummation problem to errors, regardless of how it is computed, is called its condition number and is defined as κ ( L ) = ∑ | L i | / | ∑ | L i | {\textstyle \kappa (L)=\sum |L_{i}|/|\sum |L_{i}|} . The inherent error from machine precision is notated as ε {\displaystyle \varepsilon } .
Naive summation Naive summation works simply by going over the numbers one-by-one, calculating
function sumNaive(L) var sum = 0.0 for i = 1 to L.length do sum = sum + L[i] return sum
The worst-case error for a list L of length n is O ( n ε κ ( L ) ) {\displaystyle O(n\,\varepsilon \kappa (L))} , corresponding to a case where the rounding errors all add up in the same direction. The average error, corresponding to a random walk of roundoff errors, is O ( n ε κ ( L ) ) {\displaystyle O\left({\sqrt {n}}\,\varepsilon \kappa (L)\right)} . Naive summation also has a batched version, blocked summation, which entails dividing L into a number of smaller lists, summing each of them naively, then summing the sums naively. This is usually done to allow parallel computation. It also reduces the error growth by a factor of 1/b.
Precise summation On the other extreme of speed-accuracy tradeoff, the correctly-rounded ∑ L [ i ] {\textstyle \sum L[i]} may be computed. Naively this may be done via arbitrary-precision arithmetic, but more efficient methods exist:
Shewchuk's method, an adaptive method making use of Floating-point expansions. This can take anywhere from linear (optimistic, and near typical) to quadratic time complexity and anywhere from constant to linear space complexity. Kirchner and Kulisch's superaccumulator method, which only uses integer arithmetic. a hardware implementation was described by Müller, Rüb and Rülling. Zhu and Hayes's HybridSum, which splits into different accumulators by the floating-point exponent field. Uses iFastSum (below) for the final sum. Neal's xsum, two superaccumulator methods optimized for 64-bit processors. Integer-only for the main phase of accumulation. The small superaccumulator is 67 64-bit integers covering the whole exponent range, with carry-propagation so it can act like a 2176-bit integer. The large superaccumulator, used when the array is very large, collects a partial sum for every possible combination of exponent bits and sign bit to ease the load on the small superaccumulator. Zhu and Hayes's iFastSum, a "distillation" method that performs multiple passes of compensated summation.
… excerpt ends here. Continue reading the full article.
