250 lines
8.0 KiB
Go
250 lines
8.0 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
|
// SPDX-License-Identifier: MIT
|
|
|
|
package core
|
|
|
|
import (
|
|
"fmt"
|
|
"math"
|
|
"math/big"
|
|
mrand "math/rand/v2"
|
|
"testing"
|
|
)
|
|
|
|
// The compensated-scan evidence file: the sequential float64 CumSum
|
|
// chain against the Neumaier-compensated candidate, both measured
|
|
// against an exact big.Float referent (256 bits of precision, every
|
|
// widening exact). The twins below carry the two walks side by side so
|
|
// one binary answers both the accuracy question and the cost question.
|
|
|
|
// scanLegacyFloats is the sequential chain the Float64 scan arm keeps:
|
|
// out[i] = out[i-1] + src[i], one rounding per step.
|
|
func scanLegacyFloats(src, dst []float64) {
|
|
acc := src[0]
|
|
dst[0] = acc
|
|
for i := 1; i < len(src); i++ {
|
|
acc += src[i]
|
|
dst[i] = acc
|
|
}
|
|
}
|
|
|
|
// scanNeumaierFloats is the compensated candidate: a running high part
|
|
// and a correction, the output the corrected running value. A non-finite
|
|
// partial freezes the correction, so an overflow sticks to infinity and
|
|
// a NaN poisons the tail exactly as the plain chain does.
|
|
func scanNeumaierFloats(src, dst []float64) {
|
|
sum, comp := src[0], 0.0
|
|
dst[0] = sum
|
|
for i := 1; i < len(src); i++ {
|
|
x := src[i]
|
|
t := sum + x
|
|
if !math.IsInf(t, 0) && t == t {
|
|
if math.Abs(sum) >= math.Abs(x) {
|
|
comp += (sum - t) + x
|
|
} else {
|
|
comp += (x - t) + sum
|
|
}
|
|
}
|
|
sum = t
|
|
dst[i] = sum + comp
|
|
}
|
|
}
|
|
|
|
// scanPrefixExact folds the prefixes in big.Float at the given
|
|
// precision; every input is a float64, so the widening is exact.
|
|
func scanPrefixExact(src []float64, prec uint) []*big.Float {
|
|
out := make([]*big.Float, len(src))
|
|
acc := new(big.Float).SetPrec(prec)
|
|
for i, v := range src {
|
|
acc = new(big.Float).SetPrec(prec).Add(acc, new(big.Float).SetPrec(prec).SetFloat64(v))
|
|
out[i] = acc
|
|
}
|
|
return out
|
|
}
|
|
|
|
// scanError reports the largest deviation of the computed prefixes from
|
|
// the exact referent, both absolute and relative to Σ|src|.
|
|
func scanError(got []float64, exact []*big.Float) (abs, rel float64) {
|
|
scale := 0.0
|
|
for _, v := range got {
|
|
scale += math.Abs(v)
|
|
}
|
|
for i := range got {
|
|
ref, _ := exact[i].Float64()
|
|
if e := math.Abs(got[i] - ref); e > abs {
|
|
abs = e
|
|
}
|
|
}
|
|
if scale > 0 {
|
|
rel = abs / scale
|
|
}
|
|
return abs, rel
|
|
}
|
|
|
|
// scanAlternatingData builds the cancellation pattern [1, 1e100, 1,
|
|
// -1e100] repeated: the plain chain loses every small addend it meets,
|
|
// the compensated walk carries them.
|
|
func scanAlternatingData(n int) []float64 {
|
|
src := make([]float64, n)
|
|
for i := 0; i < n; i += 4 {
|
|
src[i], src[i+1], src[i+2], src[i+3] = 1, 1e100, 1, -1e100
|
|
}
|
|
return src
|
|
}
|
|
|
|
// scanSpreadData builds a deterministic log-uniform random sample whose
|
|
// magnitudes span fifty orders, the shape of data that ages a plain
|
|
// chain fastest.
|
|
func scanSpreadData(n int, seed uint64) []float64 {
|
|
rng := mrand.New(mrand.NewPCG(seed, seed))
|
|
src := make([]float64, n)
|
|
for i := range src {
|
|
mag := math.Pow(10, -25+50*rng.Float64())
|
|
if rng.Float64() < 0.5 {
|
|
mag = -mag
|
|
}
|
|
src[i] = mag
|
|
}
|
|
return src
|
|
}
|
|
|
|
// TestScanCompensationAccuracy measures both walks against the exact
|
|
// referent and pins the compensated result's correctness on the
|
|
// cancellation pattern, where the plain chain reports a bare zero.
|
|
func TestScanCompensationAccuracy(t *testing.T) {
|
|
// The referent runs at 512 bits: the alternating pattern's values
|
|
// reach 1e100, whose mantissa needs 333 bits, and every prefix must
|
|
// stay exactly representable for the referent to judge ulp errors.
|
|
const prec = 512
|
|
cases := []struct {
|
|
name string
|
|
src []float64
|
|
}{
|
|
{"alternating 2^20", scanAlternatingData(1 << 20)},
|
|
{"spread 2^20", scanSpreadData(1<<20, 0xC0FFEE)},
|
|
{"spread 2^16", scanSpreadData(1<<16, 7)},
|
|
{"unit 2^20", func() []float64 {
|
|
rng := mrand.New(mrand.NewPCG(42, 42))
|
|
src := make([]float64, 1<<20)
|
|
for i := range src {
|
|
src[i] = rng.Float64()
|
|
}
|
|
return src
|
|
}()},
|
|
}
|
|
for _, tc := range cases {
|
|
exact := scanPrefixExact(tc.src, prec)
|
|
oldGot := make([]float64, len(tc.src))
|
|
newGot := make([]float64, len(tc.src))
|
|
scanLegacyFloats(tc.src, oldGot)
|
|
scanNeumaierFloats(tc.src, newGot)
|
|
oldAbs, oldRel := scanError(oldGot, exact)
|
|
newAbs, newRel := scanError(newGot, exact)
|
|
t.Logf("%s: legacy abs=%.3e rel=%.3e | neumaier abs=%.3e rel=%.3e | improvement %.0fx (abs)",
|
|
tc.name, oldAbs, oldRel, newAbs, newRel, oldAbs/newAbs)
|
|
// A regression guard: the compensated walk is never the worse of
|
|
// the two on these data.
|
|
if newAbs > oldAbs {
|
|
t.Errorf("%s: compensated error %.3e exceeds the chain's %.3e", tc.name, newAbs, oldAbs)
|
|
}
|
|
}
|
|
// The cancellation pin: every complete pattern's prefix returns to
|
|
// 2·k exactly, which the plain chain reports as zero. This pin fails
|
|
// against the uncompenated walk by construction.
|
|
src := []float64{1, 1e100, 1, -1e100, 1, 1e100, 1, -1e100}
|
|
exact := scanPrefixExact(src, prec)
|
|
got := make([]float64, len(src))
|
|
scanNeumaierFloats(src, got)
|
|
for i := range src {
|
|
want, _ := exact[i].Float64()
|
|
if math.Abs(got[i]-want) > 1e-9*math.Max(1, math.Abs(want)) {
|
|
t.Errorf("compensated prefix %d: got %v, want %v", i, got[i], want)
|
|
}
|
|
}
|
|
if last, _ := exact[len(src)-1].Float64(); last != 4 {
|
|
t.Fatalf("referent check: final exact prefix %v, want 4", last)
|
|
}
|
|
// The plain chain's final prefix is 0 here: the document of what the
|
|
// compensation buys.
|
|
legacy := make([]float64, len(src))
|
|
scanLegacyFloats(src, legacy)
|
|
if legacy[len(src)-1] == 4 {
|
|
t.Log("note: the plain chain happened to land exactly on this sample")
|
|
}
|
|
}
|
|
|
|
// TestScanNaNInfSemantics pins the non-finite contract of the
|
|
// compensated walk: an overflow sticks to infinity and a NaN poisons
|
|
// the tail, exactly as the plain chain reports them.
|
|
func TestScanNaNInfSemantics(t *testing.T) {
|
|
cases := [][]float64{
|
|
{math.Inf(1), 1, 2, 3},
|
|
{1, 2, math.Inf(-1), 4},
|
|
{1, math.NaN(), 3, 4},
|
|
{math.MaxFloat64, math.MaxFloat64, 1, 2},
|
|
{1, -1, math.MaxFloat64, math.MaxFloat64},
|
|
}
|
|
for _, src := range cases {
|
|
oldGot := make([]float64, len(src))
|
|
newGot := make([]float64, len(src))
|
|
scanLegacyFloats(src, oldGot)
|
|
scanNeumaierFloats(src, newGot)
|
|
for i := range src {
|
|
wantClass, gotClass := math.Signbit(oldGot[i]), math.Signbit(newGot[i])
|
|
finiteW, finiteG := oldGot[i]-oldGot[i] == 0, newGot[i]-newGot[i] == 0
|
|
if math.IsInf(oldGot[i], 0) != math.IsInf(newGot[i], 0) ||
|
|
math.IsNaN(oldGot[i]) != math.IsNaN(newGot[i]) ||
|
|
(finiteW != finiteG) || (finiteW && wantClass != gotClass) {
|
|
t.Errorf("src %v prefix %d: legacy %v, compensated %v", src, i, oldGot[i], newGot[i])
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// TestCumSumCompensationProduction pins the shipped entry point on the
|
|
// cancellation pattern: the final prefix of two full patterns is 4,
|
|
// which the uncompensated chain reports as 0 (the whole sum cancels and
|
|
// every small addend is lost on the way). This pin fails against the
|
|
// plain-chain walk by construction.
|
|
func TestCumSumCompensationProduction(t *testing.T) {
|
|
a, err := FromFloats([]float64{1, 1e100, 1, -1e100, 1, 1e100, 1, -1e100}, 8)
|
|
if err != nil {
|
|
t.Fatal(err)
|
|
}
|
|
got, err := CumSum(a, 0)
|
|
if err != nil {
|
|
t.Fatal(err)
|
|
}
|
|
if last := got.FloatAt(7); math.Abs(last-4) > 1e-9 {
|
|
t.Errorf("compensated CumSum final prefix: %v, want 4", last)
|
|
}
|
|
// The intermediate prefixes carry their compensation too: after the
|
|
// first pattern the running value is 2, not the chain's 0.
|
|
if mid := got.FloatAt(3); math.Abs(mid-2) > 1e-9 {
|
|
t.Errorf("compensated CumSum prefix 3: %v, want 2", mid)
|
|
}
|
|
}
|
|
|
|
// BenchmarkScanCompensated runs the two walks side by side over one
|
|
// line of the given length; the sub-benchmarks alternate so one binary
|
|
// answers the cost question.
|
|
func BenchmarkScanCompensated(b *testing.B) {
|
|
for _, n := range []int{1 << 10, 1 << 16, 1 << 20} {
|
|
src := scanSpreadData(n, 3)
|
|
b.Run(fmt.Sprintf("legacy/n=%d", n), func(b *testing.B) {
|
|
dst := make([]float64, n)
|
|
for b.Loop() {
|
|
scanLegacyFloats(src, dst)
|
|
}
|
|
b.SetBytes(int64(n) * 8)
|
|
})
|
|
b.Run(fmt.Sprintf("neumaier/n=%d", n), func(b *testing.B) {
|
|
dst := make([]float64, n)
|
|
for b.Loop() {
|
|
scanNeumaierFloats(src, dst)
|
|
}
|
|
b.SetBytes(int64(n) * 8)
|
|
})
|
|
}
|
|
}
|