From 045ef28d24d6a7168127163cd7e32ffe8354e90b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Petr=20Balv=C3=ADn?= <opensource@petrbalvin.org>
Date: Mon, 28 Sep 2026 21:36:09 +0200
Subject: [PATCH] fix(signal): refuse filter designs whose coefficients
 overflow

Assisted-by: GLM 5.3 Flash
---
 CHANGELOG.md                |  3 +++
 docs/API.md                 |  4 +++-
 signal/filter.go            | 26 +++++++++++++++++++++
 signal/filterdesign.go      | 13 +++++++++++
 signal/filterdesign_test.go | 46 +++++++++++++++++++++++++++++++++++++
 5 files changed, 91 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index f80784f..1f96ff9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -54,6 +54,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `WriteSVG` keeps extreme but finite data and axis ranges drawable:
   the padding, projection and tick arithmetic fall back to forms whose
   terms stay in range, so the file never carries a NaN coordinate.
+- Filter design refuses an order whose coefficient arithmetic
+  overflows the float64 range instead of shipping a numerator of
+  zeros or NaN.
 
 ## [1.0.0] - 2026-09-03
 
diff --git a/docs/API.md b/docs/API.md
index 40ec571..e6503f6 100644
--- a/docs/API.md
+++ b/docs/API.md
@@ -1054,7 +1054,9 @@ Each design returns the direct-form coefficients `b` (numerator) and `a`
 band edges prewarped to the bilinear axis and the answer exact at the mapped
 frequencies. All of them refuse an order below 1, a non-positive or infinite `fs`,
 and an edge outside `(0, fs/2)`; the band forms additionally require
-`0 < edge1 < edge2 < fs/2`.
+`0 < edge1 < edge2 < fs/2`. An order whose coefficient arithmetic
+overflows the float64 range is refused as well, never returned as a
+numerator of zeros or NaN.
 
 | Call | What it does |
 |---|---|
diff --git a/signal/filter.go b/signal/filter.go
index e5cb2f0..3b4cd4d 100644
--- a/signal/filter.go
+++ b/signal/filter.go
@@ -360,6 +360,9 @@ func butterworth(order int, fs, cutoff float64, highpass bool) (b, a []float64,
 			}
 		}
 		k := polyEvalAtMinusOne(a) / math.Pow(2, float64(order))
+		if err := designCoefficientsFinite(name, order, k, b, a); err != nil {
+			return nil, nil, err
+		}
 		for i := range b {
 			b[i] *= k
 		}
@@ -367,12 +370,35 @@ func butterworth(order int, fs, cutoff float64, highpass bool) (b, a []float64,
 	}
 	b = binomialCoeffs(order)
 	k := polyEvalAtOne(a) / math.Pow(2, float64(order))
+	if err := designCoefficientsFinite(name, order, k, b, a); err != nil {
+		return nil, nil, err
+	}
 	for i := range b {
 		b[i] *= k
 	}
 	return b, a, nil
 }
 
+// designCoefficientsFinite refuses a design whose arithmetic left the
+// float64 range. Past an order of about a thousand the gain divides by
+// an infinite 2^order and the binomial numerator overflows with it,
+// answers that would otherwise ship as a numerator of zeros or NaN
+// presented as a filter: the gain must be finite and non-zero, and
+// every coefficient of both polynomials finite.
+func designCoefficientsFinite(name string, order int, gain float64, b, a []float64) error {
+	if math.IsNaN(gain) || math.IsInf(gain, 0) || gain == 0 {
+		return base.Errf("%s: order %d overflows the coefficient arithmetic; use a lower order", name, order)
+	}
+	for _, poly := range [2][]float64{b, a} {
+		for _, v := range poly {
+			if math.IsNaN(v) || math.IsInf(v, 0) {
+				return base.Errf("%s: order %d overflows the coefficient arithmetic; use a lower order", name, order)
+			}
+		}
+	}
+	return nil
+}
+
 // mulPolyReal multiplies two real polynomials in u = z^{-1} (index m
 // is the coefficient of u^m).
 func mulPolyReal(p, q []float64) []float64 {
diff --git a/signal/filterdesign.go b/signal/filterdesign.go
index 01eafd3..e03f1e4 100644
--- a/signal/filterdesign.go
+++ b/signal/filterdesign.go
@@ -286,9 +286,22 @@ func design(sh shape, proto prototype, w1, w2 float64) (b, a []float64, err erro
 	// reference follows from the roots and is no business of the
 	// scaling.
 	scale := proto.gain / cmplx.Abs(hd)
+	// An extreme order leaves the float64 range here as it does in the
+	// Butterworth pair: an infinite or vanished scale, or a coefficient
+	// past the range, is a refusal rather than a filter of zeros or NaN.
+	if math.IsNaN(scale) || math.IsInf(scale, 0) || scale == 0 {
+		return nil, nil, base.Errf("%s: the order overflows the coefficient arithmetic; use a lower order", name)
+	}
 	for i := range b {
 		b[i] *= scale
 	}
+	for _, poly := range [2][]float64{b, a} {
+		for _, v := range poly {
+			if math.IsNaN(v) || math.IsInf(v, 0) {
+				return nil, nil, base.Errf("%s: the order overflows the coefficient arithmetic; use a lower order", name)
+			}
+		}
+	}
 	return b, a, nil
 }
 
diff --git a/signal/filterdesign_test.go b/signal/filterdesign_test.go
index f56e977..da9c7e5 100644
--- a/signal/filterdesign_test.go
+++ b/signal/filterdesign_test.go
@@ -130,6 +130,52 @@ func TestFilterDesignsAgainstReferenceValues(t *testing.T) {
 	}
 }
 
+// TestExtremeFilterOrders pins the designs against the extreme-but-valid
+// order corner: the entry gates ask only for order ≥ 1, and past a few
+// hundred the coefficient arithmetic leaves the float64 range (the
+// Butterworth gain divides by 2^order, the binomial numerator and the
+// assembled polynomials overflow with it). An order the arithmetic
+// cannot carry must come back as an error, never as a numerator of
+// zeros or NaN presented as a filter; an order it does carry must
+// answer with finite coefficients and a response whose numerator did
+// not vanish.
+func TestExtremeFilterOrders(t *testing.T) {
+	const fs = 1000.0
+	designs := []struct {
+		name  string
+		run   func() ([]float64, []float64, error)
+		probe float64
+	}{
+		{"butterworth-low-1024", func() ([]float64, []float64, error) { return ButterworthLowPass(1024, fs, 100) }, 10},
+		{"butterworth-high-1024", func() ([]float64, []float64, error) { return ButterworthHighPass(1024, fs, 100) }, 490},
+		{"butterworth-band-pass-513", func() ([]float64, []float64, error) { return ButterworthBandPass(513, fs, 100, 300) }, 200},
+		{"chebyshev1-low-1024", func() ([]float64, []float64, error) { return ChebyshevLowPass(1024, fs, 100, 1) }, 10},
+		{"chebyshev2-low-1024", func() ([]float64, []float64, error) { return InverseChebyshevLowPass(1024, fs, 100, 40) }, 10},
+		{"cauer-low-512", func() ([]float64, []float64, error) { return CauerLowPass(512, fs, 100, 1, 60) }, 10},
+	}
+	for _, d := range designs {
+		b, a, err := d.run()
+		if err != nil {
+			// A refusal naming the overflow is the honest answer.
+			continue
+		}
+		for _, v := range b {
+			if math.IsNaN(v) || math.IsInf(v, 0) {
+				t.Fatalf("%s: the numerator holds the non-finite coefficient %g", d.name, v)
+			}
+		}
+		for _, v := range a {
+			if math.IsNaN(v) || math.IsInf(v, 0) {
+				t.Fatalf("%s: the denominator holds the non-finite coefficient %g", d.name, v)
+			}
+		}
+		g := designResponse(b, a, d.probe, fs)
+		if math.IsNaN(g) || g == 0 {
+			t.Fatalf("%s: the response at %g Hz is %.12g, the arithmetic lost the design", d.name, d.probe, g)
+		}
+	}
+}
+
 // TestFilterDesignStability checks that every design's poles sit
 // inside the unit circle, the property direct-form filtering lives
 // and dies by.