chore: prepare release v1.0.1

Assisted-by: GLM 5.3 Flash

CHANGELOG.md

@@ -5,6 +5,65 @@ All notable changes to **Tensor** are documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [development]
+
+### Added
+
+-
+
+## [1.0.1] - 2026-09-28
+
+### Fixed
+
+**Linear algebra.**
+
+- `Cholesky`, `Eigen`, `EigenComplex` and the Cholesky rank-one
+  update refuse a matrix or vector holding a NaN or infinity instead
+  of returning an all-NaN factor or spectrum with a nil error, the
+  refusal the sparse solvers already make.
+- `EigenComplex` converges on a Hermitian matrix whose off-diagonal
+  entries all sit just under the per-entry skip level: the sweep no
+  longer exhausts its passes on a matrix it should have declared
+  converged.
+
+**Statistics.**
+
+- `FitHiddenMarkovModel` fits a single-observation sequence instead of
+  crashing the process: a transition row with no evidence keeps its
+  previous estimate rather than dividing by zero.
+- `Viterbi` refuses a sequence of probability zero under the model,
+  the same refusal `Forward` makes, instead of returning a meaningless
+  path beside a log probability of -Inf.
+- The Student-t tails stay accurate past the point t squared
+  overflows float64: `StudentTCDF`, the regression coefficient
+  p-values and `NoncentralTCDF` at extreme t answer the tail the
+  format still holds instead of a silent zero or an error naming a
+  NaN.
+- `LinearRegression` and `WeightedLinearRegression` keep their
+  inference alive when the squared deviations underflow to zero: the
+  standard errors, t and F statistics read factored sums of squares
+  instead of reporting infinity beside zero evidence.
+
+**Integration and optimisation.**
+
+- `TriangleMesh2D.BoundaryEdges` returns the mesh's own edges: the
+  boundary pairs sort as pairs, never as one flat index list that
+  interleaves the endpoints of unrelated edges.
+- `IntegrateFunction` names itself in its errors; the adaptive
+  quadrature no longer reports under a name absent from the surface.
+- `MinimiseDifferentialEvolution` refuses a non-finite bound instead
+  of drawing a NaN population and returning a NaN point with no
+  error.
+
+**Signal and plots.**
+
+- `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
 
 The initial release of Tensor, a scientific computing library in pure

docs/API.md

@@ -629,7 +629,7 @@ neighbours are the way in.
 | `Inv(a)` | returns the inverse of a square matrix, through the same LU kernel. |
 | `Det(a)` | returns the determinant of a square real matrix; a singular matrix gives 0, not an error. |
 | `DetComplex(a)` | the same for a square complex matrix, as a `complex128`. |
-| `Cholesky(a)` | returns the lower factor L with a = L·Lᵀ for a symmetric positive definite a; a non-positive pivot is refused. |
+| `Cholesky(a)` | returns the lower factor L with a = L·Lᵀ for a symmetric positive definite a; a non-positive pivot and a non-finite entry are refused. |
 | `CholeskyUpdate(l, x)` | returns the lower factor of A + x·xᵀ from the factor of A, by orthogonal rotations. |
 | `CholeskyDowndate(l, x)` | returns the lower factor of A − x·xᵀ, by hyperbolic rotations; leaving the positive definite cone is an error. |
 | `QR(a)` | returns q (m×m orthogonal) and r (m×n upper triangular) with a = q·r, for m ≥ n. |
@@ -641,8 +641,8 @@ neighbours are the way in.
 
 | Call | What it does |
 |---|---|
-| `Eigen(a)` | returns the eigenvalues of a real symmetric matrix, ascending, and the orthonormal eigenvectors as columns; asymmetry beyond a scale-relative 1e-12 is refused. |
+| `Eigen(a)` | returns the eigenvalues of a real symmetric matrix, ascending, and the orthonormal eigenvectors as columns; asymmetry beyond a scale-relative 1e-12 and a non-finite entry are refused. |
-| `EigenComplex(a)` | the same for a complex Hermitian matrix; values are real, vectors complex, and the same 1e-12 Hermitian tolerance applies. |
+| `EigenComplex(a)` | the same for a complex Hermitian matrix; values are real, vectors complex, the same 1e-12 Hermitian tolerance applies, and a non-finite entry is refused. |
 | `EigenGeneral(a)` | the eigenvalues and eigenvectors of a square matrix of any dtype; real and int inputs promote to `complex128`, values descend by magnitude, and a real matrix may carry conjugate pairs. |
 | `EigenGeneralised(a, b)` | solves A·v = λ·B·v for symmetric a and symmetric positive definite b, through the Cholesky factor of b; values ascend, eigenvectors are B-orthonormal columns. |
 | `SVD(a)` | the thin decomposition a = U·Σ·Vᵀ for m ≥ n: U is m×n with orthonormal columns, Σ a length-n vector descending, Vᵀ n×n orthogonal; a wide matrix is transposed first and the factors are swapped back. |
@@ -898,7 +898,7 @@ The conditions this package reports:
 - a modification that needs fill the factor does not hold:
   `SparseCholesky.Update` and `Downdate` refuse it and leave the
   factor as it was, and `CholeskyUpdate` refuses a rank-1 factor, a
-  mismatched vector or a complex input.
+  mismatched vector, a complex input or a non-finite entry.
 - a complex or wrong-length right-hand side to a sparse solve, and a
   preconditioner built for another dimension.
 - non-finite input to `NewSparseLU`, `NewSparseILU` and the sparse
@@ -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 |
 |---|---|

docs/benchmarks/release-vs-head.md

@@ -0,0 +1,76 @@
+# Benchmark report: v1.0.0 against head 045ef28
+
+One live comparison, regenerated by `just bench-report` and never
+accumulated: the newest release tag against the working tree, taken in
+one interleaved session on an idle machine. Re-run it before a release
+and commit the file the run writes.
+
+## Environment
+
+- CPU: AMD RYZEN AI MAX+ PRO 395 w/ Radeon 8060S, 32 logical cores
+- Memory: 118 GiB
+- OS: Linux 7.2.5-200.fc44.x86_64
+- Go: go version go1.27.1 linux/amd64
+- Build: the portable build, no pinned `GOAMD64` level, no `GOEXPERIMENT`
+
+## Revisions
+
+- Release: `v1.0.0`
+- Head: `045ef28`
+- The set: the benchmarks `bench_set` names that both revisions carry;
+  a name either side lacks is left out, never counted as a result.
+
+## Method
+
+Four interleaved rounds, the side order alternating between rounds, two
+counts per side per round, `-benchtime=0.5s` with `-benchmem`,
+medians over the eight samples a side collects. A factor at or above
+1.05 is faster, at or below 0.95 slower, anything between noise. Time
+of the release divides time of the head, so above one means the head is
+faster.
+
+## Results
+
+| Benchmark | v1.0.0 | head 045ef28 | factor | B/op release → head | allocs/op release → head |
+|---|---|---|---|---|---|
+| `KernelTransposeTiled256` | 87.58 µs | 72.89 µs | 1.20x | 524752 → 524752 | 4 → 4 |
+| `Prod2D` | 49.68 µs | 42.41 µs | 1.17x | 5416 → 5416 | 37 → 37 |
+| `MatMulSmall64` | 43.30 µs | 41.00 µs | 1.06x | 35040 → 35040 | 69 → 69 |
+| `Add1M` | 366.14 µs | 349.53 µs | 1.05x | 8005866 → 8005863 | 69 → 69 |
+| `CumSum2D` | 232.91 µs | 222.16 µs | 1.05x | 2098495 → 2098498 | 37 → 37 |
+| `FFT4096` | 26.88 µs | 25.52 µs | 1.05x | 65899 → 65899 | 3 → 3 |
+| `MatVec` | 29.53 µs | 28.21 µs | 1.05x | 5416 → 5416 | 37 → 37 |
+| `IntegrateHeat2D` | 936.53 µs | 909.54 µs | 1.03x | 92464 → 92464 | 159 → 159 |
+| `MinimiseLBFGSBounded` | 9.02 µs | 8.76 µs | 1.03x | 19016 → 19016 | 53 → 53 |
+| `MatMul1024` | 26.00 ms | 25.46 ms | 1.02x | 8392200 → 8391652 | 70 → 69 |
+| `MatMulOdd130` | 162.07 µs | 158.28 µs | 1.02x | 141205 → 141205 | 57 → 57 |
+| `MeanAxis2D` | 25.85 µs | 25.44 µs | 1.02x | 5464 → 5464 | 37 → 37 |
+| `Median` | 2.64 ms | 2.60 ms | 1.02x | 401414 → 401410 | 1 → 1 |
+| `Dot` | 84.49 µs | 83.68 µs | 1.01x | 1171 → 1174 | 36 → 36 |
+| `EinsumBatchedMatMul` | 314.97 µs | 312.57 µs | 1.01x | 528746 → 528745 | 78 → 78 |
+| `FFT2_512x512` | 959.31 µs | 948.74 µs | 1.01x | 4461272 → 4461249 | 168 → 168 |
+| `IntegrateRK4` | 887.97 µs | 880.62 µs | 1.01x | 3906104 → 3906102 | 24015 → 24015 |
+| `SavitzkyGolay` | 331.05 µs | 328.79 µs | 1.01x | 2099428 → 2099429 | 69 → 69 |
+| `ArgSort1M` | 10.03 ms | 10.07 ms | 1.00x | 42213251 → 42213264 | 266 → 266 |
+| `Conv1D` | 166.71 µs | 167.04 µs | 1.00x | 281028 → 280800 | 106 → 106 |
+| `CovarianceMatrix` | 101.98 µs | 102.44 µs | 1.00x | 136176 → 136176 | 23 → 23 |
+| `IntegrateBackwardEuler` | 89.86 µs | 89.96 µs | 1.00x | 352929 → 352929 | 2806 → 2806 |
+| `KernelDensity` | 241.38 µs | 241.85 µs | 1.00x | 4467 → 4462 | 69 → 69 |
+| `LogisticRegression` | 235.50 µs | 236.46 µs | 1.00x | 10656 → 10656 | 61 → 61 |
+| `Exp100k` | 182.19 µs | 183.11 µs | 0.99x | 805029 → 805029 | 69 → 69 |
+| `MatMulTall512x64x2048` | 2.09 ms | 2.10 ms | 0.99x | 8390997 → 8390963 | 69 → 69 |
+| `SumAxis3DDim1` | 2.60 µs | 2.62 µs | 0.99x | 2768 → 2768 | 4 → 4 |
+| `WelchPSD` | 197.78 µs | 198.78 µs | 0.99x | 28436 → 28603 | 75 → 75 |
+| `Norm2D` | 54.87 µs | 55.99 µs | 0.98x | 5432 → 5432 | 37 → 37 |
+| `SumAxis2D` | 26.27 µs | 26.80 µs | 0.98x | 5473 → 5473 | 37 → 37 |
+| `CWTMorlet` | 768.00 µs | 793.52 µs | 0.97x | 6427142 → 6427651 | 103 → 103 |
+| `LinearRegression` | 63.01 µs | 64.84 µs | 0.97x | 36152 → 36152 | 58 → 58 |
+| `MinimiseDifferentialEvolution` | 164.09 µs | 170.09 µs | 0.96x | 498250 → 498250 | 3845 → 3845 |
+| `SVD128` | 17.36 ms | 18.08 ms | 0.96x | 2182742 → 2181705 | 5070 → 5072 |
+| `BackwardTwoLayer` | 82.50 µs | 89.10 µs | 0.93x | 104759 → 104757 | 65 → 65 |
+| `Cholesky256` | 736.43 µs | 818.51 µs | 0.90x | 1583248 → 1583351 | 342 → 342 |
+| `Solve256` | 1.79 ms | 2.25 ms | 0.80x | 537635 → 537639 | 10 → 10 |
+
+## Summary
+
+4 of 37 benchmarks sit above the noise band (faster), 3 below it (slower), 30 inside it.

integrate/cubature.go

@@ -286,7 +286,7 @@ func IntegrateND(f func(x []float64) float64, lower, upper []float64, opts Cubat
 	total := rootVal
 	totalEst := rootEst
 	// The stopping rule scales the tolerance with the magnitude of the
-	// integral, the way Integrate combines its bounds: the error
+	// integral, the way IntegrateFunction combines its bounds: the error
 	// estimate of an integral of size 1e6 cannot fall below the
 	// rounding floor of the sum itself, so a purely absolute tolerance
 	// would burn the whole budget and report an exhausted budget

integrate/fem2d.go

@@ -114,13 +114,25 @@ func (m *TriangleMesh2D) BoundaryEdges() []int {
 		count[key(b, c)]++
 		count[key(c, a)]++
 	}
-	edges := make([]int, 0, 8)
+	sets := make([][2]int, 0, len(count))
 	for e, n := range count {
 		if n == 1 {
+			sets = append(sets, e)
+		}
+	}
+	// The pairs sort as pairs, never as one flat index list: a flat sort
+	// interleaves the endpoints of unrelated edges and hands back pairs
+	// the mesh does not carry.
+	slices.SortFunc(sets, func(x, y [2]int) int {
+		if x[0] != y[0] {
+			return x[0] - y[0]
+		}
+		return x[1] - y[1]
+	})
+	edges := make([]int, 0, 2*len(sets))
+	for _, e := range sets {
 		edges = append(edges, e[0], e[1])
 	}
-	}
-	slices.Sort(edges)
 	return edges
 }
 

integrate/fem2d_test.go

@@ -434,6 +434,34 @@ func TestTriangleMesh2DBoundaryEdges(t *testing.T) {
 	}
 }
 
+// TestTriangleMesh2DBoundaryEdgesAreMeshEdges pins the pair contract of
+// BoundaryEdges: every returned pair must be an edge the mesh actually
+// carries, and the pairs must come out sorted, which a sort of the flat
+// index list cannot deliver (it interleaves unrelated endpoints).
+func TestTriangleMesh2DBoundaryEdgesAreMeshEdges(t *testing.T) {
+	mesh, err := GridTriangleMesh2D(0, 0, 1, 1, 1, 1)
+	if err != nil {
+		t.Fatalf("GridTriangleMesh2D: %v", err)
+	}
+	edges := mesh.BoundaryEdges()
+	// The square's four sides: (0,1), (0,2), (1,3), (2,3) in sorted
+	// pair order.
+	want := []int{0, 1, 0, 2, 1, 3, 2, 3}
+	if len(edges) != len(want) {
+		t.Fatalf("boundary edge count %d, want %d", len(edges)/2, len(want)/2)
+	}
+	for p := 0; p < len(edges); p += 2 {
+		if edges[p] > edges[p+1] {
+			t.Fatalf("edge [%d,%d] is not sorted as a pair", edges[p], edges[p+1])
+		}
+	}
+	for p := range len(want) {
+		if edges[p] != want[p] {
+			t.Fatalf("boundary edges %v, want %v", edges, want)
+		}
+	}
+}
+
 // TestSolvePoissonFEM2DDuplicateDirichletNode pins the documented rule
 // for a node listed more than once: the last value is the prescribed
 // one and the node enters the assembled system exactly once, so the

integrate/quad.go

@@ -25,8 +25,8 @@ import (
 // rational substitution before the rule runs; an integrand over an
 // infinite interval has to decay to zero for this to converge.
 
-// QuadratureOptions tunes Integrate. RelTol ≤ 0 means 1e-10, AbsTol
+// QuadratureOptions tunes IntegrateFunction. RelTol ≤ 0 means 1e-10,
-// ≤ 0 means 1e-12, MaxIntervals ≤ 0 means 256.
+// AbsTol ≤ 0 means 1e-12, MaxIntervals ≤ 0 means 256.
 type QuadratureOptions struct {
 	RelTol       float64
 	AbsTol       float64
@@ -131,7 +131,7 @@ func IntegrateFunction(f func(x float64) (float64, error), a, b float64, opts Qu
 		opts.MaxIntervals = 256
 	}
 	if math.IsNaN(a) || math.IsNaN(b) {
-		return 0, 0, base.Errf("Integrate: bounds must not be NaN")
+		return 0, 0, base.Errf("IntegrateFunction: bounds must not be NaN")
 	}
 	sign := 1.0
 	if b < a {
@@ -229,7 +229,7 @@ func IntegrateFunction(f func(x float64) (float64, error), a, b float64, opts Qu
 
 	first, err := measure(lo, hi)
 	if err != nil {
-		return 0, 0, base.Errf("Integrate: %w", err)
+		return 0, 0, base.Errf("IntegrateFunction: %w", err)
 	}
 	leaves := []leaf{first}
 	budget := func() (total, errSum float64) {
@@ -242,7 +242,7 @@ func IntegrateFunction(f func(x float64) (float64, error), a, b float64, opts Qu
 	value, errSum := budget()
 	for errSum > math.Max(opts.AbsTol, opts.RelTol*math.Abs(value)) {
 		if len(leaves) >= opts.MaxIntervals {
-			return 0, 0, base.Errf("Integrate: error estimate %g exceeds the tolerance within %d subintervals",
+			return 0, 0, base.Errf("IntegrateFunction: error estimate %g exceeds the tolerance within %d subintervals",
 				errSum, opts.MaxIntervals)
 		}
 		// Bisect the leaf that contributes the most error.
@@ -255,11 +255,11 @@ func IntegrateFunction(f func(x float64) (float64, error), a, b float64, opts Qu
 		w := leaves[worst]
 		left, lerr := measure(w.l, (w.l+w.r)/2)
 		if lerr != nil {
-			return 0, 0, base.Errf("Integrate: %w", lerr)
+			return 0, 0, base.Errf("IntegrateFunction: %w", lerr)
 		}
 		right, rerr := measure((w.l+w.r)/2, w.r)
 		if rerr != nil {
-			return 0, 0, base.Errf("Integrate: %w", rerr)
+			return 0, 0, base.Errf("IntegrateFunction: %w", rerr)
 		}
 		leaves[worst] = left
 		leaves = append(leaves, right)

integrate/symplectic2.go

@@ -129,9 +129,10 @@ func verletPrelude(name string, q0, p0 *core.Array) (q, p []float64, n int, err
 
 // integerState reports whether dt is one of the integer-class state
 // dtypes the symplectic family refuses: bool and the narrow integer
-// widths follow Int into the standing "int states cannot integrate"
+// widths carry a discrete state with no place in a continuous
-// refusal, exactly as the round's follow-Int rule requires. The float
+// integrator, so they follow Int into the standing "int states cannot
-// dtypes, float16 included, keep the treatment they carry today.
+// integrate" refusal. The float dtypes, float16 included, keep the
+// treatment they carry today.
 func integerState(dt core.Dtype) bool {
 	switch dt {
 	case core.Bool, core.Int, core.Int8, core.Uint8, core.Int16, core.Uint16, core.Int32, core.Uint32:

linalg/cholupdate.go

@@ -71,14 +71,25 @@ func cholRankOne(l, x *core.Array, name string, sign int) (*core.Array, error) {
 	}
 	// The rotations below assume a lower triangular factor; a nonzero
 	// strict upper triangle would silently corrupt the sweep, so it is
-	// refused up front.
+	// refused up front, and so is a non-finite entry anywhere: the
+	// diagonal test cannot see a NaN, and the hyperbolic square of an
+	// Inf overflows mid-sweep. The sparse twin refuses the same input.
 	for i := range n {
-		for j := i + 1; j < n; j++ {
+		for j := range n {
-			if out.RawFloats()[i*n+j] != 0 {
+			w := out.RawFloats()[i*n+j]
+			if math.IsNaN(w) || math.IsInf(w, 0) {
+				return nil, base.Errf("%s: factor entry [%d,%d] is not finite", name, i, j)
+			}
+			if j > i && w != 0 {
 				return nil, base.Errf("%s: the factor must be lower triangular (nonzero entry at row %d, column %d)", name, i, j)
 			}
 		}
 	}
+	for i := range n {
+		if math.IsNaN(v[i]) || math.IsInf(v[i], 0) {
+			return nil, base.Errf("%s: entry %d is not finite", name, i)
+		}
+	}
 	// The modified matrix reads [L, x]·J·[L, x]ᵀ with J = diag(I, ±1):
 	// one sweep of (hyperbolic for −1, orthogonal for +1) rotations
 	// eliminates the vector column, and the surviving columns are the

linalg/cholupdate_test.go

@@ -5,8 +5,10 @@ package linalg
 
 import (
 	"math"
-	"sourcedock.dev/petrbalvin/tensor/internal/core"
+	"strings"
 	"testing"
+
+	"sourcedock.dev/petrbalvin/tensor/internal/core"
 )
 
 // spdSample builds a deterministic symmetric positive definite matrix:
@@ -156,3 +158,25 @@ func TestCholeskyRankOneErrors(t *testing.T) {
 		t.Fatal("expected an error for a complex vector")
 	}
 }
+
+// TestCholeskyRankOneRefusesNonFinite pins the refusal of a poisoned
+// factor or vector: the diagonal test cannot see a NaN, so the sweep
+// would answer an all-NaN factor with a nil error, where the sparse
+// twin refuses the same modification.
+func TestCholeskyRankOneRefusesNonFinite(t *testing.T) {
+	l, err := Cholesky(mustFromFloats(t, []float64{4, 0, 0, 1}, 2, 2))
+	if err != nil {
+		t.Fatalf("Cholesky: %v", err)
+	}
+	badVec := floatsToArray([]float64{math.NaN(), 0.1}, []int{2})
+	if _, err := CholeskyUpdate(l, badVec); err == nil || !strings.Contains(err.Error(), "not finite") {
+		t.Fatalf("CholeskyUpdate(NaN vector): %v", err)
+	}
+	if _, err := CholeskyDowndate(l, badVec); err == nil || !strings.Contains(err.Error(), "not finite") {
+		t.Fatalf("CholeskyDowndate(NaN vector): %v", err)
+	}
+	badFac := mustFromFloats(t, []float64{4, 0, math.Inf(1), 1}, 2, 2)
+	if _, err := CholeskyUpdate(badFac, floatsToArray([]float64{0.1, 0.1}, []int{2})); err == nil || !strings.Contains(err.Error(), "not finite") {
+		t.Fatalf("CholeskyUpdate(Inf factor): %v", err)
+	}
+}

linalg/decomp.go

@@ -48,7 +48,19 @@ func Cholesky(a *core.Array) (*core.Array, error) {
 		return nil, base.Errf("Cholesky: needs a square 2-D matrix, got shape %s", base.ShapeText(a.Shape()))
 	}
 	n := a.Shape()[0]
-	lMat, err := denseCholFactor(denseFloats(a, n, n), n)
+	aMat := denseFloats(a, n, n)
+	// A non-finite entry is neither positive nor definite, but the pivot
+	// test below cannot see it: NaN fails every comparison, so the sweep
+	// would return an all-NaN factor with a nil error. The sparse twin
+	// refuses the same input, and so does this one.
+	for i := range n {
+		for j := range n {
+			if math.IsNaN(aMat[i*n+j]) || math.IsInf(aMat[i*n+j], 0) {
+				return nil, base.Errf("Cholesky: entry [%d,%d] is not finite", i, j)
+			}
+		}
+	}
+	lMat, err := denseCholFactor(aMat, n)
 	if err != nil {
 		return nil, err
 	}

linalg/decomp2.go

@@ -706,6 +706,18 @@ func Eigen(a *core.Array) (values, vectors *core.Array, err error) {
 	}
 	n := a.Shape()[0]
 	mat := denseFloats(a, n, n)
+	// A non-finite entry slips through the symmetry test below: the
+	// difference of two NaNs never exceeds the tolerance, so a poisoned
+	// matrix reads as symmetric and the sweep hands back a NaN spectrum
+	// with a nil error. The sparse twin refuses the same input, and so
+	// does this one.
+	for i := range n {
+		for j := range n {
+			if math.IsNaN(mat[i*n+j]) || math.IsInf(mat[i*n+j], 0) {
+				return nil, nil, base.Errf("Eigen: entry [%d,%d] is not finite", i, j)
+			}
+		}
+	}
 	// A matrix outside the safe window would drive the reflector norms
 	// and the sweep's squared accumulation past the representable range;
 	// the tridiagonalisation and the QR sweep run on it scaled into the

linalg/decomp3.go

@@ -55,6 +55,18 @@ func EigenComplex(a *core.Array) (values, vectors *core.Array, err error) {
 	if !hermitianOK(h, n) {
 		return nil, nil, base.Errf("EigenComplex: matrix is not Hermitian within 1e-12 tolerance")
 	}
+	// A non-finite entry slips through the Hermitian test above: the
+	// difference of a NaN and its conjugate never exceeds the tolerance.
+	// The sparse twin refuses it in its own guard, and so does this one,
+	// before the Jacobi sweep burns its passes on a poisoned matrix.
+	for i := range n {
+		for j := range n {
+			z := h[i*n+j]
+			if math.IsNaN(real(z)) || math.IsNaN(imag(z)) || math.IsInf(real(z), 0) || math.IsInf(imag(z), 0) {
+				return nil, nil, base.Errf("EigenComplex: entry [%d,%d] is not finite", i, j)
+			}
+		}
+	}
 	// The Jacobi sweep's convergence test is a sum of squared
 	// magnitudes: outside the safe window it overflows to +Inf, which
 	// makes every sweep look already converged, or underflows to 0,
@@ -110,6 +122,17 @@ func hermitianJacobi(h []complex128, v []complex128, n int) error {
 	// the floor already diagonal and return its raw diagonal, identity
 	// eigenvectors included.
 	threshold := 1e-13 * norm
+	// The per-entry skip must sit where a matrix whose every
+	// off-diagonal entry is skippable also passes the aggregate
+	// convergence test: with P pairs the off-norm of all-skipped entries
+	// reaches skip·√P, so a skip at the plain threshold leaves a matrix
+	// whose entries are uniformly just under it frozen above the test,
+	// and the sweep exhausts its passes without turning a single
+	// rotation.
+	skip := threshold
+	if pairs := float64(n * (n - 1) / 2); pairs > 1 {
+		skip = threshold / math.Sqrt(pairs)
+	}
 	offMass := func() float64 {
 		off := 0.0
 		for p := range n {
@@ -127,7 +150,7 @@ func hermitianJacobi(h []complex128, v []complex128, n int) error {
 		for p := range n {
 			for q := p + 1; q < n; q++ {
 				z := h[p*n+q]
-				if base.AbsComplex(z) <= threshold {
+				if base.AbsComplex(z) <= skip {
 					continue
 				}
 				// Step 1: a diagonal phase turn makes the pair entry

linalg/decomp3_test.go

@@ -6,9 +6,11 @@ package linalg
 import (
 	"math"
 	"math/cmplx"
+	"strings"
+	"testing"
+
 	"sourcedock.dev/petrbalvin/tensor/internal/base"
 	"sourcedock.dev/petrbalvin/tensor/internal/core"
-	"testing"
 )
 
 func mustComplexes(t *testing.T, vals []complex128, shape ...int) *core.Array {
@@ -291,3 +293,53 @@ func subComplex(a, b []complex128) []complex128 {
 	}
 	return out
 }
+
+// TestEigenComplexNearDiagonalConverges pins the Jacobi sweep against a
+// matrix whose every off-diagonal entry sits just under the per-entry
+// skip level while the aggregate off-norm stays above the convergence
+// threshold: the skip must not freeze the sweep above its own
+// convergence test, which used to exhaust the passes and report the
+// nearly diagonal matrix as unconverged.
+func TestEigenComplexNearDiagonalConverges(t *testing.T) {
+	const n = 30
+	cv := make([]complex128, n*n)
+	for i := range n {
+		cv[i*n+i] = 1
+	}
+	for i := range n {
+		for j := i + 1; j < n; j++ {
+			cv[i*n+j] = complex(0.9e-13, 0)
+			cv[j*n+i] = complex(0.9e-13, 0)
+		}
+	}
+	a, err := core.FromComplexes(cv, n, n)
+	if err != nil {
+		t.Fatalf("FromComplexes: %v", err)
+	}
+	vals, _, err := EigenComplex(a)
+	if err != nil {
+		t.Fatalf("EigenComplex: %v", err)
+	}
+	for i := range n {
+		if math.Abs(vals.FloatAt(i)-1) > 1e-11 {
+			t.Fatalf("eigenvalue %d is %.12g, want 1", i, vals.FloatAt(i))
+		}
+	}
+}
+
+// TestEigenComplexRefusesNonFinite pins that a poisoned matrix never
+// reads as Hermitian: the mirror comparison cannot see a NaN
+// difference, so the entry is refused outright, the way the sparse
+// sibling's Hermitian check refuses it.
+func TestEigenComplexRefusesNonFinite(t *testing.T) {
+	cv := []complex128{complex(math.NaN(), 0), 0, 0, 1}
+	a := mustComplexes(t, cv, 2, 2)
+	if _, _, err := EigenComplex(a); err == nil || !strings.Contains(err.Error(), "not finite") {
+		t.Fatalf("EigenComplex(NaN): %v", err)
+	}
+	cv2 := []complex128{complex(0, math.Inf(1)), 0, 0, 1}
+	b := mustComplexes(t, cv2, 2, 2)
+	if _, _, err := EigenComplex(b); err == nil || !strings.Contains(err.Error(), "not finite") {
+		t.Fatalf("EigenComplex(Inf): %v", err)
+	}
+}

linalg/decomp_test.go

@@ -5,8 +5,10 @@ package linalg
 
 import (
 	"math"
-	"sourcedock.dev/petrbalvin/tensor/internal/core"
+	"strings"
 	"testing"
+
+	"sourcedock.dev/petrbalvin/tensor/internal/core"
 )
 
 // TestSVDReconstruction checks that A = U · Σ · Vᵀ reconstructs A
@@ -477,3 +479,18 @@ func TestHouseholderVectorIntoLargeScale(t *testing.T) {
 		t.Fatalf("beta = %v for a zero vector, want 0", empty.beta)
 	}
 }
+
+// TestEigenRefusesNonFinite pins that a poisoned matrix never reads as
+// symmetric: the guard's comparison cannot see a NaN difference, so the
+// entry is refused outright, the way the sparse sibling's symmetry
+// check refuses it.
+func TestEigenRefusesNonFinite(t *testing.T) {
+	nan := mustFromFloats(t, []float64{math.NaN(), 0, 0, 1}, 2, 2)
+	if _, _, err := Eigen(nan); err == nil || !strings.Contains(err.Error(), "not finite") {
+		t.Fatalf("Eigen(NaN): %v", err)
+	}
+	inf := mustFromFloats(t, []float64{math.Inf(1), 0, 0, 1}, 2, 2)
+	if _, _, err := Eigen(inf); err == nil || !strings.Contains(err.Error(), "not finite") {
+		t.Fatalf("Eigen(Inf): %v", err)
+	}
+}

linalg/decompositions_test.go

@@ -691,3 +691,18 @@ func TestCholeskyBlockedReconstruction(t *testing.T) {
 		}
 	}
 }
+
+// TestCholeskyRefusesNonFinite pins the refusal of a poisoned matrix:
+// the pivot test cannot see a NaN (it fails every comparison), so
+// without the gate the sweep would answer an all-NaN factor with a nil
+// error, where the sparse sibling refuses the same input.
+func TestCholeskyRefusesNonFinite(t *testing.T) {
+	nan := mustFromFloats(t, []float64{math.NaN(), 0, 0, 1}, 2, 2)
+	if _, err := Cholesky(nan); err == nil || !strings.Contains(err.Error(), "not finite") {
+		t.Fatalf("Cholesky(NaN): %v", err)
+	}
+	inf := mustFromFloats(t, []float64{0, 0, 0, math.Inf(1)}, 2, 2)
+	if _, err := Cholesky(inf); err == nil || !strings.Contains(err.Error(), "not finite") {
+		t.Fatalf("Cholesky(Inf): %v", err)
+	}
+}

optim/devolution.go

@@ -58,8 +58,12 @@ func MinimiseDifferentialEvolution(f func(*core.Array) (float64, error),
 		return nil, 0, base.Errf("%s: complex bounds are not supported", name)
 	}
 	for i := range n {
-		if !(upper.FloatAt(i) > lower.FloatAt(i)) {
+		lo, up := lower.FloatAt(i), upper.FloatAt(i)
-			return nil, 0, base.Errf("%s: bound %d runs from %g to %g", name, i, lower.FloatAt(i), upper.FloatAt(i))
+		// An infinite side has no uniform draw: the population would be
+		// born NaN and the best point would come back NaN with no error,
+		// so a non-finite bound is degenerate exactly as a crossed one is.
+		if math.IsNaN(lo) || math.IsNaN(up) || math.IsInf(lo, 0) || math.IsInf(up, 0) || !(up > lo) {
+			return nil, 0, base.Errf("%s: bound %d runs from %g to %g", name, i, lo, up)
 		}
 	}
 	pop := opts.Population

optim/devolution_test.go

@@ -112,3 +112,25 @@ func TestDEErrors(t *testing.T) {
 		t.Error("NaN objective accepted")
 	}
 }
+
+// TestDERejectsInfiniteBounds pins that an infinite side is refused: the
+// sampler draws uniformly inside the box, so an open side would populate
+// it with NaN and hand back a NaN point with a nil error.
+func TestDERejectsInfiniteBounds(t *testing.T) {
+	zero, _ := core.FromFloats([]float64{0}, 1)
+	one, _ := core.FromFloats([]float64{1}, 1)
+	constF := func(a *core.Array) (float64, error) { return 1, nil }
+	for _, tc := range []struct {
+		name   string
+		lo, hi *core.Array
+	}{
+		{"both sides open", mustBound(math.Inf(-1)), mustBound(math.Inf(1))},
+		{"upper open", zero, mustBound(math.Inf(1))},
+		{"lower open", mustBound(math.Inf(-1)), one},
+	} {
+		x, _, err := MinimiseDifferentialEvolution(constF, tc.lo, tc.hi, DifferentialEvolutionOptions{Generations: 5})
+		if err == nil {
+			t.Errorf("%s: infinite bound accepted, returned x=%g", tc.name, x.FloatAt(0))
+		}
+	}
+}

plot/plot.go

@@ -119,11 +119,37 @@ func (c Chart) WriteSVG(path string) error {
 	if !(yr[0] < yr[1]) || math.IsInf(yr[0], 0) || math.IsInf(yr[1], 0) {
 		yr = bounds(all, false)
 	}
+	// project maps a value of the range [r0, r1] onto the plot span
+	// starting at plotLo. The normal form is the arithmetic the chart
+	// has always run; a span, or an offset from the range's start, that
+	// overflows float64 falls back to the halved form, whose every term
+	// a finite input keeps finite, so a range the data spans but the
+	// arithmetic cannot still draws instead of carrying a coordinate no
+	// renderer displays.
+	project := func(v, r0, r1, plotLo, plotSpan float64) float64 {
+		u := (v - r0) / (r1 - r0)
+		if math.IsNaN(u) || math.IsInf(u, 0) {
+			u = (v/2 - r0/2) / (r1/2 - r0/2)
+		}
+		return plotLo + u*plotSpan
+	}
 	px := func(x float64) float64 {
-		return left + (x-xr[0])/(xr[1]-xr[0])*(float64(w)-left-right)
+		return project(x, xr[0], xr[1], left, float64(w)-left-right)
 	}
 	py := func(y float64) float64 {
-		return float64(h) - bottom - (y-yr[0])/(yr[1]-yr[0])*(float64(h)-top-bottom)
+		return float64(h) - bottom - project(y, yr[0], yr[1], 0, float64(h)-top-bottom)
+	}
+	// tickValue places the k-th of the five ticks. A span that overflows
+	// float64 makes the affine form NaN or Inf, so the tick falls back to
+	// the convex combination, which stays between the range's own finite
+	// ends.
+	tickValue := func(r0, r1 float64, k int) float64 {
+		t := r0 + (r1-r0)*float64(k)/4
+		if math.IsNaN(t) || math.IsInf(t, 0) {
+			f := float64(k) / 4
+			t = r0*(1-f) + r1*f
+		}
+		return t
 	}
 	var b strings.Builder
 	b.WriteString(xmlHeader)
@@ -133,7 +159,7 @@ func (c Chart) WriteSVG(path string) error {
 		left, esc(c.Title))
 	// Axes with five ticks each.
 	for k := range 5 {
-		t := xr[0] + (xr[1]-xr[0])*float64(k)/4
+		t := tickValue(xr[0], xr[1], k)
 		x := px(t)
 		fmt.Fprintf(&b, "<line x1=\"%g\" y1=\"%g\" x2=\"%g\" y2=\"%g\" stroke=\"#ccc\" stroke-width=\"1\"/>\n",
 			x, top, x, float64(h)-bottom)
@@ -141,7 +167,7 @@ func (c Chart) WriteSVG(path string) error {
 			x, float64(h)-bottom+16, tick(t))
 	}
 	for k := range 5 {
-		t := yr[0] + (yr[1]-yr[0])*float64(k)/4
+		t := tickValue(yr[0], yr[1], k)
 		y := py(t)
 		fmt.Fprintf(&b, "<line x1=\"%g\" y1=\"%g\" x2=\"%g\" y2=\"%g\" stroke=\"#ccc\" stroke-width=\"1\"/>\n",
 			left, y, float64(w)-right, y)
@@ -195,8 +221,14 @@ func bounds(pts []Point, xAxis bool) [2]float64 {
 	if hi == lo {
 		hi = lo + 1
 	}
-	pad := 0.05 * (hi - lo)
+	padLo, padHi := lo-0.05*(hi-lo), hi+0.05*(hi-lo)
-	return [2]float64{lo - pad, hi + pad}
+	if math.IsInf(padLo, 0) || math.IsInf(padHi, 0) {
+		// The padding, or the span it scales, overflows the range the
+		// finite data itself fits; the unpadded bounds keep every
+		// projection finite.
+		return [2]float64{lo, hi}
+	}
+	return [2]float64{padLo, padHi}
 }
 
 func tick(v float64) string {

plot/plot_test.go

@@ -103,6 +103,39 @@ func TestWriteSVGErrors(t *testing.T) {
 	}
 }
 
+// TestWriteSVGExtremeFiniteValuesStayFinite pins the renderer against
+// the extreme-but-finite corner: data, or an explicit axis range, whose
+// span overflows float64 must not leak NaN or Inf coordinates into the
+// drawing, because no renderer displays them and the non-finite
+// refusal at the door guarantees every input point is finite.
+func TestWriteSVGExtremeFiniteValuesStayFinite(t *testing.T) {
+	dir := t.TempDir()
+	padded := sampleChart()
+	padded.XRange = [2]float64{-1e308, 1e308}
+	padded.YRange = [2]float64{-1e308, 1e308}
+	charts := []struct {
+		name  string
+		chart Chart
+	}{
+		{"wide-data.svg", Chart{Series: []Series{{Name: "wide",
+			Points: []Point{{X: 0, Y: -1e308}, {X: 1, Y: 1e308}}}}}},
+		{"wide-range.svg", padded},
+	}
+	for _, tc := range charts {
+		path := filepath.Join(dir, tc.name)
+		if err := tc.chart.WriteSVG(path); err != nil {
+			t.Fatalf("%s: %v", tc.name, err)
+		}
+		body, err := os.ReadFile(path)
+		if err != nil {
+			t.Fatalf("%s: %v", tc.name, err)
+		}
+		if s := string(body); strings.Contains(s, "NaN") || strings.Contains(s, "Inf") {
+			t.Fatalf("%s: the rendering leaked a non-finite coordinate: %s", tc.name, s)
+		}
+	}
+}
+
 func mustFromFloats(t *testing.T, values []float64, shape ...int) *core.Array {
 	t.Helper()
 	a, err := core.FromFloats(values, shape...)

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 {

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
 }
 

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.

stats/cdf.go

@@ -269,15 +269,18 @@ func StudentTCDF(t float64, df int) (float64, error) {
 	if df < 1 {
 		return 0, base.Errf("StudentTCDF: df must be ≥ 1, got %d", df)
 	}
-	z := float64(df) / (float64(df) + t*t)
+	// One house tail: the upper-tail helper carries the asymptotic forms
-	upper, err := BetaIncomplete(z, float64(df)/2, 0.5)
+	// the heavy df ≤ 2 laws keep past t²'s overflow, where the closed
+	// form's z = df/(df+t²) collapses to 0 and the lower tail answered a
+	// silent 0 for a tail the format still holds.
+	upper, err := studentTUpperTail(math.Abs(t), df)
 	if err != nil {
 		return 0, base.Errf("StudentTCDF: %w", err)
 	}
 	if t >= 0 {
-		return 1 - upper/2, nil
+		return 1 - upper, nil
 	}
-	return upper / 2, nil
+	return upper, nil
 }
 
 // PoissonCDF returns P(N ≤ k) for N ~ Poisson(lambda), through the
@@ -701,7 +704,10 @@ func studentTUpperTail(t float64, df int) (float64, error) {
 			// huge t keeps accurate.
 			return 1 / (math.Pi * t), nil
 		case df == 2:
-			return 1 / (2 * t * t), nil
+			// The tail 1/(2t²) divided one t at a time: the literal
+			// denominator overflows past √MaxFloat64, and the quotient
+			// would flush the still-representable tail to zero.
+			return 0.5 / t / t, nil
 		}
 	}
 	z := float64(df) / (float64(df) + t*t)

stats/cdf_test.go

@@ -276,3 +276,56 @@ func TestGammaLowerLargeShape(t *testing.T) {
 		}
 	}
 }
+
+// TestStudentTCDFSquareOverflowTail pins the Student laws past the square
+// overflow: a t beyond √MaxFloat64 drives z = df/(df+t²) through Inf/Inf
+// to the NaN and 0 the closed form answers, while the heavy df = 1 and
+// df = 2 tails are still representable there. The lower tail and the
+// two-sided tail must answer their asymptotic forms, exactly as the
+// one-sided upper tail already does.
+func TestStudentTCDFSquareOverflowTail(t *testing.T) {
+	const huge = 1e155
+	// df = 1 (Cauchy): the one-sided tail is 1/(π·t), two-sided 2/(π·t).
+	one, err := StudentTCDF(-huge, 1)
+	if err != nil {
+		t.Fatalf("StudentTCDF(-1e155, 1): %v", err)
+	}
+	if want := 1 / (math.Pi * huge); math.Abs(one-want) > 1e-12*want {
+		t.Fatalf("StudentTCDF(-1e155, 1) = %.17g, want the Cauchy tail %.17g", one, want)
+	}
+	two, err := twoSidedT(huge, 1)
+	if err != nil {
+		t.Fatalf("twoSidedT(1e155, 1): %v", err)
+	}
+	if want := 2 / (math.Pi * huge); math.Abs(two-want) > 1e-12*want {
+		t.Fatalf("twoSidedT(1e155, 1) = %.17g, want the Cauchy tail %.17g", two, want)
+	}
+	// df = 2: the two-sided tail is 1 − t/√(t²+2), the one-sided half of
+	// it, both ≈ 1/t² here and still inside the subnormal range. The
+	// reference is assembled as u/(√(1+u)+1) with u = 2/t², the form
+	// that never squares t.
+	const u = 2e-310 // 2/t² at t = 1e155
+	two2, err := twoSidedT(huge, 2)
+	if err != nil {
+		t.Fatalf("twoSidedT(1e155, 2): %v", err)
+	}
+	if want := u / (math.Sqrt(1+u) + 1); math.Abs(two2-want) > 1e-6*want {
+		t.Fatalf("twoSidedT(1e155, 2) = %.17g, want the df 2 tail %.17g", two2, want)
+	}
+	one2, err := StudentTCDF(-huge, 2)
+	if err != nil {
+		t.Fatalf("StudentTCDF(-1e155, 2): %v", err)
+	}
+	if want := 0.5 * u / (math.Sqrt(1+u) + 1); math.Abs(one2-want) > 1e-6*want {
+		t.Fatalf("StudentTCDF(-1e155, 2) = %.17g, want the df 2 tail %.17g", one2, want)
+	}
+	// The upper half of the axis keeps answering 1, and a df whose tail
+	// genuinely underflows keeps answering 0: both are the honest
+	// roundings there.
+	if v, err := StudentTCDF(huge, 1); err != nil || v != 1 {
+		t.Fatalf("StudentTCDF(1e155, 1) = %v (%v), want 1", v, err)
+	}
+	if v, err := StudentTCDF(-huge, 4); err != nil || v != 0 {
+		t.Fatalf("StudentTCDF(-1e155, 4) = %v (%v), want the underflowed 0", v, err)
+	}
+}

stats/hmm.go

@@ -354,6 +354,12 @@ func (m *HiddenMarkovModel) Viterbi(observations []int) (states []int, logProbab
 			arg = k
 		}
 	}
+	if math.IsInf(best, -1) {
+		// Every path has probability zero, the structural zero an
+		// emission or transition carries: the same refusal Forward makes,
+		// not a meaningless path beside a −Inf log probability.
+		return nil, 0, base.Errf("%s: the observation sequence has probability zero under the model", name)
+	}
 	path := make([]int, len(observations))
 	path[len(observations)-1] = arg
 	for t := len(observations) - 1; t > 0; t-- {
@@ -524,6 +530,15 @@ func hmmReestimate(model *HiddenMarkovModel, observations []int, gamma, xi [][]f
 		for j := range states {
 			den += transition[k*states+j]
 		}
+		if den == 0 {
+			// No transition evidence reached this row, the single-
+			// observation sequence being the case in point: dividing the
+			// zero fills the row with NaN the constructor refuses, and
+			// the sweep after it crashed on the nil model that refusal
+			// left. The row keeps its previous estimate instead.
+			copy(transition[k*states:(k+1)*states], model.Transition[k*states:(k+1)*states])
+			continue
+		}
 		for j := range states {
 			transition[k*states+j] = math.Max(transition[k*states+j]/den, hmmFloor)
 		}

stats/hmm_test.go

@@ -343,3 +343,61 @@ func TestHiddenMarkovZeroProbabilitySequence(t *testing.T) {
 		t.Fatalf("the certain sequence answered (%g, %v), want (0, [1])", ll, filtered[0])
 	}
 }
+
+func TestHiddenMarkovFitSingleObservation(t *testing.T) {
+	// One observation carries emission and initial evidence but no
+	// transition evidence: the re-estimation divided the zero count into
+	// NaN rows, the constructor refused them, and the sweep then called
+	// a method on the nil model and crashed the process. The fit must
+	// answer with a valid model whose transitions keep their starting
+	// estimate.
+	g := core.NewGenerator(11)
+	res, err := FitHiddenMarkovModel(g, []int{0}, 2, 2)
+	if err != nil {
+		t.Fatalf("FitHiddenMarkovModel on one observation: %v", err)
+	}
+	if res.Model == nil {
+		t.Fatal("FitHiddenMarkovModel on one observation returned no model")
+	}
+	if math.IsNaN(res.LogLikelihood) || math.IsInf(res.LogLikelihood, 0) {
+		t.Fatalf("log likelihood = %g, want a finite value", res.LogLikelihood)
+	}
+	for k := range 2 {
+		sum := 0.0
+		for _, v := range res.Model.Transition[k*2 : k*2+2] {
+			if math.IsNaN(v) || v < 0 {
+				t.Fatalf("transition row %d holds %g, want probabilities", k, v)
+			}
+			sum += v
+		}
+		if math.Abs(sum-1) > 1e-9 {
+			t.Fatalf("transition row %d sums to %g, want 1", k, sum)
+		}
+	}
+}
+
+func TestViterbiZeroProbabilitySequence(t *testing.T) {
+	// Symbol 2 is emitted by no state, so every path through the sequence
+	// has probability zero. Forward and Smooth refuse such a sequence, and
+	// Viterbi documents the same refusals: a meaningless path beside a
+	// −Inf log probability is not an answer.
+	model, err := NewHiddenMarkovModel(
+		[]float64{0.5, 0.5},
+		[]float64{0.5, 0.5, 0.5, 0.5},
+		[]float64{0.5, 0.5, 0, 0.5, 0.5, 0},
+	)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if _, _, err := model.Viterbi([]int{0, 2, 1}); err == nil {
+		t.Fatal("a zero-probability sequence was accepted by Viterbi")
+	}
+	// The live symbols of the same model still decode.
+	path, lp, err := model.Viterbi([]int{0, 1})
+	if err != nil {
+		t.Fatal(err)
+	}
+	if path[0] != path[1] || lp >= 0 {
+		t.Fatalf("Viterbi([0, 1]) = (%v, %g), want a coherent path with a negative log probability", path, lp)
+	}
+}

stats/noncentral.go

@@ -304,7 +304,14 @@ func NoncentralTCDF(t float64, df int, delta float64) (float64, error) {
 		magnitude = -t
 		shift = -delta
 	}
-	x := magnitude * magnitude / (magnitude*magnitude + float64(df))
+	magnitude2 := magnitude * magnitude
+	x := magnitude2 / (magnitude2 + float64(df))
+	if math.IsInf(magnitude2, 1) {
+		// A finite t whose square overflows drove the quotient through
+		// Inf/Inf into a NaN the incomplete beta refused under its own
+		// name. The beta argument's limit there is exactly 1.
+		x = 1
+	}
 	if x == 0 {
 		// t = 0: the value collapses to Φ(−δ) exactly.
 		return NormalCDF(-delta), nil

stats/noncentral_test.go

@@ -219,6 +219,36 @@ func TestNoncentralUnderflowSurvival(t *testing.T) {
 	}
 }
 
+// TestNoncentralTCDFHugeFiniteT pins the far corner of the signed axis: a
+// finite t whose square overflows drives the beta argument to Inf/Inf, a
+// NaN the incomplete beta refused under its own name. The CDF there is 1
+// below rounding for t on the δ side and 0 above it, the same limits the
+// central law answers.
+func TestNoncentralTCDFHugeFiniteT(t *testing.T) {
+	for _, c := range []struct {
+		tv    float64
+		df    int
+		delta float64
+		want  float64
+	}{
+		{1e200, 3, 2, 1},
+		{1e155, 1, 0.5, 1},
+		{-1e200, 5, 1, 0},
+		{-1e155, 2, -3, 0},
+	} {
+		got, err := NoncentralTCDF(c.tv, c.df, c.delta)
+		if err != nil {
+			t.Fatalf("NoncentralTCDF(%g, %d, %g): %v", c.tv, c.df, c.delta, err)
+		}
+		if math.IsNaN(got) || got < 0 || got > 1 {
+			t.Fatalf("NoncentralTCDF(%g, %d, %g) = %g, want a probability", c.tv, c.df, c.delta, got)
+		}
+		if math.Abs(got-c.want) > 1e-15 {
+			t.Fatalf("NoncentralTCDF(%g, %d, %g) = %.17g, want %g", c.tv, c.df, c.delta, got, c.want)
+		}
+	}
+}
+
 // TestNoncentralTIdentityReductions pins the exact corners: δ = 0 is
 // the central Student t, t = 0 is Φ(−δ), and the two reflection
 // identities of the law hold to rounding.

stats/regression.go

@@ -177,6 +177,7 @@ func LinearRegression(x, y *core.Array) (*LinearRegressionResult, error) {
 	rss := 0.0
 	tss := 0.0
 	uncentred := 0.0
+	maxRes, maxDev, maxY := 0.0, 0.0, 0.0
 	mean := 0.0
 	if fy != nil {
 		for _, v := range fy[:n] {
@@ -212,26 +213,85 @@ func LinearRegression(x, y *core.Array) (*LinearRegressionResult, error) {
 		rss += res * res
 		tss += (yv - mean) * (yv - mean)
 		uncentred += yv * yv
+		if a := math.Abs(res); a > maxRes {
+			maxRes = a
+		}
+		if a := math.Abs(yv - mean); a > maxDev {
+			maxDev = a
+		}
+		if a := math.Abs(yv); a > maxY {
+			maxY = a
+		}
 	}
 	if !hasConstant {
 		// The null model is y = 0, so the uncentred total is what the
 		// model has to beat, and it carries n degrees of freedom.
 		tss = uncentred
 	}
+	// The factored sums of squares: a response on a scale whose squared
+	// deviations fall below the subnormal floor reads as a zero sum while
+	// its deviations are live, and the statistics below would report the
+	// evidence backwards (an exact fit the t statistics cannot support,
+	// an F of zero beside them). Each pair keeps the largest deviation as
+	// the scale and the scaled sum as the unit, so scale²·unit is the
+	// true sum wherever the plain product underflows; the unit stays 1
+	// whenever the plain sum already holds.
+	rssScale, rssUnit := 1.0, rss
+	if rss == 0 && maxRes > 0 {
+		rssScale, rssUnit = maxRes, 0.0
+		for _, res := range out.Residuals {
+			d := res / maxRes
+			rssUnit += d * d
+		}
+	}
+	tssScale, tssUnit := 1.0, tss
+	if tss == 0 {
+		devScale := maxDev
+		if !hasConstant {
+			devScale = maxY
+		}
+		if devScale > 0 {
+			tssScale = devScale
+			tssUnit = 0.0
+			for r := range n {
+				var yv, dev float64
+				if fy != nil {
+					yv = fy[r]
+				} else {
+					yv = y.FloatAt(r)
+				}
+				if !hasConstant {
+					dev = yv
+				} else {
+					dev = yv - mean
+				}
+				d := dev / devScale
+				tssUnit += d * d
+			}
+		}
+	}
 	dof := n - p
 	out.ResidualVariance = rss / float64(dof)
-	if tss == 0 {
+	tssDOF := n - 1
+	if !hasConstant {
+		tssDOF = n
+	}
+	if tss == 0 && tssScale == 1 {
 		// A constant response reproduced exactly: R² is 1 by the
 		// perfect-fit convention, not the 1 − 0/0 NaN every consumer
 		// would propagate. The same guard the F statistic below has.
 		out.RSquared = 1
 		out.AdjustedRSquared = 1
+	} else if tss == 0 {
+		// The total underflowed while the response varies: the ratio of
+		// the factored forms, the scale factors divided out one at a
+		// time. Both R² measures round back to 1 here, but the F
+		// statistic below reads the same factored pieces and does not.
+		ratio := rssUnit / tssUnit * (rssScale / tssScale) * (rssScale / tssScale)
+		out.RSquared = 1 - ratio
+		out.AdjustedRSquared = 1 - ratio*float64(tssDOF)/float64(dof)
 	} else {
 		out.RSquared = 1 - rss/tss
-		tssDOF := n - 1
-		if !hasConstant {
-			tssDOF = n
-		}
 		out.AdjustedRSquared = 1 - (rss/float64(dof))/(tss/float64(tssDOF))
 	}
 	out.DModel = p - 1
@@ -271,6 +331,24 @@ func LinearRegression(x, y *core.Array) (*LinearRegressionResult, error) {
 			}
 			out.PValues[j] = pv
 		case v == 0:
+			// The residual sum of squares may have underflowed while the
+			// residuals live: the factored standard error is representable
+			// where the squared one is not, and the t test then reports
+			// the evidence it actually holds instead of an unearned
+			// infinity.
+			if rssScale != 1 && inv[j][j] > 0 {
+				se := rssScale * math.Sqrt(rssUnit*inv[j][j]/float64(dof))
+				if se > 0 {
+					out.StandardErrors[j] = se
+					out.TStatistics[j] = beta[j] / se
+					pv, err := twoSidedT(out.TStatistics[j], dof)
+					if err != nil {
+						return nil, base.Errf("%s: %w", name, err)
+					}
+					out.PValues[j] = pv
+					continue
+				}
+			}
 			// An exact fit: the coefficient is infinitely many standard
 			// errors from zero, and the evidence is total. Reporting
 			// t = 0 next to p = 0 would contradict itself. A zero
@@ -300,6 +378,19 @@ func LinearRegression(x, y *core.Array) (*LinearRegressionResult, error) {
 			explained = 0 // rounding only, and a negative F is meaningless
 		}
 		out.FStatistic = explained / float64(out.DModel) / out.ResidualVariance
+		if (math.IsInf(out.FStatistic, 0) || math.IsNaN(out.FStatistic)) && (rssScale != 1 || tssScale != 1) {
+			// An underflowed sum of squares drove the quotient to Inf or
+			// 0/0 while the factored pieces live: F from the factored
+			// forms, every scale factor applied one division at a time so
+			// no intermediate leaves the representable range before the
+			// answer does. rssUnit 0 is the exact fit, whose F is
+			// genuinely infinite.
+			out.FStatistic = (tssUnit*tssScale/rssScale/rssScale*tssScale - rssUnit) *
+				float64(dof) / (float64(out.DModel) * rssUnit)
+			if out.FStatistic < 0 {
+				out.FStatistic = 0
+			}
+		}
 		if math.IsNaN(out.FStatistic) {
 			// 0/0: a response with no variation at all, reproduced
 			// exactly by the fit. There is no evidence of a model, so
@@ -335,14 +426,16 @@ func LinearRegression(x, y *core.Array) (*LinearRegressionResult, error) {
 // freedom, by the closed-form tail I_z(df/2, 1/2) with z = df/(df+t²).
 // The identity is used rather than 2·(1 − T_cdf(t)): near t = 0 the
 // subtraction cancels catastrophically, while the incomplete beta
-// stays accurate into the far tail where p values matter most.
+// stays accurate into the far tail where p values matter most. The
+// upper-tail helper carries it, so the asymptotic forms that hold the
+// df ≤ 2 tails past t²'s overflow serve here too: the bare closed form
+// answers a silent 0 there while the true tail is still representable.
 func twoSidedT(t float64, df int) (float64, error) {
-	z := float64(df) / (float64(df) + t*t)
+	upper, err := studentTUpperTail(math.Abs(t), df)
-	p, err := BetaIncomplete(z, float64(df)/2, 0.5)
 	if err != nil {
 		return 0, err
 	}
-	return p, nil
+	return 2 * upper, nil
 }
 
 // hasConstantColumn reports whether an (n, p) design holds a column of
@@ -482,6 +575,7 @@ func WeightedLinearRegression(x, y, w *core.Array) (*LinearRegressionResult, err
 	// regression and reports the uncentred conventions whenever the
 	// weights vary, so the four model-level fields are overwritten here.
 	sumW, sumWY, rssW := 0.0, 0.0, 0.0
+	maxResW := 0.0
 	for r := range n {
 		var wr float64
 		if fw != nil {
@@ -498,8 +592,12 @@ func WeightedLinearRegression(x, y, w *core.Array) (*LinearRegressionResult, err
 		sumW += wr
 		sumWY += wr * yv
 		rssW += wr * out.Residuals[r] * out.Residuals[r]
+		if a := math.Abs(out.Residuals[r]); a > maxResW {
+			maxResW = a
+		}
 	}
 	tssW := 0.0
+	maxDevW := 0.0
 	if hasConstant {
 		meanW := sumWY / sumW
 		for r := range n {
@@ -516,6 +614,9 @@ func WeightedLinearRegression(x, y, w *core.Array) (*LinearRegressionResult, err
 			}
 			d := yv - meanW
 			tssW += wr * d * d
+			if a := math.Abs(d); a > maxDevW {
+				maxDevW = a
+			}
 		}
 	} else {
 		// Without an intercept the null model is zero, so Σw·y² is the
@@ -533,16 +634,85 @@ func WeightedLinearRegression(x, y, w *core.Array) (*LinearRegressionResult, err
 				yv = y.FloatAt(r)
 			}
 			tssW += wr * yv * yv
+			if a := math.Abs(yv); a > maxDevW {
+				maxDevW = a
+			}
+		}
+	}
+	// The weighted sums of squares carry the same factored form the
+	// unweighted fit keeps: a response scale whose weighted squared
+	// deviations fall below the subnormal floor reads as a zero sum
+	// while the deviations live, and the F below would report zero
+	// evidence beside the t statistics' infinity.
+	rssScaleW, rssUnitW := 1.0, rssW
+	if rssW == 0 && maxResW > 0 {
+		rssScaleW = maxResW
+		rssUnitW = 0.0
+		for r := range n {
+			var wr float64
+			if fw != nil {
+				wr = fw[r]
+			} else {
+				wr = w.FloatAt(r)
+			}
+			d := out.Residuals[r] / maxResW
+			rssUnitW += wr * d * d
+		}
+	}
+	tssScaleW, tssUnitW := 1.0, tssW
+	if tssW == 0 && maxDevW > 0 {
+		tssScaleW = maxDevW
+		tssUnitW = 0.0
+		if hasConstant {
+			meanW := sumWY / sumW
+			for r := range n {
+				var wr, yv float64
+				if fw != nil {
+					wr = fw[r]
+				} else {
+					wr = w.FloatAt(r)
+				}
+				if fy != nil {
+					yv = fy[r]
+				} else {
+					yv = y.FloatAt(r)
+				}
+				d := (yv - meanW) / maxDevW
+				tssUnitW += wr * d * d
+			}
+		} else {
+			for r := range n {
+				var wr, yv float64
+				if fw != nil {
+					wr = fw[r]
+				} else {
+					wr = w.FloatAt(r)
+				}
+				if fy != nil {
+					yv = fy[r]
+				} else {
+					yv = y.FloatAt(r)
+				}
+				d := yv / maxDevW
+				tssUnitW += wr * d * d
+			}
 		}
 	}
 	tssDOF := n - 1
 	if !hasConstant {
 		tssDOF = n
 	}
-	if tssW == 0 {
+	if tssW == 0 && tssScaleW == 1 {
 		// Constant weighted response, exact fit: 1, as above.
 		out.RSquared = 1
 		out.AdjustedRSquared = 1
+	} else if tssW == 0 {
+		// The weighted total underflowed while the weighted response
+		// varies: the factored ratio, both R² measures rounding back
+		// to 1 while the F below reads the same pieces and does not.
+		ratio := rssUnitW / tssUnitW * (rssScaleW / tssScaleW) * (rssScaleW / tssScaleW)
+		out.RSquared = 1 - ratio
+		out.AdjustedRSquared = 1 - ratio*float64(tssDOF)/float64(out.DResidual)
 	} else {
 		out.RSquared = 1 - rssW/tssW
 		out.AdjustedRSquared = 1 - (rssW/float64(out.DResidual))/(tssW/float64(tssDOF))
@@ -557,6 +727,15 @@ func WeightedLinearRegression(x, y, w *core.Array) (*LinearRegressionResult, err
 			explained = 0 // rounding only, and a negative F is meaningless
 		}
 		out.FStatistic = explained / float64(out.DModel) / out.ResidualVariance
+		if (math.IsInf(out.FStatistic, 0) || math.IsNaN(out.FStatistic)) && (rssScaleW != 1 || tssScaleW != 1) {
+			// The factored F, as in the unweighted path: every scale
+			// factor divided out one step at a time.
+			out.FStatistic = (tssUnitW*tssScaleW/rssScaleW/rssScaleW*tssScaleW - rssUnitW) *
+				float64(out.DResidual) / (float64(out.DModel) * rssUnitW)
+			if out.FStatistic < 0 {
+				out.FStatistic = 0
+			}
+		}
 		if math.IsNaN(out.FStatistic) {
 			// 0/0, as in the unweighted path: nothing to test, p = 1.
 			out.FStatistic = 0

stats/regression_edge_pin_test.go

@@ -239,3 +239,112 @@ func mustMatrix(t *testing.T, vals []float64, r, c int) *core.Array {
 	}
 	return a
 }
+
+// TestLinearRegressionTinyScaleInference pins the inference of a response
+// on a scale whose squared residuals fall below the subnormal floor: the
+// plain residual and total sums of squares read zero there, and the fit
+// used to report the evidence backwards, R² of 1 with an infinite t and
+// p = 0 beside an F of zero with p = 1. The response y = [0, 0, e] over
+// x = 1, 2, 3 keeps every least-squares quantity exactly representable
+// while both sums of squares underflow: slope e/2, residual sum e²/6,
+// total 2e²/3, (XᵀX)⁻¹₁₁ = 1/2, so SE(slope) = e/(2√3), t = √3,
+// p = 1/3, F = 3 and R² = 3/4, all closed fractions.
+func TestLinearRegressionTinyScaleInference(t *testing.T) {
+	const e = 1e-200
+	x := mustMatrix(t, []float64{1, 1, 1, 2, 1, 3}, 3, 2)
+	y := mustFloats(t, []float64{0, 0, e}, 3)
+	res, err := LinearRegression(x, y)
+	if err != nil {
+		t.Fatalf("LinearRegression: %v", err)
+	}
+	if math.Abs(res.Coefficients[1]-e/2) > 1e-12*e/2 {
+		t.Fatalf("slope = %.17g, want %.17g", res.Coefficients[1], e/2)
+	}
+	wantSE := e / (2 * math.Sqrt(3))
+	se := res.StandardErrors[1]
+	if !(se > 0) || math.IsInf(se, 0) {
+		t.Fatalf("slope standard error = %g beside nonzero residuals, want %.17g", se, wantSE)
+	}
+	if math.Abs(se-wantSE) > 1e-12*wantSE {
+		t.Fatalf("slope standard error = %.17g, want %.17g", se, wantSE)
+	}
+	if math.Abs(res.TStatistics[1]-math.Sqrt(3)) > 1e-12 {
+		t.Fatalf("t = %.17g, want √3", res.TStatistics[1])
+	}
+	if math.Abs(res.PValues[1]-1.0/3) > 1e-12 {
+		t.Fatalf("p = %.17g, want 1/3", res.PValues[1])
+	}
+	if math.Abs(res.RSquared-0.75) > 1e-12 {
+		t.Fatalf("R² = %.17g, want 3/4", res.RSquared)
+	}
+	if math.Abs(res.AdjustedRSquared-0.5) > 1e-12 {
+		t.Fatalf("adjusted R² = %.17g, want 1/2", res.AdjustedRSquared)
+	}
+	if math.Abs(res.FStatistic-3) > 1e-11 {
+		t.Fatalf("F = %.17g, want 3", res.FStatistic)
+	}
+	if math.Abs(res.FPValue-1.0/3) > 1e-11 {
+		t.Fatalf("F p-value = %.17g, want 1/3", res.FPValue)
+	}
+
+	// Unit weights are the same fit, weighted statistics included.
+	w := mustFloats(t, []float64{1, 1, 1}, 3)
+	wres, err := WeightedLinearRegression(x, y, w)
+	if err != nil {
+		t.Fatalf("WeightedLinearRegression: %v", err)
+	}
+	if math.Abs(wres.TStatistics[1]-math.Sqrt(3)) > 1e-12 {
+		t.Fatalf("weighted t = %.17g, want √3", wres.TStatistics[1])
+	}
+	if math.Abs(wres.RSquared-0.75) > 1e-12 {
+		t.Fatalf("weighted R² = %.17g, want 3/4", wres.RSquared)
+	}
+	if math.Abs(wres.FStatistic-3) > 1e-11 {
+		t.Fatalf("weighted F = %.17g, want 3", wres.FStatistic)
+	}
+	if math.Abs(wres.FPValue-1.0/3) > 1e-11 {
+		t.Fatalf("weighted F p-value = %.17g, want 1/3", wres.FPValue)
+	}
+
+	// A response with a live scale beside the tiny spread keeps the same
+	// behaviour: the slope's own rounding leaves residuals near its last
+	// ulp, and the standard error must stay representable and the t
+	// finite rather than answer an exact fit the residuals contradict.
+	const a = 1e-160
+	step := math.Nextafter(3*a, math.Inf(1)) - 3*a
+	y2 := mustFloats(t, []float64{a, 2 * a, 3*a + step}, 3)
+	res2, err := LinearRegression(x, y2)
+	if err != nil {
+		t.Fatalf("LinearRegression: %v", err)
+	}
+	if se2 := res2.StandardErrors[1]; !(se2 > 0) || math.IsInf(se2, 0) {
+		t.Fatalf("slope standard error = %g beside nonzero residuals, want a representable value", se2)
+	}
+	if math.IsInf(res2.TStatistics[1], 0) {
+		t.Fatalf("t = %g beside nonzero residuals, want a finite statistic", res2.TStatistics[1])
+	}
+	if res2.FPValue > 1e-10 || res2.PValues[1] > 1e-10 {
+		t.Fatalf("p = %g, F p = %g, want the far tail both", res2.PValues[1], res2.FPValue)
+	}
+}
+
+// TestLinearRegressionTinyScaleExactLine pins the fully underflowed
+// corner: an exact line at a scale where both the residual and the total
+// sums of squares fall below the subnormal floor. The t statistics
+// already answer the exact fit with infinite evidence; the F test must
+// agree with them instead of reporting zero evidence.
+func TestLinearRegressionTinyScaleExactLine(t *testing.T) {
+	const a = 1e-300
+	x := mustMatrix(t, []float64{1, 1, 1, 2, 1, 3}, 3, 2)
+	y := mustFloats(t, []float64{a, 2 * a, 3 * a}, 3)
+	res, err := LinearRegression(x, y)
+	if err != nil {
+		t.Fatalf("LinearRegression: %v", err)
+	}
+	if res.FPValue > 1e-10 {
+		t.Fatalf("F p-value = %g on an exact line at a tiny scale, want the far tail beside the infinite t", res.FPValue)
+	}
+	if res.PValues[1] > 1e-10 {
+		t.Fatalf("slope p-value = %g, want the exact-fit report", res.PValues[1])
+	}
+}