Problem: Binary floating-point arithmetic rounds intermediate results to available precision, so a small term can disappear beside a much larger total.[3] Our run shows the effect directly: adding 1.0 to 1e16 leaves the stored value unchanged.
Baseline: The plain loop and both compensated functions use O(n) additions and constant scalar accumulator state. The input list and exact reference calculation are outside that space claim.
Solution: The Kahan implementation keeps a correction. The Neumaier implementation chooses its residual calculation from the relative magnitudes of the old total and new value. Both are compared with sum and math.fsum.
Measured result: With 10,000 repetitions of [1e16, 1.0, -1e16], the plain loop and Kahan return 0; Neumaier, sum and math.fsum return 10,000. Reversing the order makes the first two return 1, while the others still return 10,000.
Reproduce the run
Complete code
from __future__ import annotationsimport jsonimport mathimport platformimport statisticsimport sysimport timeitfrom fractions import Fractionfrom pathlib import Pathfrom typing import CallableFLOATS: dict[str, Callable[[list[float]], float]] = {}def naive_left_to_right(values: list[float]) -> float: total = 0.0 for value in values: total += value return totaldef kahan_sum(values: list[float]) -> float: total = 0.0 correction = 0.0 for value in values: adjusted = value - correction next_total = total + adjusted correction = (next_total - total) - adjusted total = next_total return totaldef neumaier_sum(values: list[float]) -> float: total = 0.0 correction = 0.0 for value in values: next_total = total + value if abs(total) >= abs(value): correction += (total - next_total) + value else: correction += (value - next_total) + total total = next_total return total + correctionFLOATS = { 'naive_loop': naive_left_to_right, 'kahan': kahan_sum, 'neumaier': neumaier_sum, 'builtin_sum': sum, 'math_fsum': math.fsum,}def exact_float_sum(values: list[float]) -> float: exact = sum((Fraction.from_float(value) for value in values), Fraction()) return float(exact)def make_values(groups: int, order: tuple[float, float, float]) -> list[float]: return list(order) * groupsdef evaluate_case(name: str, values: list[float]) -> dict[str, object]: expected = exact_float_sum(values) results: dict[str, dict[str, float]] = {} for method, function in FLOATS.items(): actual = function(values) results[method] = { 'value': actual, 'absolute_error': abs(actual - expected), } if results['math_fsum']['absolute_error'] != 0.0: raise AssertionError(f'math.fsum disagreed with exact reference in {name}') if results['neumaier']['absolute_error'] != 0.0: raise AssertionError(f'Neumaier disagreed with exact reference in {name}') return { 'order': name, 'groups': len(values) // 3, 'terms': len(values), 'exact_sum_of_float_inputs': expected, 'results': results, }def benchmark(groups: int) -> dict[str, object]: values = make_values(groups, (1e16, 1.0, -1e16)) exact = exact_float_sum(values) repeats = 7 target_terms_per_sample = 1_500_000 number = max(1, target_terms_per_sample // len(values)) methods: dict[str, dict[str, object]] = {} for method, function in FLOATS.items(): function(values) # one untimed warm-up elapsed = timeit.repeat(lambda fn=function: fn(values), number=number, repeat=repeats) per_call = [sample / number for sample in elapsed] result = function(values) methods[method] = { 'value': result, 'absolute_error': abs(result - exact), 'samples_seconds_per_sum': per_call, 'median_seconds_per_sum': statistics.median(per_call), } return { 'groups': groups, 'terms': len(values), 'order': 'large, small, negative large', 'exact_sum_of_float_inputs': exact, 'repeat_count': repeats, 'calls_per_sample': number, 'target_terms_per_sample': target_terms_per_sample, 'methods': methods, }def main() -> None: cases = [] for groups in (1, 100, 1_000, 10_000): cases.append(evaluate_case('large, small, negative large', make_values(groups, (1e16, 1.0, -1e16)))) cases.append(evaluate_case('large, negative large, small', make_values(groups, (1e16, -1e16, 1.0)))) cases.append(evaluate_case('empty input', [])) for method, function in FLOATS.items(): if function([]) != 0.0: raise AssertionError(f'{method} failed empty-input check') timings = [benchmark(groups) for groups in (100, 1_000, 10_000)] output: dict[str, object] = { 'runtime': { 'python': sys.version.split()[0], 'implementation': platform.python_implementation(), 'platform': platform.platform(), 'machine': platform.machine(), }, 'correctness_cases': cases, 'timings': timings, 'method': { 'timer': 'timeit.repeat, default perf_counter', 'warmups_per_method_and_size': 1, 'timed_repeats': 7, 'summary': 'median of seven per-call samples; every raw sample retained', 'input': 'the displayed three-value order repeated; construction excluded', 'correctness': 'Fraction.from_float exact sum, rounded once to float', 'scope': 'one process, one machine; timings are not portable rankings', }, } path = Path(__file__).with_name('benchmark-output.json') path.write_text(json.dumps(output, indent=2, ensure_ascii=False) + '\n', encoding='utf-8') # Keep stdout readable for the article while retaining the full-precision # sample vectors in benchmark-output.json. lines = [ f"runtime={output['runtime']['implementation']} {output['runtime']['python']}; " f"platform={output['runtime']['platform']}; machine={output['runtime']['machine']}", 'correctness (exact Fraction sum rounded once to float):', ] for case in cases: results = ', '.join( f"{name}={item['value']!r} (abs_error={item['absolute_error']!r})" for name, item in case['results'].items() ) lines.append(f" {case['order']}, groups={case['groups']}, exact={case['exact_sum_of_float_inputs']!r}: {results}") lines.append('timing policy: 1 warm-up, 7 samples, median; timeit default timer; GC disabled by timeit') lines.append('raw samples are microseconds per sum, rounded to 3 decimals for display; full precision retained in JSON') for row in timings: lines.append(f"terms={row['terms']}; order={row['order']}; groups={row['groups']}; calls_per_sample={row['calls_per_sample']}") builtin_median = row['methods']['builtin_sum']['median_seconds_per_sum'] for name, item in row['methods'].items(): samples = ','.join(f"{sample * 1_000_000:.3f}" for sample in item['samples_seconds_per_sum']) median = item['median_seconds_per_sum'] * 1_000_000 ratio = item['median_seconds_per_sum'] / builtin_median lines.append(f" {name}: value={item['value']!r}; abs_error={item['absolute_error']!r}; median_us={median:.3f}; ratio_vs_builtin_sum={ratio:.3f}x; samples_us=[{samples}]") print('\n'.join(lines))if __name__ == '__main__': main()
Command
python3 summation_benchmark.py
Output
$ python3 summation_benchmark.pyruntime=CPython 3.13.5; platform=Linux-6.12.47+rpt-rpi-v8-aarch64-with-glibc2.41; machine=aarch64correctness (exact Fraction sum rounded once to float): large, small, negative large, groups=1, exact=1.0: naive_loop=0.0 (abs_error=1.0), kahan=0.0 (abs_error=1.0), neumaier=1.0 (abs_error=0.0), builtin_sum=1.0 (abs_error=0.0), math_fsum=1.0 (abs_error=0.0) large, negative large, small, groups=1, exact=1.0: naive_loop=1.0 (abs_error=0.0), kahan=1.0 (abs_error=0.0), neumaier=1.0 (abs_error=0.0), builtin_sum=1.0 (abs_error=0.0), math_fsum=1.0 (abs_error=0.0) large, small, negative large, groups=100, exact=100.0: naive_loop=0.0 (abs_error=100.0), kahan=0.0 (abs_error=100.0), neumaier=100.0 (abs_error=0.0), builtin_sum=100.0 (abs_error=0.0), math_fsum=100.0 (abs_error=0.0) large, negative large, small, groups=100, exact=100.0: naive_loop=1.0 (abs_error=99.0), kahan=1.0 (abs_error=99.0), neumaier=100.0 (abs_error=0.0), builtin_sum=100.0 (abs_error=0.0), math_fsum=100.0 (abs_error=0.0) large, small, negative large, groups=1000, exact=1000.0: naive_loop=0.0 (abs_error=1000.0), kahan=0.0 (abs_error=1000.0), neumaier=1000.0 (abs_error=0.0), builtin_sum=1000.0 (abs_error=0.0), math_fsum=1000.0 (abs_error=0.0) large, negative large, small, groups=1000, exact=1000.0: naive_loop=1.0 (abs_error=999.0), kahan=1.0 (abs_error=999.0), neumaier=1000.0 (abs_error=0.0), builtin_sum=1000.0 (abs_error=0.0), math_fsum=1000.0 (abs_error=0.0) large, small, negative large, groups=10000, exact=10000.0: naive_loop=0.0 (abs_error=10000.0), kahan=0.0 (abs_error=10000.0), neumaier=10000.0 (abs_error=0.0), builtin_sum=10000.0 (abs_error=0.0), math_fsum=10000.0 (abs_error=0.0) large, negative large, small, groups=10000, exact=10000.0: naive_loop=1.0 (abs_error=9999.0), kahan=1.0 (abs_error=9999.0), neumaier=10000.0 (abs_error=0.0), builtin_sum=10000.0 (abs_error=0.0), math_fsum=10000.0 (abs_error=0.0) empty input, groups=0, exact=0.0: naive_loop=0.0 (abs_error=0.0), kahan=0.0 (abs_error=0.0), neumaier=0.0 (abs_error=0.0), builtin_sum=0 (abs_error=0.0), math_fsum=0.0 (abs_error=0.0)timing policy: 1 warm-up, 7 samples, median; timeit default timer; GC disabled by timeitraw samples are microseconds per sum, rounded to 3 decimals for display; full precision retained in JSONterms=300; order=large, small, negative large; groups=100; calls_per_sample=5000 naive_loop: value=0.0; abs_error=100.0; median_us=18.035; ratio_vs_builtin_sum=3.572x; samples_us=[18.734,18.201,17.974,18.028,17.869,18.035,18.037] kahan: value=0.0; abs_error=100.0; median_us=50.382; ratio_vs_builtin_sum=9.978x; samples_us=[50.985,51.212,51.626,50.382,50.023,49.763,49.950] neumaier: value=100.0; abs_error=0.0; median_us=79.820; ratio_vs_builtin_sum=15.809x; samples_us=[84.024,79.524,79.481,80.957,79.487,79.993,79.820] builtin_sum: value=100.0; abs_error=0.0; median_us=5.049; ratio_vs_builtin_sum=1.000x; samples_us=[5.053,5.049,5.049,5.049,5.042,5.047,5.048] math_fsum: value=100.0; abs_error=0.0; median_us=5.092; ratio_vs_builtin_sum=1.009x; samples_us=[5.042,5.126,5.092,5.093,5.078,5.111,5.091]terms=3000; order=large, small, negative large; groups=1000; calls_per_sample=500 naive_loop: value=0.0; abs_error=1000.0; median_us=176.420; ratio_vs_builtin_sum=4.071x; samples_us=[176.205,176.622,176.299,176.801,176.394,176.420,177.921] kahan: value=0.0; abs_error=1000.0; median_us=496.948; ratio_vs_builtin_sum=11.468x; samples_us=[496.948,497.076,497.037,496.311,492.248,492.199,500.582] neumaier: value=1000.0; abs_error=0.0; median_us=815.074; ratio_vs_builtin_sum=18.809x; samples_us=[815.074,870.080,809.950,871.858,801.907,868.430,802.418] builtin_sum: value=1000.0; abs_error=0.0; median_us=43.335; ratio_vs_builtin_sum=1.000x; samples_us=[46.466,43.314,43.652,47.437,42.643,42.798,43.335] math_fsum: value=1000.0; abs_error=0.0; median_us=48.109; ratio_vs_builtin_sum=1.110x; samples_us=[48.109,48.345,48.097,48.076,48.507,48.461,47.299]terms=30000; order=large, small, negative large; groups=10000; calls_per_sample=50 naive_loop: value=0.0; abs_error=10000.0; median_us=1753.962; ratio_vs_builtin_sum=3.703x; samples_us=[1751.977,1755.291,1751.795,1757.248,1755.415,1753.319,1753.962] kahan: value=0.0; abs_error=10000.0; median_us=4911.352; ratio_vs_builtin_sum=10.368x; samples_us=[4959.914,5008.631,4907.038,4911.352,4961.688,4901.555,4901.727] neumaier: value=10000.0; abs_error=0.0; median_us=7904.598; ratio_vs_builtin_sum=16.688x; samples_us=[7877.997,7855.274,7904.598,8203.830,8007.593,7913.134,7900.750] builtin_sum: value=10000.0; abs_error=0.0; median_us=473.682; ratio_vs_builtin_sum=1.000x; samples_us=[473.917,474.582,473.632,476.465,473.682,463.110,425.752] math_fsum: value=10000.0; abs_error=0.0; median_us=479.657; ratio_vs_builtin_sum=1.013x; samples_us=[483.076,479.445,479.657,478.880,485.069,480.091,472.152]exit=0
Environment
The run used CPython 3.13.5 on Linux 6.12.47, aarch64. The container did not expose the processor model.
Methodology
The exact reference adds the rational values represented by the input floats with Fraction.from_float, then converts the total to a float once. Correctness checks cover one, 100, 1,000 and 10,000 groups in both term orders, plus empty input.
Timing covers 300, 3,000 and 30,000 terms, using order A: [1e16, 1.0, -1e16]. Each sample makes enough calls to process about 1.5 million terms. The script runs one warm-up and seven timed repeats for each method and size, then reports the median per call while retaining all seven measurements. Input construction and the exact reference are outside the timed region; the Python wrapper that calls each summation function is inside it. timeit disables garbage collection by default.[4]
The displayed command runs the complete benchmark. Its output contains the correctness results and the raw timing vectors; benchmark-output.json retains full-precision values and the environment fields.
Where Kahan loses this correction
For one group, the mathematical total of 1e16 + 1.0 - 1e16 is 1. Kahan’s first two additions leave a running total of 1e16 and a correction of -1. On the final input, its adjusted value is -1e16 + 1; that rounds back to -1e16 in this run. The correction is then zero, and the result is 0.
The Neumaier function uses a different residual calculation. When the existing total is larger in magnitude, it adds (total - next_total) + value to a separate correction. For this sequence, that retains the lost 1; the final return adds it back. The code and output show the result for this exact input, not a guarantee that one compensated method wins for every sequence.
What the built-ins do
Python’s documentation says the accuracy and commutativity of sum() for floats improved in Python 3.12 on most builds.[1] That matches this CPython 3.13 run: sum() returns 10,000 for both tested orders.
math.fsum() tracks multiple intermediate partial sums.[2] The Python tutorial describes it as slower than sum() but more accurate in uncommon cases where large values cancel and leave a small result.[3] Here the seven-sample medians for 30,000 terms were 0.474 ms for sum() and 0.480 ms for math.fsum(). The plain loop took 1.754 ms, Kahan 4.911 ms and Neumaier 7.905 ms. Against sum(), those medians were 3.703×, 10.368× and 16.688× as long for the plain loop, Kahan and Neumaier; math.fsum() was 1.013×. Those timings describe these implementations on this one host.
The math.fsum() documentation notes a narrow platform caveat: on some non-Windows builds using extended-precision addition, an intermediate sum may be double-rounded and differ by one least-significant bit.[2] This run does not test other platforms.
Scope
The input is deliberately hostile to ordinary accumulation. It is a useful check for cancellation, not a sample of typical application data. No measurements here cover NaNs, infinities, overflow, mixed numeric types, decimal quantities or another processor. The Python tutorial recommends decimal when exact decimal representation is required.[3]
Sources
[1] Built-in Functions, Python 3.13 documentation
[2] math.fsum, Python 3.13 documentation
