std.sum — kompensiertes Summieren

Summen, die den Rundungsfehler mitführen statt ihn zu verlieren — und laufende Statistik nach Welford, die ohne zweiten Durchgang auskommt.

→ Mathematik-Units · std.feq · std.prob · std.numeric

import std.sum;

Schicht 0, hängt an std.math und std.alloc. 25 Funktionen, zwei Größenkonstanten (SUM_ACC_BYTES = 16, SUM_STAT_BYTES = 24). Alle Beispiele mit lyxc 1.0.21A übersetzt und ausgeführt.

Das Problem in einer Zeile

Wer 1e16, 1 und −1e16 in dieser Reihenfolge addiert, bekommt null. Die 1 geht beim ersten Schritt verloren, weil 1e16 + 1 in f64 wieder 1e16 ist — und nach dem Abziehen ist sie unwiederbringlich weg.

unit main;
import std.io;
import std.string;
import std.alloc;
import std.sum;

fn ZF(t: pchar, v: f64): void { Print(t); PrintF64(v); }

fn main(): int64 {
    var n: int64 := 3;
    var a: int64 := alloc(n * 8);
    pokef64(a,      10000000000000000.0);
    pokef64(a + 8,  1.0);
    pokef64(a + 16, -10000000000000000.0);

    ZF("Naive    : ", SumNaive(a, n));
    ZF("Kahan    : ", SumKahan(a, n));
    ZF("Neumaier : ", SumNeumaier(a, n));
    ZF("Pairwise : ", SumPairwise(a, n));

    var work: int64 := alloc(64 * 8);
    ZF("Exact    : ", SumExact(a, n, work, 64));
    return 0;
}

Naive    : 0.000000
Kahan    : 0.000000
Neumaier : 1.000000
Pairwise : 0.000000
Exact    : 1.000000

Das richtige Ergebnis ist 1.

Bemerkenswert: Kahan versagt hier ebenfalls. Das ist kein Fehler der Umsetzung, sondern die bekannte Grenze des Verfahrens. Kahan führt einen Korrekturterm mit, verrechnet ihn aber gegen den neuen Summanden — ist dieser selbst riesig (−1e16), geht die gesammelte Korrektur beim Verrechnen unter. Neumaier unterscheidet genau diesen Fall: er prüft, welcher der beiden Werte betragsmäßig größer ist, und legt den Fehler entsprechend an. Deshalb ist Neumaier und nicht Kahan die richtige Vorgabe.

Verfahren	Kosten	Wann
`SumNaive`	1 Addition	nur wenn alle Werte dieselbe Größenordnung haben
`SumPairwise`	1 Addition, rekursiv	gute Genauigkeit fast umsonst, hilft aber nicht gegen Auslöschung
`SumKahan`	~4 Operationen	klassisch, aber blind gegen große Summanden (siehe oben)
`SumNeumaier`	~5 Operationen	die Vorgabe — fängt auch den Fall oben
`SumExact`	Arbeitspuffer	fehlerfrei, Ergebnis ist die korrekt gerundete Summe

SumExact braucht einen Arbeitsbereich: SumExact(arr, n, work, workCap). Er sammelt die nicht-überlappenden Teilsummen nach Shewchuk; workCap begrenzt deren Anzahl.

Zwei Zahlen exakt addieren

Die Grundlage aller Verfahren: die Summe zweier f64-Zahlen lässt sich exakt als Paar aus Ergebnis und Fehler darstellen.

Funktion	Bedingung
`SumTwoSum(a, b, errAddr)`	keine — immer gültig, 6 Operationen
`SumFastTwoSum(a, b, errAddr)`	setzt voraus, dass der Betrag von `a` mindestens so groß ist wie der von `b` — dafür nur 3 Operationen

Beide liefern die Summe zurück und legen den Rest unter errAddr ab. SumFastTwoSum ohne die Vorbedingung liefert stillschweigend Unsinn — wer die Reihenfolge nicht garantieren kann, nimmt SumTwoSum.

Laufend summieren

Wenn die Werte nicht als Feld vorliegen, sondern nach und nach eintreffen, gibt es die Akkumulator-Form. Der Aufrufer stellt SUM_ACC_BYTES Bytes:

var acc: int64 := alloc(SUM_ACC_BYTES);
SumNeumaierInit(acc);
SumNeumaierAdd(acc, wert);        // beliebig oft
var s: f64 := SumNeumaierValue(acc);
var e: f64 := SumNeumaierError(acc);   // die mitgeführte Korrektur

Dieselbe Form gibt es für Kahan (SumKahanInit / Add / Value). SumNeumaierError gibt Auskunft darüber, wie viel kompensiert wurde — ein brauchbares Maß dafür, ob das naive Summieren an dieser Stelle getragen hätte.

Über Felder rechnen außerdem SumDot(a, b, n) (kompensiertes Skalarprodukt) und SumMean(arr, n).

Statistik nach Welford

Mittelwert und Streuung in einem Durchgang, ohne die Werte zu behalten:

unit main;
import std.io;
import std.string;
import std.alloc;
import std.sum;

fn ZF(t: pchar, v: f64): void { Print(t); PrintF64(v); }

fn main(): int64 {
    var m: int64 := 5;
    var b: int64 := alloc(m * 8);
    var i: int64 := 0;
    while (i < m) {
        pokef64(b + i * 8, 2.0 + (i as f64));
        i := i + 1;
    }

    var st: int64 := alloc(SUM_STAT_BYTES);
    SumStatInit(st);
    i := 0;
    while (i < m) {
        var v: f64 := peekf64(b + i * 8);
        SumStatAdd(st, v);
        i := i + 1;
    }

    Print("Count : "); PrintLn(IntToStr(SumStatCount(st)));
    ZF("Mean  : ", SumStatMean(st));
    ZF("VarP  : ", SumStatVarianceP(st));
    ZF("VarS  : ", SumStatVarianceS(st));
    return 0;
}

Count : 5
Mean  : 4.000000
VarP  : 2.000000
VarS  : 2.500000

Für die Werte 2…6: Mittelwert 4, Varianz der Grundgesamtheit 2, Varianz der Stichprobe 2,5. Der Unterschied ist der Nenner — n gegen n−1. Wer eine Stichprobe hat und VarianceP nimmt, unterschätzt die Streuung systematisch.

Warum nicht die Schulformel? E[x²] − E[x]² verlangt zwei große, fast gleiche Zahlen voneinander abzuziehen — bei Werten weit von null verliert das Ergebnis alle Stellen und kann sogar negativ werden. Welford führt Mittelwert und Abweichungsquadrat schrittweise mit und hat das Problem nicht.

SumStatMerge(a, b, out) führt zwei getrennt gerechnete Statistiken zusammen — für parallel ausgewertete Teilmengen, ohne die Werte erneut zu lesen.

Der Akkumulator ist roher Speicher. SumStatInit muss vor dem ersten Add laufen, sonst wird auf uninitialisiertem Speicher gerechnet. Die Größen SUM_ACC_BYTES und SUM_STAT_BYTES sind Konstanten der Unit — nicht selbst ausrechnen, sondern verwenden.

Wo es noch naiv summiert wird

Innerhalb der Mathematik-Units sind std.dist, std.sparse, std.linalg und std.prob auf kompensiertes Summieren umgestellt. Noch nicht umgestellt sind std.tensor (3 Stellen), std.convolve (8), std.wavelet (2) und std.fft (1) — bei convolve und fft betrifft es die inneren Faltungsschleifen. Wer dort über sehr viele Werte mit stark verschiedenen Größenordnungen rechnet, sollte das Ergebnis gegen eine SumExact-Rechnung prüfen.

Letzte Aktualisierung: 2026-08-16 · alle Beispiele mit lyxc 1.0.21A übersetzt und ausgeführt; Summen- und Statistikwerte gegen eine unabhängige Referenzrechnung geprüft.