// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

package linalg

import (
	"math"
	"math/cmplx"
	"os"
	"strings"
	"testing"

	"sourcedock.dev/petrbalvin/tensor/internal/core"
)

// Regression pins at the extremes of the float64 range: the dense and
// sparse decompositions, solvers and eigensolvers must keep their
// contracts at magnitudes from 1e-300 to 1e300, where a raw square
// overflows or vanishes.

// scales sweeps the magnitudes the decompositions have to survive:
// ordinary 1 as the control that pins the untouched arithmetic, the two
// edges of the squared-arithmetic window (a raw square leaves the normal
// range at about 1.3e154 and 1.5e-162), and the extremes of the float64
// range in both directions.
var scales = []float64{1, 1e150, 1e155, 1e200, 1e300, 1e-150, 1e-160, 1e-200, 1e-300}

// scale multiplies every entry of a flat matrix by s.
func scale(vals []float64, s float64) []float64 {
	out := make([]float64, len(vals))
	for i, v := range vals {
		out[i] = v * s
	}
	return out
}

// scaleC multiplies every complex entry by the real s.
func scaleC(vals []complex128, s float64) []complex128 {
	out := make([]complex128, len(vals))
	for i, v := range vals {
		out[i] = complex(real(v)*s, imag(v)*s)
	}
	return out
}

// iJ is the 4x4 matrix I+J: 2 on the diagonal, 1 elsewhere. Its
// spectrum is closed form, 1 three times and 5 once, which is the
// reference every 1e+-200 eigen test below checks against.
func iJ(n int) []float64 {
	out := make([]float64, n*n)
	for i := range n {
		for j := range n {
			if i == j {
				out[i*n+j] = 2
			} else {
				out[i*n+j] = 1
			}
		}
	}
	return out
}

// pinDenseSolve solves a·x = b by Gaussian elimination with partial
// pivoting. It is the independent brute-force reference the sparse
// solvers are checked against, deliberately written without any of the
// library's own machinery.
func pinDenseSolve(t *testing.T, a, b []complex128, n int) []complex128 {
	t.Helper()
	m := append([]complex128(nil), a...)
	x := append([]complex128(nil), b...)
	for k := range n {
		piv := k
		for i := k + 1; i < n; i++ {
			if cmplx.Abs(m[i*n+k]) > cmplx.Abs(m[piv*n+k]) {
				piv = i
			}
		}
		if m[piv*n+k] == 0 {
			t.Fatalf("dense reference: singular matrix at column %d", k)
		}
		if piv != k {
			for j := range n {
				m[k*n+j], m[piv*n+j] = m[piv*n+j], m[k*n+j]
			}
			x[k], x[piv] = x[piv], x[k]
		}
		for i := k + 1; i < n; i++ {
			f := m[i*n+k] / m[k*n+k]
			for j := k; j < n; j++ {
				m[i*n+j] -= f * m[k*n+j]
			}
			x[i] -= f * x[k]
		}
	}
	for k := n - 1; k >= 0; k-- {
		for j := k + 1; j < n; j++ {
			x[k] -= m[k*n+j] * x[j]
		}
		x[k] /= m[k*n+k]
	}
	return x
}

// hermitianStencil builds an n×n Hermitian tridiagonal stencil:
// 2 on the diagonal and conjugate mirrored imaginary couplings. It is
// the mirrored-stencil shape the complex sparse solvers are exercised on.
func hermitianStencil(n int) []complex128 {
	a := make([]complex128, n*n)
	for i := range n {
		a[i*n+i] = 2
		if i+1 < n {
			a[i*n+i+1] = complex(0, 0.5)
			a[(i+1)*n+i] = complex(0, -0.5)
		}
	}
	return a
}

// cSRFromDense builds a complex SparseCOO from a dense flat matrix,
// dropping the exact zeros the way a caller's assembly would.
func cSRFromDense(t *testing.T, a []complex128, n int) *core.SparseCOO {
	t.Helper()
	var idx []int64
	var vals []complex128
	for i := range n {
		for j := range n {
			if a[i*n+j] == 0 {
				continue
			}
			idx = append(idx, int64(i), int64(j))
			vals = append(vals, a[i*n+j])
		}
	}
	ind, err := core.FromInts(idx, len(vals), 2)
	if err != nil {
		t.Fatalf("FromInts: %v", err)
	}
	val, err := core.FromComplexes(vals, len(vals))
	if err != nil {
		t.Fatalf("FromComplexes: %v", err)
	}
	coo, err := core.NewSparseCOO(ind, val, []int{n, n})
	if err != nil {
		t.Fatalf("NewSparseCOO: %v", err)
	}
	return coo
}

// TestHouseholderVectorReflectsAtExtremeScale pins the reflector itself
// (report T1/F7): for a non-zero x, beta must be positive and finite at
// every magnitude, not +0 (which callers read as "no reflection") or
// +Inf (whose product with a zero dot is NaN), and H must still map x to
// -sign(x0)·‖x‖·e1. The reflection is applied to x/s, so the check's own
// arithmetic stays in range at every scale.
func TestHouseholderVectorReflectsAtExtremeScale(t *testing.T) {
	for _, s := range scales {
		for _, x0 := range []float64{s, -s} {
			x := []float64{x0, s, s}
			dst := make([]float64, 3)
			hh := householderVectorInto(dst, x)
			if !(hh.beta > 0) || math.IsInf(hh.beta, 0) {
				t.Fatalf("scale %g, x0 %g: beta = %v, want a positive finite reflector", s, x0, hh.beta)
			}
			for i, v := range hh.v {
				if math.IsInf(v, 0) || math.IsNaN(v) {
					t.Fatalf("scale %g, x0 %g: v[%d] = %v, want finite", s, x0, i, v)
				}
			}
			xs := []float64{x0 / s, 1, 1}
			dot := 0.0
			for i := range xs {
				dot += hh.v[i] * xs[i]
			}
			w := hh.beta * dot
			sign := 1.0
			if x0 < 0 {
				sign = -1
			}
			for i := range xs {
				want := 0.0
				if i == 0 {
					want = -sign * math.Sqrt(3)
				}
				if got := xs[i] - hh.v[i]*w; math.Abs(got-want) > 1e-13 {
					t.Fatalf("scale %g, x0 %g: (H x)[%d] = %g, want %g", s, x0, i, got, want)
				}
			}
		}
		// The zero vector still answers the identity reflector: beta == 0
		// is the legitimate "nothing to reflect" signal and must stay
		// distinguishable from the overflow the fix removes.
		if hh := householderVectorInto(make([]float64, 3), []float64{0, 0, 0}); hh.beta != 0 {
			t.Fatalf("scale %g: zero vector beta = %v, want 0", s, hh.beta)
		}
	}
}

// TestEigenSpectrumAtExtremeScale pins the closed-form spectrum of
// (I+J)·s: 1,1,1,5 times s, with the eigenvectors orthonormal and the
// residual at rounding level. Before the fix, 1e155 returned the raw
// diagonal {2,2,2,2}e155 (the reflector's beta came out as +0 and was
// skipped) and 1e-200 returned all NaN (beta came out as +Inf).
func TestEigenSpectrumAtExtremeScale(t *testing.T) {
	const n = 4
	want := []float64{1, 1, 1, 5}
	for _, s := range scales {
		a := mustFloats(t, scale(iJ(n), s), n, n)
		vals, vecs, err := Eigen(a)
		if err != nil {
			t.Fatalf("Eigen at scale %g: %v", s, err)
		}
		for i := range n {
			if got := vals.FloatAt(i) / s; math.Abs(got-want[i]) > 1e-12 {
				t.Fatalf("scale %g: eigenvalue %d = %g, want %g", s, i, got, want[i])
			}
		}
		// Residual of the eigenpair on the O(1) shift: ‖(I+J)v - (λ/s)v‖,
		// which is the residual of the original problem divided by s.
		base := iJ(n)
		for k := range n {
			lam := vals.FloatAt(k) / s
			worst := 0.0
			for i := range n {
				acc := 0.0
				for j := range n {
					acc += base[i*n+j] * vecs.FloatAt(j*n+k)
				}
				worst = math.Max(worst, math.Abs(acc-lam*vecs.FloatAt(i*n+k)))
			}
			if worst > 1e-12 {
				t.Fatalf("scale %g: residual of eigenpair %d = %g, want <= 1e-12", s, k, worst)
			}
		}
		// VᵀV = I: a scaling mistake that dropped the factor would still
		// leave orthonormal columns, so this is a cheap second net.
		for j := range n {
			for k := range n {
				acc := 0.0
				for i := range n {
					acc += vecs.FloatAt(i*n+j) * vecs.FloatAt(i*n+k)
				}
				want := 0.0
				if j == k {
					want = 1
				}
				if math.Abs(acc-want) > 1e-12 {
					t.Fatalf("scale %g: (VᵀV)[%d,%d] = %g, want %g", s, j, k, acc, want)
				}
			}
		}
	}
}

// TestSVDReconstructsAtExtremeScale pins A = U·Σ·Vᵀ against the
// hand-built 2×2 ([[0,s],[s,0]] has both singular values s) and against
// a dense 4x4 with a closed-form scale-free reconstruction. Before the
// fix, 1e155 returned sigma = (+Inf, +Inf) and the dense case all NaN.
func TestSVDReconstructsAtExtremeScale(t *testing.T) {
	for _, s := range scales {
		a := mustFloats(t, []float64{0, s, s, 0}, 2, 2)
		u, sigma, vt, err := SVD(a)
		if err != nil {
			t.Fatalf("SVD at scale %g: %v", s, err)
		}
		for i := range 2 {
			if got := sigma.FloatAt(i) / s; math.Abs(got-1) > 1e-13 {
				t.Fatalf("scale %g: sigma[%d] = %g, want 1 (in units of s)", s, i, got)
			}
		}
		// [[0,1],[1,0]] = U·(Σ/s)·Vᵀ in scaled units.
		for i := range 2 {
			for j := range 2 {
				acc := 0.0
				for k := range 2 {
					acc += u.FloatAt(i*2+k) * (sigma.FloatAt(k) / s) * vt.FloatAt(k*2+j)
				}
				want := 0.0
				if i != j {
					want = 1
				}
				if math.Abs(acc-want) > 1e-13 {
					t.Fatalf("scale %g: (UΣVᵀ)[%d,%d] = %g, want %g", s, i, j, acc, want)
				}
			}
		}
	}
	// The dense case the report saw as all NaN at 1e155.
	const n = 4
	base := []float64{
		4, 1, 2, 0.5,
		1, 3, 0.25, 1,
		2, 0.25, 5, 1,
		0.5, 1, 1, 2,
	}
	for _, s := range []float64{1, 1e155, 1e-155, 1e200} {
		a := mustFloats(t, scale(base, s), n, n)
		u, sigma, vt, err := SVD(a)
		if err != nil {
			t.Fatalf("SVD dense at scale %g: %v", s, err)
		}
		worst := 0.0
		for i := range n {
			for j := range n {
				// A/s = U·(Σ/s)·Vᵀ.
				acc := 0.0
				for k := range n {
					acc += u.FloatAt(i*n+k) * (sigma.FloatAt(k) / s) * vt.FloatAt(k*n+j)
				}
				worst = math.Max(worst, math.Abs(acc-base[i*n+j]))
			}
		}
		if worst > 1e-13 {
			t.Fatalf("scale %g: dense reconstruction error %g, want <= 1e-13", s, worst)
		}
	}
}

// TestEigenGeneralSpectrumAtExtremeScale pins EigenGeneral on (I+J)·s,
// where the spectrum is closed form and the matrix is far from
// defective. Before the fix, 1e-200 silently returned {2,2,2,2}e-200 and
// 1e160 failed with a spurious "QR iteration failed to converge".
func TestEigenGeneralSpectrumAtExtremeScale(t *testing.T) {
	const n = 4
	want := []float64{1, 1, 1, 5}
	for _, s := range scales {
		a := mustFloats(t, scale(iJ(n), s), n, n)
		vals, vecs, err := EigenGeneral(a)
		if err != nil {
			t.Fatalf("EigenGeneral at scale %g: %v", s, err)
		}
		// Values come back descending by magnitude: 5 then 1,1,1.
		got := make([]float64, n)
		for i := range n {
			got[i] = cmplx.Abs(vals.ComplexAt(i)) / s
		}
		if math.Abs(got[0]-5) > 1e-12 {
			t.Fatalf("scale %g: |λ| max = %g, want 5", s, got[0])
		}
		for i := 1; i < n; i++ {
			if math.Abs(got[i]-1) > 1e-12 {
				t.Fatalf("scale %g: |λ| %d = %g, want %g", s, i, got[i], want[i])
			}
		}
		// Residual of every eigenpair on the O(1) shift.
		base := iJ(n)
		for k := range n {
			lam := vals.ComplexAt(k) / complex(s, 0)
			worst := 0.0
			for i := range n {
				acc := complex(0, 0)
				for j := range n {
					acc += complex(base[i*n+j], 0) * vecs.ComplexAt(j*n+k)
				}
				worst = math.Max(worst, cmplx.Abs(acc-lam*vecs.ComplexAt(i*n+k)))
			}
			if worst > 1e-12 {
				t.Fatalf("scale %g: residual of eigenpair %d = %g, want <= 1e-12", s, k, worst)
			}
		}
	}
	// A mixed-scale matrix: a huge diagonal with couplings 1e-260 of it.
	// The k = 1 reflector's column is then about 1e-200 while the matrix
	// max is 1e60, so vMax*vMax underflows to zero, beta comes out as
	// +Inf and the update turns the reduction into NaN, even though the
	// matrix itself is inside the safe window. The spectrum is that of
	// the huge diagonal to rounding.
	const big, rel = 1e60, 1e-260
	mixed := []float64{
		big, 0, 0, 0,
		0, big, big * rel, big * rel,
		0, big * rel, big, big * rel,
		0, big * rel, big * rel, big,
	}
	valsM, _, err := EigenGeneral(mustFloats(t, mixed, n, n))
	if err != nil {
		t.Fatalf("EigenGeneral of the mixed-scale matrix: %v", err)
	}
	for i := range n {
		got := valsM.ComplexAt(i)
		if math.IsNaN(real(got)) || math.IsNaN(imag(got)) {
			t.Fatalf("mixed-scale eigenvalue %d = %v, want a finite value", i, got)
		}
		if scale := cmplx.Abs(got) / big; math.Abs(scale-1) > 1e-12 {
			t.Fatalf("mixed-scale eigenvalue %d = %v, want magnitude %g", i, got, big)
		}
	}
	// The reflector-side face of the same defect, reached with an
	// ordinary matrix: a unit diagonal with a column at 1e-160. That
	// column is above the magnitude floor, so the reflector is built, but
	// its raw squares are subnormal, beta overflows to +Inf and the whole
	// reduction comes back NaN where the matrix is perfectly valid. Only
	// the reflector's rescaling survives it, so this case pins that
	// mechanism on its own.
	const tiny = 1e-160
	near := []float64{
		1, 0, 0, 0,
		0, 1, tiny, tiny,
		0, tiny, 1, tiny,
		0, tiny, tiny, 1,
	}
	valsN, _, err := EigenGeneral(mustFloats(t, near, n, n))
	if err != nil {
		t.Fatalf("EigenGeneral of the subnormal-column matrix: %v", err)
	}
	for i := range n {
		got := valsN.ComplexAt(i)
		if math.IsNaN(real(got)) || math.IsNaN(imag(got)) {
			t.Fatalf("subnormal-column eigenvalue %d = %v, want a finite value", i, got)
		}
		if math.Abs(cmplx.Abs(got)-1) > 1e-12 {
			t.Fatalf("subnormal-column eigenvalue %d = %v, want magnitude 1", i, got)
		}
	}
}

// TestSchurComplexAndMatrixSqrtAtExtremeScale pins the Schur contract
// A = Q·T·Qᴴ with T upper triangular and the diagonal of T the closed-form
// circulant spectrum, plus the defining property of the principal square
// root, r·r = A. Before the fix, SchurComplex, MatrixSqrt and MatrixLog
// all failed with "the Schur iteration did not converge" at 1e-200 and
// 1e160 on perfectly valid input.
func TestSchurComplexAndMatrixSqrtAtExtremeScale(t *testing.T) {
	const n = 4
	base := []float64{
		4, 1, 2, 3,
		3, 4, 1, 2,
		2, 3, 4, 1,
		1, 2, 3, 4,
	}
	// Eigenvalues of the symmetric circulant with first row [4,1,2,3]:
	// 10, 2, 2+2i, 2-2i, computed from the closed form Σ c_j ω^{jk}.
	want := []complex128{10, 2, complex(2, 2), complex(2, -2)}
	for _, s := range scales {
		a := mustFloats(t, scale(base, s), n, n)
		tm, q, err := SchurComplex(a)
		if err != nil {
			t.Fatalf("SchurComplex at scale %g: %v", s, err)
		}
		// A/s = Q·(T/s)·Qᴴ.
		worst := 0.0
		for i := range n {
			for j := range n {
				acc := complex(0, 0)
				for k := range n {
					for l := range n {
						acc += q.ComplexAt(i*n+k) * (tm.ComplexAt(k*n+l) / complex(s, 0)) * cmplx.Conj(q.ComplexAt(j*n+l))
					}
				}
				worst = math.Max(worst, cmplx.Abs(acc-complex(base[i*n+j], 0)))
			}
		}
		if worst > 1e-13 {
			t.Fatalf("scale %g: Schur reconstruction error %g, want <= 1e-13", s, worst)
		}
		// T upper triangular, and its diagonal the closed-form spectrum.
		for i := 1; i < n; i++ {
			for j := 0; j < i; j++ {
				if v := cmplx.Abs(tm.ComplexAt(i*n+j)) / s; v > 1e-13 {
					t.Fatalf("scale %g: T[%d,%d] = %g, want upper triangular", s, i, j, v)
				}
			}
		}
		got := make([]complex128, n)
		for i := range n {
			got[i] = tm.ComplexAt(i*n+i) / complex(s, 0)
		}
		if !spectrumMatches(got, want, 1e-12) {
			t.Fatalf("scale %g: Schur diagonal %v, want %v", s, got, want)
		}
		// The principal square root of [[7,10],[15,22]]·s satisfies
		// r·r = A; the closed form of the O(1) matrix is [[1,2],[3,4]]²,
		// so the check is (r/√s)² = [[7,10],[15,22]].
		ms, err := MatrixSqrt(mustFloats(t, scale([]float64{7, 10, 15, 22}, s), 2, 2))
		if err != nil {
			t.Fatalf("MatrixSqrt at scale %g: %v", s, err)
		}
		r := math.Sqrt(s)
		rr := []float64{
			ms.FloatAt(0)*ms.FloatAt(0) + ms.FloatAt(1)*ms.FloatAt(2),
			ms.FloatAt(0)*ms.FloatAt(1) + ms.FloatAt(1)*ms.FloatAt(3),
			ms.FloatAt(2)*ms.FloatAt(0) + ms.FloatAt(3)*ms.FloatAt(2),
			ms.FloatAt(2)*ms.FloatAt(1) + ms.FloatAt(3)*ms.FloatAt(3),
		}
		sq := []float64{7, 10, 15, 22}
		for i := range 4 {
			if got := rr[i] / (r * r); math.Abs(got-sq[i]) > 1e-12 {
				t.Fatalf("scale %g: (r·r)[%d] = %g, want %g", s, i, got, sq[i])
			}
		}
	}
}

// spectrumMatches reports whether every want finds a distinct got
// within tol: a spectrum is a set, and the order of near-equal complex
// values is not part of the contract.
func spectrumMatches(got, want []complex128, tol float64) bool {
	used := make([]bool, len(got))
	for _, w := range want {
		found := false
		for i, g := range got {
			if !used[i] && cmplx.Abs(g-w) <= tol {
				used[i] = true
				found = true
				break
			}
		}
		if !found {
			return false
		}
	}
	return true
}

// TestSVDComplexReconstructsAtExtremeScale pins A = U·Σ·Vᴴ on the
// hand-built [[0,t],[t,0]], whose singular values are both t, and on
// diag(1e160, 1e-160). Before the fix, t = 1e-200 silently returned
// sigma = (0,0) (the reflector norm underflowed to zero and the
// below-diagonal mass was zeroed with it) and t = 1e160 returned NaN.
func TestSVDComplexReconstructsAtExtremeScale(t *testing.T) {
	for _, t0 := range scales {
		a := mustFromComplexes(t, []complex128{
			0, complex(t0, 0),
			complex(t0, 0), 0,
		}, 2, 2)
		u, sigma, vh, err := SVDComplex(a)
		if err != nil {
			t.Fatalf("SVDComplex at scale %g: %v", t0, err)
		}
		for i := range 2 {
			if got := sigma.FloatAt(i) / t0; math.Abs(got-1) > 1e-13 {
				t.Fatalf("scale %g: sigma[%d] = %g, want 1 (in units of t)", t0, i, got)
			}
		}
		for i := range 2 {
			for j := range 2 {
				acc := complex(0, 0)
				for k := range 2 {
					acc += u.ComplexAt(i*2+k) * complex(sigma.FloatAt(k)/t0, 0) * vh.ComplexAt(k*2+j)
				}
				want := complex(0, 0)
				if i != j {
					want = 1
				}
				if d := cmplx.Abs(acc - want); d > 1e-13 {
					t.Fatalf("scale %g: (UΣVᴴ)[%d,%d] error %g, want 0", t0, i, j, d)
				}
			}
		}
	}
	// The mixed-spectrum case the report saw as NaN, NaN.
	a := mustFromComplexes(t, []complex128{complex(1e160, 0), 0, 0, complex(1e-160, 0)}, 2, 2)
	_, sigma, _, err := SVDComplex(a)
	if err != nil {
		t.Fatalf("SVDComplex diag(1e160, 1e-160): %v", err)
	}
	if got := sigma.FloatAt(0) / 1e160; math.Abs(got-1) > 1e-13 {
		t.Fatalf("diag(1e160, 1e-160): sigma[0] = %g, want 1e160", sigma.FloatAt(0))
	}
	if got := sigma.FloatAt(1) / 1e-160; math.Abs(got-1) > 1e-13 {
		t.Fatalf("diag(1e160, 1e-160): sigma[1] = %g, want 1e-160", sigma.FloatAt(1))
	}
}

// TestEigenComplexSpectrumAtExtremeScale pins the closed-form spectrum of
// the Hermitian [[2s,is],[-is,2s]], which is {s, 3s}. Before the fix the
// Jacobi convergence test's raw sum of squares overflowed to +Inf at
// 1e200 (every sweep looked converged) and underflowed to 0 at 1e-200,
// so both returned the unrotated diagonal {2s, 2s} silently.
func TestEigenComplexSpectrumAtExtremeScale(t *testing.T) {
	for _, s := range scales {
		a := mustFromComplexes(t, []complex128{
			complex(2*s, 0), complex(0, s),
			complex(0, -s), complex(2*s, 0),
		}, 2, 2)
		vals, vecs, err := EigenComplex(a)
		if err != nil {
			t.Fatalf("EigenComplex at scale %g: %v", s, err)
		}
		if got := vals.FloatAt(0) / s; math.Abs(got-1) > 1e-13 {
			t.Fatalf("scale %g: eigenvalue 0 = %g, want 1 (in units of s)", s, got)
		}
		if got := vals.FloatAt(1) / s; math.Abs(got-3) > 1e-13 {
			t.Fatalf("scale %g: eigenvalue 1 = %g, want 3 (in units of s)", s, got)
		}
		// Residual of both eigenpairs on the O(1) shift [[2,i],[-i,2]].
		base := []complex128{2, complex(0, 1), complex(0, -1), 2}
		for k := range 2 {
			lam := complex(vals.FloatAt(k)/s, 0)
			worst := 0.0
			for i := range 2 {
				acc := complex(0, 0)
				for j := range 2 {
					acc += base[i*2+j] * vecs.ComplexAt(j*2+k)
				}
				worst = math.Max(worst, cmplx.Abs(acc-lam*vecs.ComplexAt(i*2+k)))
			}
			if worst > 1e-13 {
				t.Fatalf("scale %g: residual of eigenpair %d = %g, want <= 1e-13", s, k, worst)
			}
		}
	}
}

// TestSpEigenComplexAtExtremeScale pins the Hermitian Lanczos on the
// mirrored stencil against the closed-form spectrum of the ordinary
// matrix: 2 + cos(k·π/(n+1)) for k = 1..n, whose two largest values are
// the top Ritz pairs. Before the scaling the recurrence collapsed at
// 1e200 (the projected tridiagonal's squared magnitudes overflowed) and
// returned NaN, and at 1e-200 the values came back as zero.
func TestSpEigenComplexAtExtremeScale(t *testing.T) {
	const n = 8
	want := []float64{
		2 + math.Cos(math.Pi/float64(n+1)),
		2 + math.Cos(2*math.Pi/float64(n+1)),
	}
	base := hermitianStencil(n)
	for _, s := range scales {
		coo := cSRFromDense(t, scaleC(base, s), n)
		vals, vecs, err := SpEigenComplex(coo, 2, core.NewGenerator(7))
		if err != nil {
			t.Fatalf("SpEigenComplex at scale %g: %v", s, err)
		}
		for i := range 2 {
			got := vals.FloatAt(i)
			if math.IsNaN(got) {
				t.Fatalf("scale %g: Ritz value %d = NaN", s, i)
			}
			if rel := got / s; math.Abs(rel-want[i]) > 1e-12 {
				t.Fatalf("scale %g: Ritz value %d = %g, want %g (in units of s)", s, i, rel, want[i])
			}
		}
		// The Ritz vectors solve the O(1) stencil to rounding.
		for k := range 2 {
			lam := complex(vals.FloatAt(k)/s, 0)
			worst := 0.0
			for i := range n {
				acc := complex(0, 0)
				for j := range n {
					acc += base[i*n+j] * vecs.ComplexAt(j*2+k)
				}
				worst = math.Max(worst, cmplx.Abs(acc-lam*vecs.ComplexAt(i*2+k)))
			}
			if worst > 1e-12 {
				t.Fatalf("scale %g: Ritz residual %d = %g, want <= 1e-12", s, k, worst)
			}
		}
	}
}

// TestSpSolveComplexAtExtremeScale pins both complex sparse solvers
// against a brute-force dense elimination of the mirrored-stencil system,
// at ordinary scale (the control the oracle pins) and at the extremes.
// Before the fix, a huge b made bNorm +Inf so every convergence check
// passed and both solvers returned after one step (residual 0.5s), and a
// tiny b made bNorm 0 so the zero vector was returned as the exact
// solution.
func TestSpSolveComplexAtExtremeScale(t *testing.T) {
	const n = 8
	// Reference: the ordinary-scale system solved densely.
	base := hermitianStencil(n)
	rhs := make([]complex128, n)
	for i := range n {
		rhs[i] = complex(1, 0.25)
	}
	want := pinDenseSolve(t, base, rhs, n)

	for _, s := range scales {
		coo := cSRFromDense(t, scaleC(base, s), n)
		b := mustFromComplexes(t, scaleC(rhs, s), n)
		x, err := SpSolveComplexCG(coo, b, 1e-12, 400)
		if err != nil {
			t.Fatalf("SpSolveComplexCG at scale %g: %v", s, err)
		}
		// (sA)x = s·b has exactly the ordinary-scale solution, so the
		// returned x must equal the dense reference as it stands: the
		// factor cancels and must not be applied again.
		worst := 0.0
		for i := range n {
			worst = math.Max(worst, cmplx.Abs(x.ComplexAt(i)-want[i]))
		}
		if worst > 1e-9 {
			t.Fatalf("scale %g: CG solution differs from the dense reference by %g", s, worst)
		}
		// The same system through BiCGSTAB (Hermitian input is a valid
		// special case of its contract).
		coo2 := cSRFromDense(t, scaleC(base, s), n)
		xb, err := SpSolveComplexBiCGSTAB(coo2, b, 1e-12, 400)
		if err != nil {
			t.Fatalf("SpSolveComplexBiCGSTAB at scale %g: %v", s, err)
		}
		worst = 0.0
		for i := range n {
			worst = math.Max(worst, cmplx.Abs(xb.ComplexAt(i)-want[i]))
		}
		if worst > 1e-9 {
			t.Fatalf("scale %g: BiCGSTAB solution differs from the dense reference by %g", s, worst)
		}
	}

	// The non-Hermitian shape BiCGSTAB exists for, same reference method.
	nonHerm := make([]complex128, n*n)
	for i := range n {
		nonHerm[i*n+i] = complex(2, 1)
		if i+1 < n {
			nonHerm[i*n+i+1] = complex(1, 0.3)
			nonHerm[(i+1)*n+i] = complex(0, -0.25)
		}
	}
	wantNH := pinDenseSolve(t, nonHerm, rhs, n)
	for _, s := range scales {
		coo := cSRFromDense(t, scaleC(nonHerm, s), n)
		b := mustFromComplexes(t, scaleC(rhs, s), n)
		x, err := SpSolveComplexBiCGSTAB(coo, b, 1e-12, 400)
		if err != nil {
			t.Fatalf("non-Hermitian BiCGSTAB at scale %g: %v", s, err)
		}
		worst := 0.0
		for i := range n {
			worst = math.Max(worst, cmplx.Abs(x.ComplexAt(i)-wantNH[i]))
		}
		if worst > 1e-9 {
			t.Fatalf("non-Hermitian scale %g: solution differs from the dense reference by %g", s, worst)
		}
	}
}

// TestSpEigenComplexErrorNamesK pins the malformed diagnostic (report
// F14): the k-range error had three verbs and two arguments and rendered
// as "%!s(int=2): k must be in [1, 5], got %!d(MISSING)".
func TestSpEigenComplexErrorNamesK(t *testing.T) {
	coo := cSRFromDense(t, hermitianStencil(5), 5)
	_, _, err := SpEigenComplex(coo, 6, nil)
	if err == nil {
		t.Fatal("SpEigenComplex with k > n: want an error")
	}
	msg := err.Error()
	if strings.Contains(msg, "%!") {
		t.Fatalf("malformed error message: %q", msg)
	}
	if !strings.Contains(msg, "SpEigenComplex: k must be in [1, 5], got 6") {
		t.Fatalf("error message = %q, want it to name the entry point and the range", msg)
	}
	if _, _, err := SpEigenComplex(coo, 0, nil); err == nil || strings.Contains(err.Error(), "%!") {
		t.Fatalf("k = 0 error = %v, want a well-formed range error", err)
	}
}

// TestSymmetryGuardsAreRelativeAtSmallScale pins the guards of Eigen,
// EigenComplex and the complex sparse Hermitian check on a matrix whose
// asymmetry is a tenth of its own scale: an absolute floor of 1e-15
// approved such a matrix and the symmetric path then answered for a
// matrix that is neither A nor Aᵀ, while the wording of the refusal is
// unchanged. The exactly symmetric matrix of the same scale must still
// be accepted, so the rule is relative and not simply stricter.
func TestSymmetryGuardsAreRelativeAtSmallScale(t *testing.T) {
	const s = 1e-14
	// Relative asymmetry 1e-1.
	asym := mustFloats(t, []float64{s, 0, 1e-15, s}, 2, 2)
	if _, _, err := Eigen(asym); err == nil {
		t.Fatal("Eigen: asymmetric at 10% of scale 1e-14 was accepted, want a refusal")
	} else if !strings.Contains(err.Error(), "Eigen: matrix is not symmetric within 1e-12 tolerance") {
		t.Fatalf("Eigen refusal = %q, want the unchanged wording", err)
	}
	// A genuinely symmetric matrix of the same tiny scale still passes
	// the guard: if it did not, the test would pin a stricter rule than
	// the contract asks for (the tiny p.d. case must reach the solver).
	sym := mustFloats(t, []float64{s, 1e-15, 1e-15, s}, 2, 2)
	if _, _, err := Eigen(sym); err != nil {
		t.Fatalf("Eigen of a symmetric matrix at scale %g: %v", s, err)
	}

	nonHerm := mustFromComplexes(t, []complex128{
		complex(s, 0), complex(1e-15, 0),
		0, complex(s, 0),
	}, 2, 2)
	if _, _, err := EigenComplex(nonHerm); err == nil {
		t.Fatal("EigenComplex: non-Hermitian at 10% of scale 1e-14 was accepted, want a refusal")
	} else if !strings.Contains(err.Error(), "EigenComplex: matrix is not Hermitian within 1e-12 tolerance") {
		t.Fatalf("EigenComplex refusal = %q, want the unchanged wording", err)
	}
	herm := mustFromComplexes(t, []complex128{
		complex(s, 0), complex(0, 1e-15),
		complex(0, -1e-15), complex(s, 0),
	}, 2, 2)
	if _, _, err := EigenComplex(herm); err != nil {
		t.Fatalf("EigenComplex of a Hermitian matrix at scale %g: %v", s, err)
	}

	// The sparse complex twin: both mirrored entries are stored, and the
	// asymmetry between i·1e-16 and i·1e-17 is about a percent of the
	// 1e-14 scale yet below the 1e-15 absolute floor, so the pre-fix guard
	// approved it and Lanczos ran on a matrix that is not Hermitian. The
	// refusal wording is pinned as well.
	coo := cSRFromDense(t, []complex128{
		complex(s, 0), complex(0, 1e-16),
		complex(0, 1e-17), complex(s, 0),
	}, 2)
	if _, _, err := SpEigenComplex(coo, 1, core.NewGenerator(3)); err == nil {
		t.Fatal("SpEigenComplex: non-Hermitian at 10% of scale 1e-14 was accepted, want a refusal")
	} else if !strings.Contains(err.Error(), "SpEigenComplex: matrix is not Hermitian within 1e-12 tolerance") {
		t.Fatalf("SpEigenComplex refusal = %q, want the unchanged wording", err)
	}
	// The exactly Hermitian matrix of the same scale must still be
	// accepted: the rule is relative, not simply stricter.
	hermSparse := cSRFromDense(t, []complex128{
		complex(s, 0), complex(0, 1e-16),
		complex(0, -1e-16), complex(s, 0),
	}, 2)
	if _, _, err := SpEigenComplex(hermSparse, 1, core.NewGenerator(3)); err != nil {
		t.Fatalf("SpEigenComplex of a Hermitian matrix at scale %g: %v", s, err)
	}
}

// TestNoArrowSymbolsInSource pins the source convention: the arrow
// code points must not appear as assignment notation in comments,
// which the house style forbids.
func TestNoArrowSymbolsInSource(t *testing.T) {
	for _, name := range []string{
		"decomp2.go", "decomp3.go", "eigenreal.go",
		"schur.go", "csvd.go", "sparsecomplex.go",
	} {
		src, err := os.ReadFile(name)
		if err != nil {
			t.Fatalf("read %s: %v", name, err)
		}
		for i, line := range strings.Split(string(src), "\n") {
			for _, r := range []rune{'←', '→', '↑', '↓', '⇒'} {
				if strings.ContainsRune(line, r) {
					t.Fatalf("%s:%d contains U+%04X: %s", name, i+1, r, line)
				}
			}
		}
	}
}

// stencil builds the n×n tridiagonal 2 / 0.5 stencil. Its spectrum
// is closed form, 2 + cos(k·π/(n+1)) for k = 1..n, which is the reference
// the sparse eigensolvers below are checked against.
func stencil(n int) []float64 {
	out := make([]float64, n*n)
	for i := range n {
		out[i*n+i] = 2
		if i+1 < n {
			out[i*n+i+1] = 0.5
			out[(i+1)*n+i] = 0.5
		}
	}
	return out
}

// topStencilValues is the two largest eigenvalues of an n×n 2 / 0.5
// stencil, derived from that closed form.
func topStencilValues(n int) []float64 {
	return []float64{
		2 + math.Cos(math.Pi/float64(n+1)),
		2 + math.Cos(2*math.Pi/float64(n+1)),
	}
}

// realCOO builds a real SparseCOO from a dense flat matrix,
// dropping the exact zeros the way a caller's assembly would.
func realCOO(t *testing.T, a []float64, n int) *core.SparseCOO {
	t.Helper()
	var idx []int64
	var vals []float64
	for i := range n {
		for j := range n {
			if a[i*n+j] == 0 {
				continue
			}
			idx = append(idx, int64(i), int64(j))
			vals = append(vals, a[i*n+j])
		}
	}
	ind, err := core.FromInts(idx, len(vals), 2)
	if err != nil {
		t.Fatalf("FromInts: %v", err)
	}
	val, err := core.FromFloats(vals, len(vals))
	if err != nil {
		t.Fatalf("FromFloats: %v", err)
	}
	coo, err := core.NewSparseCOO(ind, val, []int{n, n})
	if err != nil {
		t.Fatalf("NewSparseCOO: %v", err)
	}
	return coo
}

// ritzResidual returns max_k ‖(A/s)·u_k − (λ_k/s)·u_k‖ for the k
// Ritz pairs of a sparse matrix held as a COO, with A/s formed from the
// stored values so nothing overflows at either extreme of the scale
// sweep. vals and vecs are the returned values and the (n, k) vectors.
func ritzResidual(a *core.SparseCOO, vals, vecs []complex128, n, k int, s float64) float64 {
	idx := a.Indices.RawInts()
	nnz := a.Indices.Shape()[0]
	worst := 0.0
	for kk := range k {
		lam := vals[kk] / complex(s, 0)
		for i := range n {
			acc := complex(0, 0)
			for p := range nnz {
				if int(idx[p*2]) != i {
					continue
				}
				v := complex(0, 0)
				if a.Values.Dtype() == core.Complex {
					v = a.Values.ComplexAt(p)
				} else {
					v = complex(a.Values.FloatAt(p), 0)
				}
				acc += (v / complex(s, 0)) * vecs[int(idx[p*2+1])*k+kk]
			}
			if d := cmplx.Abs(acc - lam*vecs[i*k+kk]); d > worst {
				worst = d
			}
		}
	}
	return worst
}

// ritzResidualReal is ritzResidual for a real Ritz pair set.
func ritzResidualReal(a *core.SparseCOO, vals []float64, vecs *core.Array, n, k int, s float64) float64 {
	cv := make([]complex128, n*k)
	for i := range n {
		for j := range k {
			cv[i*k+j] = complex(vecs.FloatAt(i*k+j), 0)
		}
	}
	cvals := make([]complex128, k)
	for j := range k {
		cvals[j] = complex(vals[j], 0)
	}
	return ritzResidual(a, cvals, cv, n, k, s)
}

// TestSparseSymmetryGuardIsRelativeAtSmallScale pins the fourth absolute
// floor (report F16, sparseigen.go checkSymmetric, the guard behind
// SpEigen, SpSolve and SpExpApply): an asymmetry that is a percent of a
// 1e-14 matrix, yet below the old 1e-15 floor, must be refused with the
// same wording as its dense and complex siblings, and the exactly
// symmetric matrix of that scale must still be accepted.
func TestSparseSymmetryGuardIsRelativeAtSmallScale(t *testing.T) {
	const s = 1e-14
	// Both mirrored entries are stored; the difference is 9e-17, below
	// the removed 1e-15 floor and above the purely relative 1e-26.
	asym := realCOO(t, []float64{s, 1e-16, 1e-17, s}, 2)
	b := mustFloats(t, []float64{s, s}, 2)
	if _, _, err := SpEigen(asym, 2, core.NewGenerator(7)); err == nil {
		t.Fatal("SpEigen: asymmetric at a percent of scale 1e-14 was accepted, want a refusal")
	} else if !strings.Contains(err.Error(), "SpEigen: matrix is not symmetric within 1e-12 tolerance") {
		t.Fatalf("SpEigen refusal = %q, want the sibling wording", err)
	}
	if _, err := SpSolve(asym, b, 1e-12, 100); err == nil {
		t.Fatal("SpSolve: asymmetric at a percent of scale 1e-14 was accepted, want a refusal")
	} else if !strings.Contains(err.Error(), "SpSolve: matrix is not symmetric within 1e-12 tolerance") {
		t.Fatalf("SpSolve refusal = %q, want the sibling wording", err)
	}
	if _, err := SpExpApply(asym, b, 0); err == nil {
		t.Fatal("SpExpApply: asymmetric at a percent of scale 1e-14 was accepted, want a refusal")
	} else if !strings.Contains(err.Error(), "SpExpApply: matrix is not symmetric within 1e-12 tolerance") {
		t.Fatalf("SpExpApply refusal = %q, want the sibling wording", err)
	}

	// The exactly symmetric matrix of the same scale passes all three.
	sym := realCOO(t, []float64{s, 1e-16, 1e-16, s}, 2)
	if _, _, err := SpEigen(sym, 2, core.NewGenerator(7)); err != nil {
		t.Fatalf("SpEigen of a symmetric matrix at scale %g: %v", s, err)
	}
	if _, err := SpSolve(sym, b, 1e-12, 100); err != nil {
		t.Fatalf("SpSolve of a symmetric matrix at scale %g: %v", s, err)
	}
	if _, err := SpExpApply(sym, b, 0); err != nil {
		t.Fatalf("SpExpApply of a symmetric matrix at scale %g: %v", s, err)
	}
}

// TestSpEigenSpectrumAtExtremeScale pins the real Hermitian Lanczos on the
// 2 / 0.5 stencil against the closed-form spectrum 2+cos(k·π/(n+1)). The
// projected tridiagonal's squared accumulation overflows in the symmetric
// sweep above about 1.3e154, which deflates the whole block at once and
// returns the raw Rayleigh quotients as the spectrum.
func TestSpEigenSpectrumAtExtremeScale(t *testing.T) {
	const n = 8
	want := topStencilValues(n)
	for _, s := range scales {
		coo := realCOO(t, scale(stencil(n), s), n)
		vals, vecs, err := SpEigen(coo, 2, core.NewGenerator(7))
		if err != nil {
			t.Fatalf("SpEigen at scale %g: %v", s, err)
		}
		for i := range 2 {
			got := vals.FloatAt(i) / s
			if math.IsNaN(got) || math.IsInf(got, 0) {
				t.Fatalf("scale %g: Ritz value %d = %v, want %g", s, i, vals.FloatAt(i), want[i])
			}
			if math.Abs(got-want[i]) > 1e-12 {
				t.Fatalf("scale %g: Ritz value %d = %g, want %g (in units of s)", s, i, got, want[i])
			}
		}
		if res := ritzResidualReal(coo, vals.RawFloats(), vecs, n, 2, s); res > 1e-12 {
			t.Fatalf("scale %g: Ritz residual %g, want <= 1e-12", s, res)
		}
	}
}

// TestSpEigenGeneralSpectrumAtExtremeScale pins both general sparse
// eigensolvers on the same closed-form spectrum. The Arnoldi recurrence
// computes its norms as a raw sum of squares on the unscaled operator, so
// above about 1.3e154 the projected coupling became +Inf and below about
// 1.5e-162 it became zero, which the exhaustion test read as a collapsed
// Krylov block: both extremes returned wrong Ritz values silently.
func TestSpEigenGeneralSpectrumAtExtremeScale(t *testing.T) {
	const n = 8
	want := topStencilValues(n)
	// The Hermitian stencil for the complex entry: 2 on the diagonal,
	// conjugate mirrored imaginary couplings, same closed-form spectrum.
	cb := hermitianStencil(n)
	// A real symmetric matrix is a valid, though not typical, input to
	// the general entry, and its real spectrum makes the reference
	// unambiguous.
	cScales := []float64{1, 1e150, 1e200, 1e-150, 1e-200}
	for _, s := range cScales {
		coo := realCOO(t, scale(stencil(n), s), n)
		vals, vecs, err := SpEigenGeneral(coo, 2, core.NewGenerator(5))
		if err != nil {
			t.Fatalf("SpEigenGeneral at scale %g: %v", s, err)
		}
		cv := vals.RawComplexes()
		for i := range 2 {
			got := cmplx.Abs(cv[i]) / s
			if math.Abs(got-want[i]) > 1e-12 {
				t.Fatalf("SpEigenGeneral scale %g: |λ| %d = %g, want %g", s, i, got, want[i])
			}
		}
		if res := ritzResidual(coo, cv, vecs.RawComplexes(), n, 2, s); res > 1e-12 {
			t.Fatalf("SpEigenGeneral scale %g: Ritz residual %g, want <= 1e-12", s, res)
		}

		ccoo := cSRFromDense(t, scaleC(cb, s), n)
		cvals, cvecs, err := SpEigenGeneralComplex(ccoo, 2, core.NewGenerator(5))
		if err != nil {
			t.Fatalf("SpEigenGeneralComplex at scale %g: %v", s, err)
		}
		ccv := cvals.RawComplexes()
		for i := range 2 {
			got := cmplx.Abs(ccv[i]) / s
			if math.Abs(got-want[i]) > 1e-12 {
				t.Fatalf("SpEigenGeneralComplex scale %g: |λ| %d = %g, want %g", s, i, got, want[i])
			}
		}
		if res := ritzResidual(ccoo, ccv, cvecs.RawComplexes(), n, 2, s); res > 1e-12 {
			t.Fatalf("SpEigenGeneralComplex scale %g: Ritz residual %g, want <= 1e-12", s, res)
		}
	}
}