Computing average without overflows and with half-decent precision

Question

At some point I needed to compute the average of a big collection of integers. I knew the size of the collection prior computation, but a naive average computation was prone to integer overflow while repeatedly adding every elements divided by the size was prone to loss of precision. To mitigate overflow and loss of precision (speed wasn't an issue), I came up with the following algorithm:

#include <cstdint>
#include <iterator>
#include <type_traits>

template<typename RandomAccessIterator>
auto average(RandomAccessIterator first, RandomAccessIterator last)
    -> long double
{
    static_assert(std::is_integral<
        typename std::iterator_traits<RandomAccessIterator>::value_type
    >::value, "average only available for built-in integer types");

    auto size = std::distance(first, last);
    if (size == 0) return 0.0L;

    std::intmax_t accumulator = 0;
    long double res = 0.0L;

    while (first != last) {
        std::intmax_t tmp;
        if (__builtin_add_overflow(*first, accumulator, &tmp)) {
            res += static_cast<long double>(accumulator) / static_cast<long double>(size);
            accumulator = *first;
        } else {
            accumulator = tmp;
        }
        ++first;
    }

    res += static_cast<long double>(accumulator) / static_cast<long double>(size);
    return res;
}

Basically, here is what the program does:

• It initializes the result and an accumulator with 0.
• It accumulates the values from the collection in an std::intmax_t accumulator as long as said accumulator does not overflow.
• When the accumulator is on the verge of overflowing, the accumulator is divided by the size of the collection and added to the result. It then takes the value that would have made it overflow. The loop continues.
• When then is nothing more to treat in the collection, the current value of the accumulator is divided by the size of the collection and added to the result.

I think that the goal of avoiding overflows is reached: the accumulation is guarded against overflow, and the division of every element to add by the size ensures that the result itself can't overflow either.

Regarding the precision, adding integers in an std::intmax_t as long as possible ensures that no precision is lost before the division. In the best case, a single division will occur. Moreover, the algorithm should divide bit integers most of the time, which may help the precision too. I picked long double to have the best precision guaranteed by the standard too.

Now, there are a few known shortcomings, to name a few:

• Returning 0.0L is probably not the best solution when the collection is empty, but I'm not sure that returning NaN would be better either. One way to circumvent that would be not to handle the case, and document it as a precondition.
• __builtin_add_overflow is a GCC builtin, not a standard C++ function. Unfortunately, I don't know any standard equivalent to it, and manually checking for overflows is hard to do correctly.

The algorithm obviously uses O(1) extra memory. Now, did I properly manage to correctly avoid overflow while mitigating precision loss? Is there anything else to be said about this code?

Answer 1

1. First, having a forward-iterator is good enough, don't ask for more. Yes, it might mean you have to go over the range twice (once to measure its size), but that's a minor irritant.
   Just please, verify it is a forward-iterator like this:

   std::forward_iterator_tag = std::iterator_traits<ForwardIt>::iterator_category{};
   

2. Use the right loop for the job. It's more readable, and more concise.

3. Mark any variable you can const to help comprehension and error-checking.

4. Next, you really have to improve your overflow-handling, as res += static_cast<long double>(accumulator) / static_cast<long double>(size); might loose a small amount of precision.
   Stay in the integer-realm by adding / subtracting enough multiples of size to change the sign of accumulator.

5. Don't cast to long double unless you can guarantee perfect accuracy.

6. You can easily supplant the builtin with a bit of logic:
   Test the sign of accumulate and depending on the result check whether MAXVALUE - accumulate > *first or MINVALUE - accumulate < *first.
   Not sure whether the compiler can optimize that as well as the builtin though...

7. Consider whether making it SFINAE-friendly by avoiding hard errors might be worth the more verbose error-messages.

Modified code:

#include <cstdint>
#include <iterator>
#include <type_traits>
#include <limits>

template<typename ForwardIt>
long double average(ForwardIt first, ForwardIt last)
{
    static_assert(std::is_integral<typename std::iterator_traits<ForwardIt>::value_type>(),
        "average only available for built-in integer types");
    std::forward_iterator_tag = std::iterator_traits<ForwardIt>::iterator_category();

    const auto size = std::distance(first, last);
    if (!size) return 0;

    std::intmax_t accumulator = 0;
    long double res = 0;
    do {
        using snl = std::numeric_limits<decltype(accumulator)>;
        if (accumulator < 0 && snl::min() - accumulator > *first) {
            res += accumulator / size - 1;
            accumulator = accumulator % size + size;
        } else if (accumulator > 0 && snl::max() - accumulator < *first) {
            res += accumulator / size + 1;
            accumulator = accumulator % size - size;
        }
        accumulator += *first;
    } while (++first != last);
    res += accumulator / size;
    res += static_cast<decltype(res)>(accumulator % size) / size;
    return res;
}