# Demonstrate effects of summation order

When we add lots of floating-point values, we can lose precision if we blindly use std::accumulate() or similar, when the running total becomes much larger than the terms we're adding to it.

I wrote a small program to demonstrate those effects using four different strategies, each briefly commented in the code.

#include <algorithm>
#include <chrono>
#include <functional>
#include <iostream>
#include <iomanip>
#include <numeric>
#include <queue>
#include <ranges>
#include <string_view>
#include <vector>

template<typename T>
void demo(const char* type)
{
// We'll add a range of numbers that use the full precision of floating-point type.
// Dividing integers by 3 gives repeating binary fraction.
constexpr auto n = 5'000'000u;
auto affine = [](auto x){return static_cast<T>(x)/3;};

std::ranges::sized_range auto const numbers
= std::views::iota(0u, n) | std::views::transform(affine);
auto const expected = affine(0.5 * n * (n - 1));

// Here's the ways in which we'll add all the numbers
auto const methods = std::vector<std::pair<std::string_view, std::function<T()>>>{
{"forward accumulate",
[&numbers](){
return std::ranges::fold_left(numbers, T{}, std::plus<T>{});
}},
{"reverse accumulate",
[&numbers](){
// Add values largest first (we expect this to be less
// accurate than forward accumulation).
return std::ranges::fold_right(numbers, T{}, std::plus<T>{});
}},
{"tree pairwise",
[&numbers](){
// Recursively split the input in half, so both sides of
// the addition contain the same number of terms, ±1.
auto const sum = [](this auto& sum, auto a, auto z) {
if (a == z) return T{};
if (a+1 == z) return *a;
auto m = a + (z - a) / 2;
return sum(a, m) + sum(m, z);
};
return sum(numbers.begin(), numbers.end());
}},
[&numbers](){
// Treat input as a set; replace smallest available two
// numbers with their sum, until there's nothing left to
if (std::ranges::empty(numbers)) { return T{}; }
std::priority_queue q(std::greater<T>{}, numbers | std::ranges::to<std::vector>());
for (;;) {
auto a = q.top();
q.pop();
if (q.empty()) { return a; }
auto b = q.top();
q.pop();
q.push(a + b);
}
}},
};

auto const namelen
= std::ranges::max(methods | std::views::transform([](auto& x){ return x.first.size(); }));

for (auto [name, test]: methods) {
std::cout << std::setw(static_cast<int>(namelen)) << name << '<' << type << ">: "
<< std::flush;
auto start = clock::now();
auto result = test();
auto duration = clock::now() - start;
if (result == expected) {
std::cout << "(exact),       ";
} else {
std::cout << std::showpos << std::fixed << std::setprecision(9)
<< (100 * (result - expected) / expected) <<  "%, "
<< std::noshowpos;
}
std::cout << duration_cast<std::chrono::milliseconds>(duration) << '\n' << std::flush;
}
std::cout << '\n';
}

#define demo(type) demo<type>(#type)
int main()
{
demo(float);
demo(double);
demo(long double);
}


When I ran this (compiled by GCC 14 -O3 targeting Intel i7-3770), I got this output:

forward accumulate<float>: +0.156298667%, 4ms
reverse accumulate<float>: -0.306790322%, 3ms
tree pairwise<float>: +0.000006291%, 16ms

forward accumulate<double>: (exact),       18ms
reverse accumulate<double>: (exact),       18ms
tree pairwise<double>: -0.000000000%, 19ms

forward accumulate<long double>: (exact),       24ms
reverse accumulate<long double>: (exact),       23ms
tree pairwise<long double>: +0.000000000%, 25ms


Several observations arise:

• Accuracy is a real problem for float type - all methods give acceptable results for this test using double or wider.
• As expected, largest-first accumulation performs much worse than smallest-first.
• The divide-and-conquer method ("tree pairwise") doesn't optimise as well on float, but does improve accuracy immensely.
• Trying to be clever with the ordering slightly helps accuracy, but costs a lot of time with the heap manipulation, so is rarely likely to be worth doing.

Note that the methods here are demonstration quality and only intended to operate on the sequence provided. In particular, the adaptive pairwise version only works when all inputs are positive. I've made no attempt to generalise this one, because I've already decided it's not worth developing any further.

I'm most interested in reviews of the test methodology and the presentation of the results.

I would be interested in distributions that better expose inaccuracy in double; I was able to get a small error (+0.000000001%) from reverse accumulate by changing numbers to be a series of squares (still divided by three to occupy the whole mantissa):

    std::ranges::sized_range auto const numbers
= std::views::iota(0u, n)
| std::views::transform([](auto x){ return 1. * x * x;})
| std::views::transform(affine);
auto const expected = affine(1./6 * n * (n - 1) * (n + n - 1));


Recommendations arising from the results are also welcome, even though that's not strictly review.

• "presentation of the results" When sweating the details on FP, often we care about 1 ULP, and decimal output can be hard to visually parse, e.g. 0.1 and 0.3 are slightly off, in different directions, but they roundtrip cleanly. So I encourage also displaying floating results in hex form.
– J_H
Feb 4 at 16:02
• Well, I don't even display the result, only the relative error. It was intended more as something I could quote when reviewing poor uses of std::accumulate(), rather than attempting bit-accurate results. Feb 4 at 16:08
• Related: Simd matmul program gives different numerical results - SIMD vectorization (and/or multiple accumulators to hide FP latency) is a step in the direction of pairwise summation. And see Efficient stable sum of a sorted array in AVX2 - Kahan summation (using two floats or two doubles) can be fairly fast and give twice the precision. Of course for float, summing to double might be faster than Kahan summation, even if you're vectorizing so that means half the elements per vector. Feb 5 at 8:10

# Tree pairwise can be written in a more functional style

The "tree pairwise" function can be rewritten to use ranges and views to get the same result, avoiding the need to deal with two iterators:

{"tree pairwise", [&numbers](){
auto const sum = [](this auto& sum, auto range) -> T {
if (range.size() == 1)
return range.front();

return std::ranges::fold_left(
range
| std::views::chunk((range.size() + 1) / 2)
| std::views::transform(sum),
T{}, std::plus{});
};
return sum(numbers);
}},


Is it an improvement? Maybe not. It would be nicer if we had a std::ranges::accumulate() and maybe something named std::ranges::split_n() to split a range in a given number of subranges, then you could write:

return std::ranges::accumulate(range
| std::views::split_n(2)
| std::views::transform(sum));


Of course you could make you own implementations of those.

# Use std::print()

Formatting the output can be greatly simplified by using std::print(), or if you cannot use it yet, at least make use of std::format().

# Unnecessary use of macros

There is no need to use macros. Of course you could simply write the type name manually:

demo<float>("float");


But to avoid the duplication, I see two other ways to get a string containing the name of T inside demo():

1. Use typeid(T).name(). Depending on the platform, this might return the expected name, or a mangled version. In the latter case you could demangle it, but that would not be platform-independent.

2. Use std::source_location::current().function_name(). When called from inside a template function, it should include the actual types of the template parameters. Of course this might contain more information than you need, and the exact format is implementation-defined, but you could do something like:

std::println("{}:", std::source_location::current().function_name());
for (auto [name, test]: methods) {
…
std::println("- {:{}}: ", name, name_len);
…
}


I found it amusing to sweep through multiple affine transforms. With denominators of 11 and 17 we see identical error in both directions!

The random permutation is especially devastating. It offers multiple opportunities for large- and small- magnitude numbers to clump together inducing noticeable error. And then those errors accumulate. We don't really see cancelling errors of opposite sign leading to "mean error is zero". Though with this seed and denominator 7, it did happen to work out.

asc    3   0x0.0p+0                 0.0
dec    3   0x0.0p+0                 0.0
rnd    3   -0x1.8400000000000p-3    -0.189453125

asc    5   0x0.0p+0                 0.0
dec    5   -0x1.0000000000000p-11   -0.00048828125
rnd    5   -0x1.a800000000000p-4    -0.103515625

asc    7   0x0.0p+0                 0.0
dec    7   -0x1.0000000000000p-12   -0.000244140625
rnd    7   0x0.0p+0                 0.0

asc   11   0x1.0000000000000p-12    0.000244140625
dec   11   0x1.0000000000000p-12    0.000244140625
rnd   11   -0x1.1f00000000000p-4    -0.070068359375

asc   13   0x0.0p+0                 0.0
dec   13   0x1.0000000000000p-13    0.0001220703125
rnd   13   -0x1.9c00000000000p-6    -0.025146484375

asc   17   0x1.0000000000000p-13    0.0001220703125
dec   17   0x1.0000000000000p-13    0.0001220703125
rnd   17   0x1.b800000000000p-8     0.0067138671875

asc   19   0x0.0p+0                 0.0
dec   19   0x0.0p+0                 0.0
rnd   19   0x1.5100000000000p-5     0.0411376953125

asc   23   0x0.0p+0                 0.0
dec   23   0x1.0000000000000p-14    6.103515625e-05
rnd   23   0x1.0000000000000p-10    0.0009765625

asc   29   -0x1.0000000000000p-14   -6.103515625e-05
dec   29   -0x1.0000000000000p-14   -6.103515625e-05
rnd   29   -0x1.9a00000000000p-6    -0.0250244140625

asc   31   0x0.0p+0                 0.0
dec   31   -0x1.0000000000000p-14   -6.103515625e-05
rnd   31   -0x1.3400000000000p-7    -0.0093994140625

asc   37   0x0.0p+0                 0.0
dec   37   0x0.0p+0                 0.0
rnd   37   -0x1.5d00000000000p-5    -0.0426025390625

asc   41   0x0.0p+0                 0.0
dec   41   0x1.0000000000000p-14    6.103515625e-05
rnd   41   0x1.6000000000000p-10    0.0013427734375

from typing import Generator, Iterable

from numpy.random import permutation, seed

def _get_fractions(k: int, r: Iterable[int]) -> Generator[float, None, None]:
for i in r:
yield i / k

def _display(label: str, k: int, result: float) -> None:
print(f"{label} {k:4}   {float.hex(result):22s}   {result}")

def add_doubles(n: int = 5_000_000) -> None:
seed(42)

# dividing by K ensures we have lots of mantissa bits set
for k in [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]:
expected = n * (n - 1) / 2 / k
print()
_display("asc", k, expected - sum(_get_fractions(k, range(n))))
_display("dec", k, expected - sum(_get_fractions(k, reversed(range(n)))))
_display("rnd", k, expected - sum(permutation(range(n)) / k))

if __name__ == "__main__":