ClickHouse Sums Floats with Kahan Compensation and a 128-Bit Decimal Accumulator

TL;DR: Floating-point addition loses the small bits. Kahan summation keeps them: alongside the running sum there is a compensation that holds exactly what was rounded away, and the next add subtracts it back in. The state is two Float64s instead of one, the merge has to re-derive the compensation from two partial sums (and the source says do not simplify this expression), and the result is 1.0 where the naive sum gives 0.9999999999999999. For Decimal, the problem is different — the accumulator is simply widened to Decimal128.

Every language has a sum. Most of them get it subtly wrong on floats, and they get it wrong silently, which is the worst way to be wrong.

sql
SELECT sum(0.1), sumKahan(0.1) FROM numbers(10);

text
┌───────────sum(0.1)─┬─sumKahan(0.1)─┐
│ 0.9999999999999999 │             1 │
└────────────────────┴───────────────┘

Ten copies of 0.1 — a number that is not exactly representable in binary floating point — should sum to exactly 1.0. The naive sum gives 0.9999999999999999. The difference is 1.1 × 10⁻¹⁶, and it accumulated over ten additions, each of which rounded.

The naive sum is not a bug. It is the fastest possible summation, and for most workloads the error is irrelevant. ClickHouse’s answer is to ship both: sum for speed, sumKahan for when you care.

Why the naive sum loses bits

The mechanism is not exotic. 0.1 in Float64 is 0.1000000000000000055511151231257827. Adding it to itself ten times, each addition rounds the result back to the nearest representable Float64. The rounding error at each step is a fraction of an ulp, and those errors accumulate.

The classical fix, published by William Kahan in 1965, is to not throw the error away. Keep a second number that remembers it.

ClickHouse’s implementation is four lines:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:298
template <typename T>
struct AggregateFunctionSumKahanData
{
    static_assert(is_floating_point<T>,
        "It doesn't make sense to use Kahan Summation algorithm for non floating point types");

    T sum{};
    T compensation{};

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:307
    template <typename Value>
    void ALWAYS_INLINE addImpl(Value value, T & out_sum, T & out_compensation)
    {
        auto compensated_value = static_cast<T>(value) - out_compensation;
        auto new_sum = out_sum + compensated_value;
        out_compensation = (new_sum - out_sum) - compensated_value;
        out_sum = new_sum;
    }

Read it line by line.

compensated_value = value - compensation. The compensation holds the negative of the low-order bits that were rounded away from the previous sum. Subtracting it from the incoming value folds those bits back in before the addition, so the addition sees a more accurate operand.

new_sum = out_sum + compensated_value. The actual addition. The rounding still happens here — it is unavoidable in Float64 — but the operand is now correct.

out_compensation = (new_sum - out_sum) - compensated_value. This is the recovery step, and it is the clever one. By the Sterbenz lemma, (new_sum - out_sum) is exact whenever out_sum and new_sum are within a factor of two of each other (which is the normal case). So the difference new_sum - out_sum is the exact value of compensated_value plus the rounding error. Subtracting compensated_value from it leaves exactly the rounding error — the bits that got lost.

So compensation is not an estimate. It is the precise value of what was thrown away, to within the precision of a Float64.

out_sum = new_sum. The running sum is the ordinary, rounded sum. The extra precision lives entirely in compensation.

The net effect: sum is the correctly-rounded-ish result and compensation is the residue. If you wanted, you could return sum + compensation for a slightly better answer — but the algorithm’s cleverness is that it doesn’t need to: the residue is carried forward into the next addition, so it never gets lost.

The state is two doubles

The cost is visible in the state size. AggregateFunctionSumKahanData<T> holds T sum and T compensation — for Float64, sixteen bytes instead of eight.

That matters because the state is what gets serialized, merged, and stored in an AggregatingMergeTree:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:417
    static constexpr size_t serialized_size_bound = sizeof(T) * 2;

Compare with the ordinary sum, whose bound is a single value:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:269
    static constexpr size_t serialized_size_bound = sizeof(AccumulateResult);

The merge is the hard part

Merging two compensated sums is not “add the sums and add the compensations”. Float addition is not associative, and the whole Kahan machinery depends on the exact sequence of operations. The implementation says so, in capitals:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:401
    void ALWAYS_INLINE mergeImpl(T & to_sum, T & to_compensation, T from_sum, T from_compensation)
    {
        auto raw_sum = to_sum + from_sum;
        auto rhs_compensated = raw_sum - to_sum;
        /// Kahan summation is tricky because it depends on non-associativity of float arithmetic.
        /// Do not simplify this expression if you are not sure.
        auto compensations = ((from_sum - rhs_compensated) + (to_sum - (raw_sum - rhs_compensated))) + compensation + from_compensation;
        to_sum = raw_sum + compensations;
        to_compensation = compensations - (to_sum - raw_sum);
    }

The trick is the same one as in addImpl, applied in reverse. The merge needs to compute the exact value of to_sum + from_sum, and it does so by recovering the rounding error of the raw addition and adding the two compensations to it.

The line

cpp
auto compensations = ((from_sum - rhs_compensated) + (to_sum - (raw_sum - rhs_compensated))) + compensation + from_compensation;

is doing four things at once:

rhs_compensated = raw_sum - to_sum — the (exact, by Sterbenz) part of from_sum that survived the rounding of raw_sum.
from_sum - rhs_compensated — the part of from_sum that was rounded away.
to_sum - (raw_sum - rhs_compensated) — the corresponding lost bits of to_sum.
+ compensation + from_compensation — the two carried residues.

And the last two lines re-derive the new compensation in the same style as addImpl:

cpp
to_sum = raw_sum + compensations;
to_compensation = compensations - (to_sum - raw_sum);

The comment Do not simplify this expression if you are not sure is a warning from someone who already simplified it and got a bug. Every sub-expression here is load-bearing.

The vectorised path unrolls by 4

sumKahan also has a batch path, and it is deliberately less aggressive than the ordinary sum’s:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:321
    /// Vectorized version
    template <typename Value>
    void NO_INLINE addMany(const Value * __restrict ptr, size_t start, size_t end)
    {
        /// Less than in ordinary sum, because the algorithm is more complicated and too large loop unrolling is questionable.
        /// But this is just a guess.
        constexpr size_t unroll_count = 4;
        T partial_sums[unroll_count]{};
        T partial_compensations[unroll_count]{};

Compare with the ordinary sum, which unrolls by sixteen for Float64:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:94
    /// Vectorized version
    MULTITARGET_FUNCTION_X86_V4(
    MULTITARGET_FUNCTION_HEADER(
    template <typename Value>
    void NO_SANITIZE_UNDEFINED NO_INLINE
    ), addManyImpl, MULTITARGET_FUNCTION_BODY((const Value * __restrict ptr, size_t start, size_t end) /// NOLINT
    {
        ptr += start;
        size_t count = end - start;
        const auto * end_ptr = ptr + count;

        if constexpr (is_floating_point<T>)
        {
            /// Compiler cannot unroll this loop, do it manually.
            /// (at least for floats, most likely due to the lack of -fassociative-math)

            /// Something around the number of SSE registers * the number of elements fit in register.
            constexpr size_t unroll_count = 128 / sizeof(T);
            AccumulateResult partial_sums[unroll_count]{};

128 / sizeof(T) — 128 bytes is “something around the number of SSE registers * the number of elements fit in register”. For Float64 that is 16 accumulators; for Float32, 32.

And the reason for the manual unrolling is spelled out:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:107
            /// Compiler cannot unroll this loop, do it manually.
            /// (at least for floats, most likely due to the lack of -fassociative-math)

Without -fassociative-math, the compiler is not allowed to reorder or group float additions, so it cannot vectorise the reduction. Doing it by hand — sixteen independent accumulators, then combining them — is legal and fast, because each accumulator chain is a serial dependency but the sixteen chains are independent.

The comment clang cannot vectorize the loop if accumulator is class member instead of local variable is the other half:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:127
        /// clang cannot vectorize the loop if accumulator is class member instead of local variable.
        AccumulateResult local_sum{};
        while (ptr < end_ptr)
        {
            Impl::add(local_sum, toAccumulateResult(static_cast<T>(*ptr)));
            ++ptr;
        }
        Impl::add(sum, local_sum);

A local local_sum that the compiler can see is a pure local, and then a single Impl::add(sum, local_sum) at the end. If the accumulator were a class member, clang would have to assume it might be aliased or observed elsewhere, and it would refuse to vectorise.

Note that this manual unrolling changes the summation order — sixteen partial sums instead of one chain. For the ordinary sum that is fine (it is approximate anyway). For sumKahan, the unroll count is 4, and each of the four lanes carries its own compensation, so the four lanes are four independent Kahan chains that are merged at the end:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:337
        while (ptr < unrolled_end)
        {
            for (size_t i = 0; i < unroll_count; ++i)
                addImpl(ptr[i], partial_sums[i], partial_compensations[i]);
            ptr += unroll_count;
        }

        for (size_t i = 0; i < unroll_count; ++i)
            mergeImpl(sum, compensation, partial_sums[i], partial_compensations[i]);

The four lanes are merged with the full mergeImpl — the tricky non-associativity-aware one — not with a simple +=. So the vectorisation does not silently break the compensation.

The `x86-64-v4` dispatch

The ordinary sum has a multi-target implementation:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:138
    /// Vectorized version
    template <typename Value>
    void NO_INLINE addMany(const Value * __restrict ptr, size_t start, size_t end)
    {
#if USE_MULTITARGET_CODE
        if (isArchSupported(TargetArch::x86_64_v4))
        {
            addManyImpl_x86_64_v4(ptr, start, end);
            return;
        }
#endif

        addManyImpl(ptr, start, end);
    }

MULTITARGET_FUNCTION_X86_V4 compiles the same body twice — once for the baseline ISA and once with AVX-512 — and picks at runtime. x86_64_v4 is ClickHouse’s name for the AVX-512 microarchitecture level.

The Decimal path: widen, don’t compensate

For Decimal, the problem is completely different, and the solution is completely different. A Decimal64 is a 64-bit integer with a scale. Adding two of them cannot lose precision — the only way to lose it is to overflow. So the fix is not compensation but a wider accumulator:

cpp
// src/AggregateFunctions/AggregateFunctionSum.cpp:20
template <typename T>
struct SumSimple
{
    /// @note It uses slow Decimal128/256 (cause we need such a variant). sumWithOverflow is faster for Decimal32/64
    using ResultType = std::conditional_t<is_decimal<T>,
                                        std::conditional_t<std::is_same_v<T, Decimal256>, Decimal256, Decimal128>,
                                        NearestFieldType<T>>;
    using AggregateDataType = AggregateFunctionSumData<ResultType>;
    using Function = AggregateFunctionSum<T, ResultType, AggregateDataType, AggregateFunctionTypeSum>;
};

So sum(Decimal32) accumulates in a Decimal128, and sum(Decimal256) accumulates in a Decimal256. The comment slow Decimal128/256 (cause we need such a variant) is the price: a 128-bit addition is more expensive than a 64-bit one, but it never overflows for any realistic row count.

The same widening happens for plain integers and for Float32, via NearestFieldType:

cpp
// src/Core/Field.h:205
template <> struct NearestFieldTypeImpl<UInt16> { using Type = UInt64; };
template <> struct NearestFieldTypeImpl<UInt32> { using Type = UInt64; };

cpp
// src/Core/Field.h:240
template <> struct NearestFieldTypeImpl<Float32> { using Type = Float64; };
template <> struct NearestFieldTypeImpl<Float64> { using Type = Float64; };

sum(UInt16) accumulates in a UInt64, so 65535 + 1 does not wrap. And sum(Float32) accumulates in a Float64 — a Float32 has 24 bits of mantissa, a Float64 has 53, and a sum of many Float32s can exceed the Float32 range even when every individual value is small.

sumWithOverflow is the alternative: it keeps the input type exactly, and the overflow happens for real:

cpp
// src/AggregateFunctions/AggregateFunctionSum.cpp:31
template <typename T>
struct SumSameType
{
    using ResultType = T;
    using AggregateDataType = AggregateFunctionSumData<ResultType>;
    using Function = AggregateFunctionSum<T, ResultType, AggregateDataType, AggregateFunctionTypeSumWithOverflow>;
};

The difference is visible in the result type, not just the value:

sql
CREATE TABLE employees
(
    id UInt32,
    name String,
    monthly_salary UInt16 -- selected so that the sum of values produces an overflow
)
ENGINE = Memory;

INSERT INTO employees VALUES
    (1, 'John', 20000),
    (2, 'Jane', 18000),
    (3, 'Bob', 12000),
    (4, 'Alice', 10000),
    (5, 'Charlie', 8000);

SELECT
    sum(monthly_salary) AS no_overflow,
    sumWithOverflow(monthly_salary) AS overflow,
    toTypeName(no_overflow),
    toTypeName(overflow)
FROM employees;

text
┌─no_overflow─┬─overflow─┬─toTypeName(no_overflow)─┬─toTypeName(overflow)─┐
│       68000 │     2464 │ UInt64                  │ UInt16               │
└─────────────┴──────────┴──────────────────────────┴──────────────────────┘

sum widens UInt16 to UInt64 (via NearestFieldType) and gives the right answer, 68000. sumWithOverflow keeps the type at UInt16, and 68000 mod 65536 = 2464. The ordinary sum is the one you want; the WithOverflow variant exists for the case where you need the result to have the same type as the input — a schema-stability requirement, not a numerical one.

And sumKahan does not apply to Decimal at all:

cpp
// src/AggregateFunctions/AggregateFunctionSum.cpp:39
template <typename T>
struct SumKahan
{
    using ResultType = Float64;
    using AggregateDataType = AggregateFunctionSumKahanData<ResultType>;
    using Function = AggregateFunctionSum<T, ResultType, AggregateDataType, AggregateFunctionTypeSumKahan>;
};

cpp
// src/AggregateFunctions/AggregateFunctionSum.cpp:49
template <typename T> using AggregateFunctionSumKahan =
    std::conditional_t<is_decimal<T>, typename SumSimple<T>::Function, typename SumKahan<T>::Function>;

For a Decimal argument, sumKahan silently becomes the exact sum. The same happens for Interval, which is backed by Int64:

cpp
// src/AggregateFunctions/AggregateFunctionSum.cpp:63
    else if (isInterval(data_type))
    {
        /// `Interval` is backed by `Int64` and is not a numeric type, so it is not part of the
        /// dispatch above. The sum of a set of intervals of one unit is an interval of the same
        /// unit, so only the underlying number of units is summed and the unit of the argument is
        /// kept: values of different units are already brought to their common unit, where one
        /// exists, before the aggregation.
        /// Kahan compensation applies to floating point values only, so, as for `Decimal`, the
        /// exact integer summation is used instead of degrading the count of units to `Float64`.
        using IntervalFunction = std::conditional_t<
            std::is_same_v<typename Function<Int64>::ResultType, Int64>,
            Function<Int64>,
            typename SumSimple<Int64>::Function>;
        res = std::make_shared<IntervalFunction>(argument_types, data_type);
    }

The last sentence of that comment is important: the alternative would have been to sum the units in Float64 (because Function<Int64>::ResultType for Kahan is Float64), which would have introduced exactly the error the user was trying to avoid. Instead, the exact integer summation is used.

The same reasoning appears in the user-facing documentation:

cpp
// src/AggregateFunctions/AggregateFunctionSum.cpp:212
    FunctionDocumentation::Description description_kahan = R"(
Calculates the sum of the numbers with [Kahan compensated summation algorithm](https://en.wikipedia.org/wiki/Kahan_summation_algorithm).
Slower than [`sum`](/reference/functions/aggregate-functions/sum) function.
The compensation works only for [Float](/reference/data-types/float) types.
For the [Decimal](/reference/data-types/decimal) and [Interval](/reference/data-types/special-data-types/interval) types, which are backed by integers, the exact summation of `sum` is used instead.
    )";

The masked trick for conditional float sums

There is one more piece of cleverness in the ordinary sum, and it is a genuinely neat bit of bit-twiddling. When the aggregate has a condition (sumIf) or a null map, the values that should not be added must be zeroed. For integers, multiplication by 0 or 1 does it:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:164
        if constexpr ((is_integer<T> || is_decimal<T>) && !is_over_big_int<T>)
        {
            /// For integers we can vectorize the operation if we replace the null check using a multiplication (by 0 for null, 1 for not null)
            /// https://quick-bench.com/q/MLTnfTvwC2qZFVeWHfOBR3U7a8I
            AccumulateResult local_sum{};
            while (ptr < end_ptr)
            {
                AccumulateResult multiplier = !*condition_map == add_if_zero;
                Impl::add(local_sum, toAccumulateResult(static_cast<T>(*ptr)) * multiplier);
                ++ptr;
                ++condition_map;
            }
            Impl::add(sum, local_sum);
            return;
        }

multiplier is 0 or 1, so the multiply is a conditional select that the compiler can vectorise. Multiplying by 1.0 is exact in IEEE-754, so the “keep” case is unaffected.

For floats you cannot do that: x * 0.0 is 0.0 for finite x but NaN for inf, and x * 1.0 is not the identity for NaN or for signed zero. So the float path reinterprets the value as an integer and masks its bits:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:195
        else if constexpr (is_floating_point<T> && (sizeof(Value) == 4 || sizeof(Value) == 8))
        {
            /// For floating point we use a similar trick as above, except that now we reinterpret the floating point number as an unsigned
            /// integer of the same size and use a mask instead (0 to discard, 0xFF..FF to keep)
            using EquivalentInteger = typename std::conditional_t<sizeof(Value) == 4, UInt32, UInt64>;

            constexpr size_t unroll_count = 128 / sizeof(T);
            AccumulateResult partial_sums[unroll_count]{};

            const auto * unrolled_end = ptr + (count / unroll_count * unroll_count);

            while (ptr < unrolled_end)
            {
                for (size_t i = 0; i < unroll_count; ++i)
                {
                    EquivalentInteger value;
                    memcpy(&value, &ptr[i], sizeof(Value));
                    value &= (!condition_map[i] != add_if_zero) - 1;
                    Value d;
                    memcpy(&d, &value, sizeof(Value));
                    Impl::add(partial_sums[i], toAccumulateResult(static_cast<T>(d)));
                }
                ptr += unroll_count;
                condition_map += unroll_count;
            }

The mask is (!condition_map[i] != add_if_zero) - 1. In C++, true - 1 is 0 and false - 1 is -1, which as an unsigned integer is 0xFFFFFFFFFFFFFFFF. So the mask is all-zeros when the value must be discarded and all-ones when it must be kept. Anding the bit pattern with zero gives +0.0 (all bits zero), which adds nothing to the sum.

Bit-masking the encoding is exact, unlike multiplying. NaN masked to zero becomes +0.0; inf masked to zero becomes +0.0; -0.0 kept as -0.0 stays -0.0. Nothing is corrupted, because the mask is applied to the bits before the value is ever interpreted as a float.

Where the state is stored, and what is serialised

The sum state is written out with a fixed size bound, which lets the framework write states back-to-back into a buffer:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:271
    /// `out` is either a WriteBuffer or a raw `char *` cursor; both are advanced past the state.
    template <typename Out>
    void write(Out & out) const
    {
        writeBinaryLittleEndian(sum, out);
    }

    void read(ReadBuffer & buf)
    {
        readBinaryLittleEndian(sum, buf);
    }

and the Kahan variant writes both halves:

cpp
// src/AggregateFunctions/AggregateFunctionSum.h:419
    /// `out` is either a WriteBuffer or a raw `char *` cursor; both are advanced past the state.
    template <typename Out>
    void write(Out & out) const
    {
        writeBinary(sum, out);
        writeBinary(compensation, out);
    }

    void read(ReadBuffer & buf)
    {
        readBinary(sum, buf);
        readBinary(compensation, buf);
    }

Note the difference in the two calls. The ordinary sum’s accumulator can be an integer — UInt64, Int64, a Decimal128 — and integers have an endianness, so it is written with writeBinaryLittleEndian, which byte-swaps if the machine is big-endian:

cpp
// src/Common/transformEndianness.h:34
template <std::endian ToEndian, std::endian FromEndian = std::endian::native, typename T>
requires is_decimal<T>
inline void transformEndianness(T & x)
{
    transformEndianness<ToEndian, FromEndian>(x.value);
}

The Kahan accumulator is always a Float64 — the static_assert at the top of AggregateFunctionSumKahanData enforces it — so it is written with plain writeBinary, in the machine’s native order.

The trade-offs

What this buys you.

An exact answer where one is possible. sum(Decimal64) accumulates in Decimal128 and cannot lose a single digit. sumKahan(Float64) loses at most the rounding of the final addition.
A guarantee you can state. “Every reported count is within error of the truth” is a contract; “the sum is probably about right” is not.
Speed where precision does not matter. sum is the vectorised, multi-accumulator, x86-64-v4-dispatched version, and it stays that way. The user chooses.
A merge that respects the algorithm. The Kahan merge is not a naive += of the partial sums; it re-derives the compensation, so a distributed sumKahan gives the same answer as a single-threaded one.

What it costs you.

Twice the state. Sixteen bytes instead of eight for a Float64 sum, and twice as many bytes on the wire and in an AggregatingMergeTree column.
A serial dependency. Kahan’s addImpl cannot be reassociated — the compensation depends on the previous sum. The unroll-by-4 mitigates this with four independent chains, but a single chain would still be latency-bound. The ordinary sum with sixteen accumulators has no such constraint.
The merge is genuinely hard. The mergeImpl expression is opaque, order-sensitive, and carries a warning not to touch it. Every future change to the algorithm has to keep the merge consistent with the add.
Silent fallbacks. sumKahan on a Decimal is quietly the exact sum, not a Kahan sum. That is the right behaviour, but it means the function name is a promise about intent, not about the arithmetic that ran.
Masking is not free. The conditional float path does a memcpy in and out of an integer for every element. It is still faster than a branch, but it is not free.

The essential idea is that precision is a state you can carry. The naive sum throws the rounding error away at every step; Kahan keeps it in a second register and folds it into the next addition. The cost is one extra Float64 per group and a merge routine that has to be exactly as careful as the add.