Summation algorithm

A summation algorithm is an algorithm that computes the sum of a finite list of numbers ${\textstyle \sum L[i]}$. It is especially relevant in floating-point arithmetic where the associative property ${\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..^[1]

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.^{[2]: 102 [3]}

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.^[4]: §2.1 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 ${\displaystyle a+b}$ versus its floating-point version ${\displaystyle a\oplus b}$, an algorithm called 2Sum gives both ${\displaystyle s:=a\oplus b}$ and the exact error ${\displaystyle t=a+b-(a\oplus b)}$. 2Sum forms the basis of many accurate summation algorithms.^[2]

The intrinsic sensitivity of asummation problem to errors, regardless of how it is computed, is called its condition number and is defined as ${\textstyle \kappa (L)=\sum |L_{i}|/|\sum |L_{i}|}$. The inherent error from machine precision is notated as ${\displaystyle \varepsilon }$.^[5]

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 ${\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 ${\displaystyle O\left({\sqrt {n}}\,\varepsilon \kappa (L)\right)}$.^[6]

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.^[7]

On the other extreme of speed-accuracy tradeoff, the correctly-rounded ${\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.^[8]
• Kirchner and Kulisch's superaccumulator method, which only uses integer arithmetic.^[9] a hardware implementation was described by Müller, Rüb and Rülling.^[10]
  • Zhu and Hayes's HybridSum, which splits into different accumulators by the floating-point exponent field. Uses iFastSum (below) for the final sum.^[11]
  • Neal's xsum, two superaccumulator methods optimized for 64-bit processors.^[12] 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.^[13]
• Zhu and Hayes's iFastSum, a "distillation" method that performs multiple passes of compensated summation.^[11]

One idea for reducing the occurrence of roundoff errors is by making sure most additions happen between values of similar magnitudes. One way of doing so is rearranging the addition into a binary tree of operations, known as pairwise summation. This method takes no additional addition operations compared to the naive method. The worst-case error is ${\displaystyle O(\log(n)\,\varepsilon \kappa (L))}$ and the average-case error is ${\displaystyle O({\sqrt {\log n}}\,\varepsilon \kappa (L))}$.^[14]

Pairwise summation can run in parallel to a factor of ${\displaystyle \log n}$ given light modification.^[15]

Another way to reducing the occurrence of roundoff errors is to make use of the aforementioned 2Sum routine, a compensated summation:^[16]

function KahanBabushkaNeumaierSum2(input)
    var sum = 0.0
    var c = 0.0
    for i = 1 to input.length do
        var y = input[i] + c
        (sum,c) = 2Sum(sum,y)
    next i

    return sum

The typical error bound is ${\displaystyle O(\varepsilon \kappa (L))}$.^[6]

Compensated summation can be made parallel by using a conversion similar to the change from naive to block summation: split the list into several (m) pieces, compute a small number of sub-sums using the compensated method in parallel, then sum them up together using the regular compensated method.^[17]

Depending on the nature of the data, a mixture of summation methods can be used. For example, the FABsum ("fast and accurate block summation") uses a block structure, also dividing te input into m pieces. The m sub-sums are to be generated using a fast method such as the naive one. The sub-sums are then added together using a more accurate method. The error bound is ${\displaystyle (b+1)u+O(u^{2})}$, where b is the block size (i.e. ${\displaystyle m=\left\lceil n/b\right\rceil }$) and u is the unit roundoff (${\displaystyle O(\varepsilon \kappa (L))}$). As a result, a reasonable block size such as 128 or 256 can end up doing the great majority of the summations in a "fast" way, while achieving a constant error magnitude comparable to compensated summation.^[18]

Python includes the precise summation routine as math.fsum.^[19] A few ports are available.^[20] ECMAScript 2026 includes the analogous Math.sumPrecise().^[12]

A compensated summation routine (Neumaier) is imcluded in CPython as sum since version 3.12 and higher.^[21]

In the Julia language, the default implementation of the sum function does pairwise summation for high accuracy with good performance.^[22]

• Floating-point Summation, Dr. Dobb's Journal September, 1996