====== std.sum — kompensiertes Summieren ======

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

→ [[lyx_-_programmiersprache:units:mathematik|Mathematik-Units]] · [[lyx_-_programmiersprache:units:feq|std.feq]] · [[lyx_-_programmiersprache:units:prob|std.prob]] · [[lyx_-_programmiersprache:units:numeric|std.numeric]]

<code lyx>
import std.sum;
</code>

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.

<code lyx>
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;
}
</code>

<code>
Naive    : 0.000000
Kahan    : 0.000000
Neumaier : 1.000000
Pairwise : 0.000000
Exact    : 1.000000
</code>

Das richtige Ergebnis ist **1**.

<box tip>
> **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.
</box>

^ 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:

<code lyx>
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
</code>

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:

<code lyx>
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;
}
</code>

<code>
Count : 5
Mean  : 4.000000
VarP  : 2.000000
VarS  : 2.500000
</code>

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.

<box tip>
> **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.
</box>

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

<box gh>
> **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.
</box>

----

===== 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.