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

package linalg

// SVD, Eigen, Pinverse, MatrixRank and Cond. Operates on
// float64 in and out; complex matrices are rejected.
//
// SVD: two-sided Jacobi rotation (Kogbetliantz / Press et al. NR3 §11.4).
// Both U and V are computed simultaneously by applying Jacobi
// rotations column-by-column on A and row-by-row on AᵀA; the
// algorithm is quadratic in the number of sweeps but converges
// robustly on ill-conditioned and rank-deficient matrices and is
// straightforward to keep correct. The matrix sizes tensor calls
// (n ≪ 256) keep this competitive with the QR approach.
//
// Eigen: implicit-shift QR on a symmetric tridiagonal matrix
// (Householder reduction first). Quadratic in sweeps but standard
// and stable. Both eigenvalues and eigenvectors are returned.

import (
	"math"
	"slices"
	"sourcedock.dev/petrbalvin/tensor/internal/base"
	"sourcedock.dev/petrbalvin/tensor/internal/core"
	"sourcedock.dev/petrbalvin/tensor/internal/engine"
)

// svdMinWork is the measured parallel gate for one reflector
// application or matrix product slice, counted in element touches
// (items dispatched times the per-item inner length). Below it the
// work finishes before the worker dispatch pays off and runs on the
// calling goroutine. The bulge chase itself (symmetricQr,
// tridiagQrStep) stays serial: each rotation reads what the previous
// one wrote, and even the rounding-level entries far from the band
// feed a bulge a few rotations later, so no part of the tMat sweep may
// be reordered or narrowed. Only the eigenvector accumulation is
// order-safe: its mixes touch qMat alone, so they are recorded during
// the chase and replayed in batches over disjoint row ranges.
const svdMinWork = 16384

// svdMinItemsPerWorker is the floor on items per worker in the SVD
// dispatches.
const svdMinItemsPerWorker = 8

// qrFlushCap bounds the recorded rotation buffer between batched
// applications, keeping the replay scratch at a fixed size for any n.
const qrFlushCap = 65536

// qrRotation records one Givens rotation accumulated on the right of
// the eigenvector matrix: coordinates k and k+1 rotate by
// G = [[c, −s], [s, c]].
type qrRotation struct {
	k    int
	c, s float64
}

// qrSweep collects the rotations of one symmetricQr run and applies
// them to the eigenvector accumulator in batches.
//
// Deferring the mixes is bit-identical to applying them inside the
// chase. A qMat cell (i, k) is touched only by rotations whose band
// position is k or k-1, each computing c·old + s·old' and -s·old +
// c·old' from the cell's current value and the rotation's (c, s); the
// pair (c, s) is read off tMat alone. Recording the rotations and
// replaying them in the same chronological order therefore feeds every
// cell exactly the arithmetic sequence the inline loop produced, and a
// flush only cuts that sequence into consecutive chunks. The replay
// walks disjoint row ranges per worker, so the parallel batch and the
// serial fallback write the same bits.
type qrSweep struct {
	qMat []float64
	n    int
	rots []qrRotation
	// inline mirrors the pre-deferral behaviour: below qrDeferMin the
	// chase applies each rotation to the accumulator as it runs,
	// because the buffer, the flush and the worker dispatch cost more
	// than the n-row mix they would save.
	inline bool
}

// qrDeferMin is the accumulator size from which the eigenvector mixes
// defer to batched replay. Measured on the square sweep: at 128 the
// deferral loses a third of the run to the buffer and the dispatch,
// at 192 it measures even, and at 256 the batched replay wins 1.6
// times, with the gap widening as the mix grows cubically.
const qrDeferMin = 192

func newQrSweep(qMat []float64, n int) *qrSweep {
	return &qrSweep{qMat: qMat, n: n, rots: make([]qrRotation, 0, 256),
		inline: n < qrDeferMin}
}

// record applies one rotation inline below qrDeferMin and otherwise
// appends it, flushing the buffer at the cap so the scratch stays
// bounded for large n. Both modes feed the cell the same arithmetic in
// the same chronological position, so the result does not depend on
// the mode.
func (qa *qrSweep) record(k int, c, s float64) {
	if qa.inline {
		qMat, n := qa.qMat, qa.n
		for i := range n {
			qi := qMat[i*n+k]
			qi1 := qMat[i*n+k+1]
			qMat[i*n+k] = c*qi + s*qi1
			qMat[i*n+k+1] = -s*qi + c*qi1
		}
		return
	}
	qa.rots = append(qa.rots, qrRotation{k: k, c: c, s: s})
	if len(qa.rots) >= qrFlushCap {
		qa.flush()
	}
}

// flush applies the recorded rotations to qMat, one worker per row
// range, and empties the buffer. Below the dispatch gate the batch
// runs on the calling goroutine.
func (qa *qrSweep) flush() {
	rots := qa.rots
	if len(rots) == 0 {
		return
	}
	n := qa.n
	qMat := qa.qMat
	apply := func(start, end int) {
		for i := start; i < end; i++ {
			row := i * n
			for _, r := range rots {
				qi := qMat[row+r.k]
				qi1 := qMat[row+r.k+1]
				qMat[row+r.k] = r.c*qi + r.s*qi1
				qMat[row+r.k+1] = -r.s*qi + r.c*qi1
			}
		}
	}
	if len(rots)*n >= svdMinWork {
		engine.ParallelMin(n, 1, func(start, end int) {
			apply(start, end)
		})
	} else {
		apply(0, n)
	}
	qa.rots = rots[:0]
}

// scaleWinLo and scaleWinHi bound the magnitude window the squared
// arithmetic of the decompositions is exact and finite in. Every square,
// sum of squares, product and fourth power of an entry inside it stays a
// normal float64: 2^200 squared is 2^400 and its fourth power 2^800, far
// inside the 1.8e308 ceiling, and 2^-300 squared is 2^-600, far above
// the 2^-1022 normal floor. Outside the window the raw accumulations
// underflow to zero, which silently skips reflectors, or overflow to
// +Inf, which zeroes a reflector's beta and poisons the update with NaN.
const (
	scaleWinLo = -300 // 2^-300 ≈ 4.9e-91
	scaleWinHi = 200  // 2^200  ≈ 1.6e60
)

// windowScale returns the exact power of two that moves maxAbs into
// [2^scaleWinLo, 2^scaleWinHi], and 1 when it already lies inside. A
// decomposition multiplies its matrix by the factor and its spectrum by
// 1/f: that is the same problem in scaled units, because eigenvalues and
// singular values scale with the matrix while eigenvectors and the
// orthogonal factors do not. Powers of two carry no rounding, so an
// in-window input keeps every bit of the original arithmetic and a
// rescaled one keeps the accuracy its exponent leaves it.
func windowScale(maxAbs float64) float64 {
	if maxAbs <= 0 || math.IsNaN(maxAbs) || math.IsInf(maxAbs, 0) {
		return 1
	}
	// Frexp splits maxAbs into m·2^e with m in [0.5, 1).
	_, e := math.Frexp(maxAbs)
	switch {
	case e-1 > scaleWinHi:
		return math.Ldexp(1, scaleWinHi-e)
	case e-1 < scaleWinLo:
		return math.Ldexp(1, scaleWinLo+1-e)
	}
	return 1
}

// maxMagF64 returns the largest absolute value in a real slice.
func maxMagF64(s []float64) float64 {
	m := 0.0
	for _, v := range s {
		if a := math.Abs(v); a > m {
			m = a
		}
	}
	return m
}

// maxMagComplex returns the largest magnitude in a complex slice. The
// magnitudes come from math.Hypot, so an entry beyond ~1e154 keeps a
// finite magnitude where a sum of raw squares overflows.
func maxMagComplex(s []complex128) float64 {
	m := 0.0
	for _, z := range s {
		if a := cmplxAbs(z); a > m {
			m = a
		}
	}
	return m
}

// scaleFloats multiplies every entry of s by f; f is an exact power of
// two, so the entries come back bit for bit from unscaleFloats.
func scaleFloats(s []float64, f float64) {
	for i := range s {
		s[i] *= f
	}
}

// unscaleFloats divides every entry of s by f, an exact power of two.
func unscaleFloats(s []float64, f float64) {
	for i := range s {
		s[i] /= f
	}
}

// scaleComplexes multiplies every entry of s by the real power of two f.
func scaleComplexes(s []complex128, f float64) {
	for i := range s {
		s[i] = complex(real(s[i])*f, imag(s[i])*f)
	}
}

// unscaleComplexes divides every entry of s by the real power of two f.
func unscaleComplexes(s []complex128, f float64) {
	for i := range s {
		s[i] = complex(real(s[i])/f, imag(s[i])/f)
	}
}

// SVD returns the thin singular value decomposition a = U · Σ · Vᵀ
// of an m×n matrix a with m ≥ n. For m < n the matrix is
// transposed first and the orthogonal factors are swapped on
// return. U is m×n with orthonormal columns (Uᵀ U = I_n); Σ is a
// 1-D length-n array of non-negative singular values sorted
// descending; Vᵀ is n×n orthogonal.
//
// Algorithm:
//  1. Bidiagonalise A, giving A = U₁ · B · V₁ᵀ via Householder.
//  2. Form C = Bᵀ B (n × n symmetric tridiagonal).
//  3. Symmetric implicit-shift QR with Wilkinson shift on C,
//     accumulating rotations into V = V₁ · Q.
//  4. Σ = √(eigenvalues of C). Reconstruct U = A · V / Σ and
//     orthogonalise via Gram-Schmidt to absorb numerical drift.
//
// This is the standard double-step SVD via Bᵀ B + symmetric QR:
// robust, well-understood, and reuses the symmetric tridiagonal
// eigensolver already shipped for `Eigen`.
//
// The double step prices the small end of the spectrum: the singular
// values come from the eigenvalues of BᵀB, so their accuracy on very
// small singular values is that of a squared condition number, fine
// for rank decisions, weaker for resolving near-null directions.
// `SVDComplex` states the same contract for its Aᴴ·A route.
func SVD(a *core.Array) (u, sigma, vt *core.Array, err error) {
	if a.Dtype() == core.Complex {
		return nil, nil, nil, base.Errf("SVD: complex matrices are not supported")
	}
	if a.NDim() != 2 {
		return nil, nil, nil, base.Errf("SVD: needs a 2-D matrix, got shape %s", base.ShapeText(a.Shape()))
	}
	m, n := a.Shape()[0], a.Shape()[1]
	if m == 0 || n == 0 {
		return nil, nil, nil, base.Errf("SVD: zero-sized matrix, got shape %s", base.ShapeText(a.Shape()))
	}
	transposed := false
	Amat := denseFloats(a, a.Shape()[0], a.Shape()[1])
	// Outside the safe window every squared intermediate of the
	// bidiagonalisation, of C = BᵀB and of the σ recovery would leave
	// the normal range. The matrix is moved into the window for the
	// computation and the singular values, which carry the scale, are
	// moved back once the factorisation is done; the orthogonal factors
	// are scale-free.
	ws := windowScale(maxMagF64(Amat))
	if ws != 1 {
		scaleFloats(Amat, ws)
	}
	origM, origN := m, n
	if m < n {
		transposed = true
		Amat = base.TransposeFlat(Amat, m, n)
		m, n = n, m
	}
	bDiag, bSuper, v1 := bidiagonalise(Amat, m, n)
	// Form C = Bᵀ B (n × n symmetric tridiagonal): AᵀA = V₁ C V₁ᵀ.
	// The diagonal carries the superdiagonal's contribution:
	// C[i,i] = B[i,i]² + B[i−1,i]²; the off-diagonal is
	// C[i,i+1] = B[i,i]·B[i,i+1].
	cMat := make([]float64, n*n)
	for i := 0; i < n; i++ {
		cMat[i*n+i] = bDiag[i] * bDiag[i]
		if i > 0 {
			cMat[i*n+i] += bSuper[i-1] * bSuper[i-1]
		}
		if i+1 < n {
			cMat[i*n+i+1] = bDiag[i] * bSuper[i]
			cMat[(i+1)*n+i] = cMat[i*n+i+1]
		}
	}
	// Eigendecompose C with the shipped symmetric solver: C = W Λ Wᵀ,
	// so the right singular vectors are V = V₁·W with Σ = √Λ. C is
	// (n × n) symmetric and tiny (n ≪ 1024 in the typical workspace),
	// so the eigensolver is much cheaper than a QR sweep per call.
	// C is exactly symmetric by construction, so the internals of
	// `Eigen` run directly.
	tMat, qMat := householderTridiag(cMat, n)
	qT := base.TransposeFlat(qMat, n, n)
	if err := symmetricQr(tMat, qT, n); err != nil {
		return nil, nil, nil, base.Errf("SVD: %w", err)
	}
	eVals := make([]float64, n)
	for i := range n {
		eVals[i] = tMat[i*n+i]
	}
	eIdx := sortAscIndices(eVals)
	cVals := make([]float64, n) // ascending eigenvalues, as Eigen returns them
	eVecs := make([]float64, n*n)
	for j := range n {
		cVals[j] = eVals[eIdx[j]]
		for i := range n {
			eVecs[i*n+j] = qT[i*n+eIdx[j]]
		}
	}
	vMat := make([]float64, n*n) // V = V₁·W
	// Rows are independent; each element keeps its serial dot over k.
	if n*n*n >= svdMinWork {
		engine.ParallelMin(n, 1, func(start, end int) {
			for i := start; i < end; i++ {
				for j := range n {
					s := 0.0
					for k := range n {
						s += v1[i*n+k] * eVecs[k*n+j]
					}
					vMat[i*n+j] = s
				}
			}
		})
	} else {
		for i := range n {
			for j := range n {
				s := 0.0
				for k := range n {
					s += v1[i*n+k] * eVecs[k*n+j]
				}
				vMat[i*n+j] = s
			}
		}
	}
	// Eigen returns ascending order; the SVD contract wants Σ
	// descending, so permute the columns of V to match.
	idx := sortDescIndices(cVals)
	permuteCols(vMat, n, n, idx)
	// Recover Σ directly from A as σₖ = ‖A·vₖ‖: the BᵀB round trip
	// squares the condition number, so reading √λ back would lose half
	// the digits on ill-conditioned inputs. Columns are independent and
	// each norm keeps its serial accumulation order over i, then j.
	sorted := make([]float64, n)
	if m*n*n >= svdMinWork {
		engine.ParallelMin(n, 1, func(start, end int) {
			for k := start; k < end; k++ {
				norm := 0.0
				for i := range m {
					acc := 0.0
					for j := range n {
						acc += Amat[i*n+j] * vMat[j*n+k]
					}
					norm += acc * acc
				}
				sorted[k] = math.Sqrt(norm)
			}
		})
	} else {
		for k := range n {
			norm := 0.0
			for i := range m {
				acc := 0.0
				for j := range n {
					acc += Amat[i*n+j] * vMat[j*n+k]
				}
				norm += acc * acc
			}
			sorted[k] = math.Sqrt(norm)
		}
	}
	// The recovered Σ can disagree with the eigenvalue order by
	// rounding on near-degenerate values; settle on one descending
	// order before building U, permuting V along with it.
	if ord := sortDescIndices(sorted); !isIdentityPerm(ord) {
		permuteCols(vMat, n, n, ord)
		perm := make([]float64, n)
		copy(perm, sorted)
		for i := range n {
			sorted[i] = perm[ord[i]]
		}
	}
	// Thin U = A·V·Σ⁻¹ from the working (m ≥ n) matrix, with A and Σ in
	// the same scaled units.
	uThin := reconstructUThin(Amat, vMat, sorted, m, n)
	// The singular values carry the matrix's scale and the orthogonal
	// factors do not, so only Σ is moved back out of the scaled units
	// the factorisation ran in.
	if ws != 1 {
		unscaleFloats(sorted, ws)
	}
	U := floatsToArray(uThin, []int{m, n})
	S := floatsToArray(sorted, []int{n})
	Vt := floatsToArray(base.TransposeFlat(vMat, n, n), []int{n, n})
	if transposed {
		// The input had m < n; the factorisation ran on the transpose
		// Aᵀ = U' Σ V'ᵀ, so the original reads A = V' Σ U'ᵀ:
		// U_orig = V' itself (origM × origM) and Vᵀ_orig = (thin U')ᵀ.
		return floatsToArray(vMat, []int{origM, origM}),
			S,
			floatsToArray(base.TransposeFlat(uThin, m, n), []int{origM, origN}),
			nil
	}
	return U, S, Vt, nil
}

// bidiagonalise reduces A (m×n with m ≥ n, row-major) to upper
// bidiagonal form via Golub-Kahan Householder reflectors: B =
// H_p…H₁·A·K₁…K_q, so A = U₁·B·V₁ᵀ with V₁ = K₁K₂…K_q accumulated by
// right-multiplying the reflectors into its columns. The left factor U₁
// is not accumulated: the only caller recovers U from A·V after the
// sweep, and carrying the m×m accumulator cost an eye(m) matrix and an
// m-row update per left reflector without feeding any returned value.
// The bidiagonal and the right accumulator are exactly the values the
// accumulation carried: neither reads uMat.
//
// Reflectors apply strictly in sweep order; within one reflector every
// touched column or row is updated exactly once from the frozen
// reflector vector and beta, so the per-line work (dot along the line,
// then the same-operand rank-1 subtraction) dispatches over the lines
// of bFull and the accumulator in one parallel range. The per-element
// sequence is the serial one, so the factorisation is bit-identical
// across crew sizes.
func bidiagonalise(aMat []float64, m, n int) (bDiag, bSuper []float64, vMat []float64) {
	bDiag = make([]float64, n)
	if n > 1 {
		bSuper = make([]float64, n-1)
	}
	bFull := make([]float64, m*n)
	copy(bFull, aMat)
	vMat = eye(n)
	// Reflector scratch reused across the sweep; each step uses the
	// prefix it needs.
	vec := make([]float64, m)
	vecR := make([]float64, n)
	for k := range n {
		// Left reflector: rows k..m-1 of column k, zeroing below (k, k).
		lv := vec[:m-k]
		for i := range lv {
			lv[i] = bFull[(k+i)*n+k]
		}
		hh := householderVectorInto(lv, lv)
		if hh.beta != 0 {
			ln := m - k
			bColumn := func(j int) {
				dot := 0.0
				for i := range ln {
					dot += hh.v[i] * bFull[(k+i)*n+j]
				}
				w := hh.beta * dot
				for i := range ln {
					bFull[(k+i)*n+j] -= hh.v[i] * w
				}
			}
			if (n-k)*ln >= svdMinWork {
				engine.ParallelMin(n-k, svdMinItemsPerWorker, func(start, end int) {
					for j := start; j < end; j++ {
						bColumn(k + j)
					}
				})
			} else {
				for j := k; j < n; j++ {
					bColumn(j)
				}
			}
		}
		if k+1 >= n {
			break
		}
		// Right reflector: columns k+1..n-1 of row k, zeroing right of
		// (k, k+1).
		rv := vecR[:n-k-1]
		for j := range rv {
			rv[j] = bFull[k*n+(k+1+j)]
		}
		hhR := householderVectorInto(rv, rv)
		if hhR.beta != 0 {
			ln := len(rv)
			bRow := func(i int) {
				dot := 0.0
				for j := range ln {
					dot += hhR.v[j] * bFull[i*n+(k+1+j)]
				}
				w := hhR.beta * dot
				for j := range ln {
					bFull[i*n+(k+1+j)] -= hhR.v[j] * w
				}
			}
			vRow := func(j int) {
				dot := 0.0
				for i := range ln {
					dot += hhR.v[i] * vMat[j*n+(k+1+i)]
				}
				w := hhR.beta * dot
				for i := range ln {
					vMat[j*n+(k+1+i)] -= hhR.v[i] * w
				}
			}
			// V₁ = K₁K₂…: likewise from the right, columns k+1..n-1.
			// Rows k..m-1 of bFull and rows of vMat: disjoint lines.
			items := (m - k) + n
			if items*ln >= svdMinWork {
				engine.ParallelMin(items, svdMinItemsPerWorker, func(start, end int) {
					for it := start; it < end; it++ {
						if it < m-k {
							bRow(k + it)
						} else {
							vRow(it - (m - k))
						}
					}
				})
			} else {
				for i := k; i < m; i++ {
					bRow(i)
				}
				for j := range n {
					vRow(j)
				}
			}
		}
	}
	for i := range n {
		bDiag[i] = bFull[i*n+i]
		if i+1 < n {
			bSuper[i] = bFull[i*n+i+1]
		}
	}
	return bDiag, bSuper, vMat
}

// smallOff is the deflation threshold for the tridiagonal QR sweep.
func smallOff(d1, d2 float64) float64 {
	return base.EpsF * (math.Abs(d1) + math.Abs(d2))
}

// reconstructUThin builds thin U = A_orig · V · diag(1/Σ) when Σ is
// non-zero, and falls back to the kernel projection for small
// σᵢ. aOrig is (post-transpose) input shape (m, n). vMat is (n, n).
// sVals are the singular values in matched column order. Result
// has shape (m, n).
//
// The candidate projection A_orig·V is computed once into a scratch,
// dispatched over the rows of aOrig: every element is the same serial
// dot over j (ascending) the per-column loop used to run, so the
// values are bit-identical and only their computing order moved. The
// modified Gram-Schmidt pass stays serial: column k reads every
// earlier column, and each projection's dot is an order-bound
// reduction over the column being updated.
func reconstructUThin(aOrig, vMat, sVals []float64, m, n int) []float64 {
	biggest := 0.0
	for _, v := range sVals {
		if v > biggest {
			biggest = v
		}
	}
	thresh := base.EpsF * float64(m) * biggest
	// The Gram-Schmidt pass tests the norm of a column that the division
	// by σ makes unit-length, so its floor is an absolute one relative to
	// that unit and not to the largest singular value: with the singular
	// values themselves as the reference, a matrix whose spectrum is
	// above ~1e16 would route every healthy column to the kernel rebuild
	// below and lose the factorisation.
	unitTol := base.EpsF * float64(m)
	// The un-orthonormalised candidate is A_orig·V/Σ on the columns
	// where σ is non-tiny (and A_orig·V otherwise), re-orthonormalised
	// via modified Gram-Schmidt below so that UᵀU = I_n holds even when
	// the reduction leaves small numerical drift. A column whose
	// residual is essentially zero (σ at rounding level) is rebuilt from
	// the best orthogonalised coordinate candidate e_j: the largest
	// residual is at least 1/√m, because a unit vector of the
	// orthogonal complement has a coordinate at least that large.
	uMat := make([]float64, m*n)
	// Scratch reused across columns: the working column and, for kernel
	// columns, the candidate buffers (fully overwritten before use).
	best := make([]float64, m)
	cand := make([]float64, m)
	// proj holds the un-orthonormalised projection A_orig·V, one value
	// per (row, singular vector) pair, each from the same serial dot.
	proj := engine.GetFloat64Buf(m * n)
	defer engine.PutFloat64Buf(proj)
	if m*n*n >= svdMinWork {
		engine.ParallelMin(m, 1, func(start, end int) {
			for i := start; i < end; i++ {
				row := aOrig[i*n : i*n+n]
				out := i * n
				for k := range n {
					s := 0.0
					for j := range n {
						s += row[j] * vMat[j*n+k]
					}
					proj[out+k] = s
				}
			}
		})
	} else {
		for i := range m {
			for k := range n {
				s := 0.0
				for j := range n {
					s += aOrig[i*n+j] * vMat[j*n+k]
				}
				proj[i*n+k] = s
			}
		}
	}
	for k := range n {
		if sVals[k] > thresh {
			inv := 1 / sVals[k]
			for i := range m {
				uMat[i*n+k] = proj[i*n+k] * inv
			}
		} else {
			for i := range m {
				uMat[i*n+k] = proj[i*n+k]
			}
		}
		// Subtract projections on the earlier columns.
		for j := range k {
			dot := 0.0
			for i := range m {
				dot += uMat[i*n+j] * uMat[i*n+k]
			}
			for i := range m {
				uMat[i*n+k] -= dot * uMat[i*n+j]
			}
		}
		nrm := 0.0
		for i := range m {
			nrm += uMat[i*n+k] * uMat[i*n+k]
		}
		nrm = math.Sqrt(nrm)
		if nrm > unitTol {
			inv := 1 / nrm
			for i := range m {
				uMat[i*n+k] *= inv
			}
			continue
		}
		// Kernel column: orthogonalise every coordinate candidate and
		// keep the healthiest residual.
		bestRes := -1.0
		for j := range m {
			for i := range m {
				cand[i] = 0
			}
			cand[j] = 1
			for range 2 {
				for j2 := range k {
					d := 0.0
					for i := range m {
						d += uMat[i*n+j2] * cand[i]
					}
					for i := range m {
						cand[i] -= d * uMat[i*n+j2]
					}
				}
			}
			rn := 0.0
			for i := range m {
				rn += cand[i] * cand[i]
			}
			if rn > bestRes {
				bestRes = rn
				copy(best, cand)
			}
		}
		if bestRes > 0 {
			inv := 1 / math.Sqrt(bestRes)
			for i := range m {
				uMat[i*n+k] = best[i] * inv
			}
		}
	}
	return uMat
}

// Eigen returns the eigenvalues and orthonormal eigenvectors of a
// real symmetric n×n matrix a. Eigenvalues are returned in a 1-D
// float array, sorted in ascending order; the eigenvectors are
// the columns of an n×n orthogonal array. Asymmetric matrices are
// not supported.
func Eigen(a *core.Array) (values, vectors *core.Array, err error) {
	if a.Dtype() == core.Complex {
		return nil, nil, base.Errf("Eigen: complex matrices are not supported")
	}
	if a.NDim() != 2 || a.Shape()[0] != a.Shape()[1] {
		return nil, nil, base.Errf("Eigen: needs a square 2-D matrix, got shape %s", base.ShapeText(a.Shape()))
	}
	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
	// window and the eigenvalues are scaled back. The eigenvectors are
	// the same for every positive multiple of the matrix.
	ws := windowScale(maxMagF64(mat))
	if ws != 1 {
		scaleFloats(mat, ws)
	}
	if !isSymmetric(mat, n) {
		return nil, nil, base.Errf("Eigen: matrix is not symmetric within 1e-12 tolerance")
	}
	tMat, qMat := householderTridiag(mat, n)
	// The reduction maintains T = qMat·A·qMatᵀ (reflectors applied on
	// the left of qMat), so the eigenvector matrix of A after the QR
	// sweep is qMatᵀ·G_total. The sweep accumulates rotations on the
	// right of its accumulator, hence the transpose before the call.
	qT := base.TransposeFlat(qMat, n, n)
	if err := symmetricQr(tMat, qT, n); err != nil {
		return nil, nil, base.Errf("Eigen: %w", err)
	}
	vals := make([]float64, n)
	for i := range n {
		vals[i] = tMat[i*n+i]
	}
	idx := sortAscIndices(vals)
	out := make([]float64, n)
	for i := range n {
		out[i] = vals[idx[i]]
	}
	sortedCols := make([]float64, n*n)
	for j := range n {
		for i := range n {
			sortedCols[i*n+j] = qT[i*n+idx[j]]
		}
	}
	if ws != 1 {
		unscaleFloats(out, ws)
	}
	return floatsToArray(out, []int{n}), floatsToArray(sortedCols, []int{n, n}), nil
}

// Pinverse returns the Moore-Penrose pseudoinverse of a 2-D matrix,
// built from the SVD by inverting singular values strictly greater
// than ε. When eps ≤ 0 the default is max(m, n) · max(Σ) · base.EpsF.
func Pinverse(a *core.Array, eps float64) (*core.Array, error) {
	if a.Dtype() == core.Complex {
		return nil, base.Errf("Pinverse: complex matrices are not supported")
	}
	if a.NDim() != 2 {
		return nil, base.Errf("Pinverse: needs a 2-D matrix, got shape %s", base.ShapeText(a.Shape()))
	}
	m, n := a.Shape()[0], a.Shape()[1]
	if m == 0 || n == 0 {
		return nil, base.Errf("Pinverse: zero-sized matrix, got shape %s", base.ShapeText(a.Shape()))
	}
	u, sigma, vt, err := SVD(a)
	if err != nil {
		return nil, err
	}
	sVals := sigma.RawFloats()
	// The thin factorisation carries min(m, n) singular values: U is
	// (m, r) and Vᵀ is (r, n), whatever the input aspect ratio.
	r := min(m, n)
	maxS := 0.0
	for _, s := range sVals {
		if s > maxS {
			maxS = s
		}
	}
	if eps <= 0 {
		dim := max(n, m)
		eps = float64(dim) * maxS * base.EpsF
	}
	dInv := make([]float64, r)
	for i, s := range sVals {
		if s > eps {
			dInv[i] = 1 / s
		}
	}
	uMat := denseFloats(u, m, r)
	vtMat := denseFloats(vt, r, n)
	out := make([]float64, n*m)
	// A⁺ = V · Σ⁻¹ · Uᵀ = (V · Σ⁻¹) · Uᵀ, shape (n, m):
	// out[i,j] = Σ_k V[i,k] · Σ⁻¹[k] · U[j,k], with V[i,k] = Vᵀ[k,i].
	for i := range n {
		for j := range m {
			s := 0.0
			for k := range r {
				s += vtMat[k*n+i] * dInv[k] * uMat[j*r+k]
			}
			out[i*m+j] = s
		}
	}
	return floatsToArray(out, []int{n, m}), nil
}

// MatrixRank returns the number of singular values of a 2-D matrix
// strictly greater than eps. With eps ≤ 0 the default is
// max(m, n) · max(Σ) · base.EpsF.
func MatrixRank(a *core.Array, eps float64) (int, error) {
	if a.Dtype() == core.Complex {
		return 0, base.Errf("MatrixRank: complex matrices are not supported")
	}
	if a.NDim() != 2 {
		return 0, base.Errf("MatrixRank: needs a 2-D matrix, got shape %s", base.ShapeText(a.Shape()))
	}
	_, sigma, _, err := SVD(a)
	if err != nil {
		return 0, err
	}
	return countAboveThreshold(sigma.RawFloats(), eps, a.Shape()[0], a.Shape()[1]), nil
}

// Cond returns the 2-norm condition number σ_max / σ_min of a 2-D
// matrix. If any singular value is zero (or at-or-below a positive
// eps) the result is +Inf, mirroring the standard convention.
func Cond(a *core.Array, eps float64) (float64, error) {
	if a.Dtype() == core.Complex {
		return 0, base.Errf("Cond: complex matrices are not supported")
	}
	if a.NDim() != 2 {
		return 0, base.Errf("Cond: needs a 2-D matrix, got shape %s", base.ShapeText(a.Shape()))
	}
	_, sigma, _, err := SVD(a)
	if err != nil {
		return 0, err
	}
	return condValue(sigma.RawFloats(), eps, a.Shape()[0], a.Shape()[1]), nil
}

// householderTridiag reduces a symmetric matrix to tridiagonal form via
// Householder reflectors, accumulating the orthogonal factor into qMat.
//
// The sweep order is untouched. Within one reflector the LEFT
// application to T (columns) and the accumulation into qMat (columns)
// run over disjoint buffers, so they share one dispatch; the RIGHT
// application to T (rows) must follow, because it reads the entries the
// left application has just written; the drift symmetrisation follows
// it serially. Each line is updated exactly once from the frozen
// reflector with the serial dot order, so the result is bit-identical.
func householderTridiag(aMat []float64, n int) (tMat, qMat []float64) {
	tMat = make([]float64, n*n)
	copy(tMat, aMat)
	qMat = eye(n)
	// Reflector scratch reused across the sweep.
	v := make([]float64, n)
	for k := 0; k < n-2; k++ {
		vv := v[:n-k-1]
		for i := 0; i < n-k-1; i++ {
			vv[i] = tMat[(k+1+i)*n+k]
		}
		u := householderVectorInto(vv, vv)
		if u.beta == 0 {
			continue
		}
		ln := n - k - 1
		tColumn := func(j int) {
			dot := 0.0
			for i := range ln {
				dot += u.v[i] * tMat[(k+1+i)*n+j]
			}
			w := u.beta * dot
			for i := range ln {
				tMat[(k+1+i)*n+j] -= u.v[i] * w
			}
		}
		qColumn := func(j int) {
			dot := 0.0
			for i := range ln {
				dot += u.v[i] * qMat[(k+1+i)*n+j]
			}
			w := u.beta * dot
			for i := range ln {
				qMat[(k+1+i)*n+j] -= u.v[i] * w
			}
		}
		// Apply H from both sides to T: left pass over the columns of T
		// together with the qMat accumulation (disjoint buffers).
		items := (n - k) + n
		if items*ln >= svdMinWork {
			engine.ParallelMin(items, svdMinItemsPerWorker, func(start, end int) {
				for it := start; it < end; it++ {
					if it < n-k {
						tColumn(k + it)
					} else {
						qColumn(it - (n - k))
					}
				}
			})
		} else {
			for j := k; j < n; j++ {
				tColumn(j)
			}
			for j := range n {
				qColumn(j)
			}
		}
		// Right pass over the rows of T: reads the left pass's writes,
		// so it stays strictly after it.
		tRow := func(i int) {
			dot := 0.0
			for j := range ln {
				dot += u.v[j] * tMat[i*n+(k+1)+j]
			}
			w := u.beta * dot
			for j := range ln {
				tMat[i*n+(k+1)+j] -= u.v[j] * w
			}
		}
		if (n-k)*ln >= svdMinWork {
			engine.ParallelMin(n-k, svdMinItemsPerWorker, func(start, end int) {
				for i := start; i < end; i++ {
					tRow(k + i)
				}
			})
		} else {
			for i := k; i < n; i++ {
				tRow(i)
			}
		}
		// Symmetrise against numerical drift.
		for i := k + 1; i < n; i++ {
			for j := i + 1; j < n; j++ {
				avg := (tMat[i*n+j] + tMat[j*n+i]) / 2
				tMat[i*n+j] = avg
				tMat[j*n+i] = avg
			}
		}
	}
	return tMat, qMat
}

// symmetricQr runs the implicit-shift symmetric QR algorithm with
// the Wilkinson shift on a symmetric tridiagonal tMat. Returns the
// updated tridiagonal (whose diagonal is the eigenvalues) and the
// eigenvector matrix Q, or an error when the sweep exhausts its
// passes without deflating every subdiagonal entry: an unconverged
// spectrum is never returned as an answer.
func symmetricQr(tMat, qMat []float64, n int) error {
	maxIters := max(30*n, 30)
	// Deflation floor relative to the matrix norm: a purely
	// neighbour-relative threshold never triggers when both diagonal
	// entries are small, even though the off-diagonal is already at
	// rounding level of the whole transform (backward stability
	// guarantees nothing better than eps·‖T‖ anyway).
	scale := 0.0
	for i := range n {
		scale += tMat[i*n+i] * tMat[i*n+i]
		if i+1 < n {
			e := tMat[i*n+i+1]
			scale += 2 * e * e
		}
	}
	tolAbs := base.EpsF * math.Sqrt(scale)
	negligible := func(i int) bool { // tests the subdiagonal entry (i, i-1)
		e := math.Abs(tMat[i*n+(i-1)])
		// An exactly zero entry is deflated by definition: a zero
		// matrix has tolAbs 0 and relative floors of 0, and a strict
		// inequality would never deflate it.
		return e == 0 || e < smallOff(tMat[i*n+i], tMat[(i-1)*n+(i-1)]) || e < tolAbs
	}
	acc := newQrSweep(qMat, n)
	for range maxIters {
		h := n - 1
		for h > 0 && negligible(h) {
			tMat[h*n+(h-1)] = 0
			tMat[(h-1)*n+h] = 0
			h--
		}
		if h <= 0 {
			break
		}
		// The active block containing h reaches up to the first
		// negligible subdiagonal entry (its decoupling boundary):
		// expanding l past that boundary would leave a zero entry
		// inside the window, where the bulge chase would abort and
		// strand the block below it.
		l := h - 1
		for l > 0 && !negligible(l) {
			l--
		}
		if l > 0 {
			tMat[l*n+(l-1)] = 0
			tMat[(l-1)*n+l] = 0
		}
		if h <= l {
			continue
		}
		if h-l == 1 {
			// A 2×2 block is diagonalised exactly by one Jacobi
			// rotation; the shift chase leaves a rounding-level
			// off-diagonal that can sit forever above the deflation
			// threshold, so finish it off directly.
			jacobi2x2(tMat, acc, l, n)
			tMat[l*n+h] = 0
			tMat[h*n+l] = 0
			continue
		}
		d := tMat[(h-1)*n+(h-1)]
		e := tMat[(h-1)*n+h]
		f := tMat[h*n+h]
		shift := wilkinsonShift2x2(d, e, f)
		tridiagQrStep(tMat, acc, l, h, n, shift)
	}
	// An exhausted sweep must not pass its tridiagonal off as a
	// spectrum: every sibling iterator here (hessenbergQr, the
	// Golub-Reinsch SVD, schurQR) errors on exhaustion, and so does
	// this one.
	for i := 1; i < n; i++ {
		if !negligible(i) {
			return base.Errf("the QR sweep failed to deflate a %d×%d matrix in %d passes", n, n, maxIters)
		}
	}
	acc.flush()
	return nil
}

// applyGivens applies the rotation G = [[c, −s], [s, c]] on coordinates
// (k, k+1) as a similarity transform of tMat and records it for the
// batched accumulation on the right of the eigenvector matrix.
func applyGivens(tMat []float64, acc *qrSweep, k, n int, c, s float64) {
	for j := range n {
		tk := tMat[k*n+j]
		tk1 := tMat[(k+1)*n+j]
		tMat[k*n+j] = c*tk + s*tk1
		tMat[(k+1)*n+j] = -s*tk + c*tk1
	}
	for i := range n {
		ti := tMat[i*n+k]
		ti1 := tMat[i*n+k+1]
		tMat[i*n+k] = c*ti + s*ti1
		tMat[i*n+k+1] = -s*ti + c*ti1
	}
	acc.record(k, c, s)
}

// jacobi2x2 diagonalises the symmetric 2×2 block at (k, k+1) exactly
// with one Jacobi rotation, recording the rotation into acc.
func jacobi2x2(tMat []float64, acc *qrSweep, k, n int) {
	a := tMat[k*n+k]
	b := tMat[k*n+k+1]
	d := tMat[(k+1)*n+(k+1)]
	if b == 0 {
		return
	}
	tau := (d - a) / (2 * b)
	// Smaller root of t² + ((a−d)/b)·t − 1 = 0: t = −sign(τ)/(|τ|+√(1+τ²)).
	t := -1.0 / (math.Abs(tau) + math.Sqrt(1+tau*tau))
	if tau < 0 {
		t = -t
	}
	c := 1 / math.Sqrt(1+t*t)
	s := t * c
	applyGivens(tMat, acc, k, n, c, s)
}
func wilkinsonShift2x2(d, e, f float64) float64 {
	delta := (d - f) / 2
	if delta == 0 {
		return f - math.Abs(e)
	}
	sq := math.Sqrt(delta*delta + e*e)
	if delta > 0 {
		return f - e*e/(delta+sq)
	}
	return f + e*e/(sq-delta)
}

// tridiagQrStep performs one implicit-shift QR step on the symmetric
// tridiagonal block [l, h] of tMat, recording the rotations into acc
// for the batched eigenvector accumulation. Every rotation is a
// similarity transform applied to full rows and columns, so tMat stays
// exactly similar to its input; the pair (c, s) comes from the
// subdiagonal and the bulge the previous rotation left behind, the
// bulge chase of Golub & Van Loan §8.5.
//
// The tMat passes are strictly sequential and must stay so: rotation
// k+1 annihilates the bulge rotation k wrote, and the rounding-level
// entries far from the band drift inward one diagonal per pass, so
// narrowing or reordering any pass would move bits. Only the qMat mix
// is free of that chain; it is replayed from the recording, in order,
// by acc.flush.
func tridiagQrStep(tMat []float64, acc *qrSweep, l, h, n int, shift float64) {
	// The first rotation is read from the leading entry of (T − μI);
	// later ones annihilate the chased bulge.
	x := tMat[l*n+l] - shift
	y := tMat[l*n+(l+1)]
	for k := l; k < h; k++ {
		if k > l {
			x = tMat[(k-1)*n+k]     // subdiagonal entry
			y = tMat[(k-1)*n+(k+1)] // bulge to annihilate
		}
		r := math.Hypot(x, y)
		if r == 0 {
			return // nothing to rotate; the band is already clean
		}
		c, s := x/r, y/r
		// Rows k, k+1: left multiplication by Gᵀ, G = [[c, −s], [s, c]].
		for j := range n {
			tk := tMat[k*n+j]
			tk1 := tMat[(k+1)*n+j]
			tMat[k*n+j] = c*tk + s*tk1
			tMat[(k+1)*n+j] = -s*tk + c*tk1
		}
		// Columns k, k+1: right multiplication by G.
		for i := range n {
			ti := tMat[i*n+k]
			ti1 := tMat[i*n+k+1]
			tMat[i*n+k] = c*ti + s*ti1
			tMat[i*n+k+1] = -s*ti + c*ti1
		}
		// Eigenvector accumulation on the right: Q = Q·G, replayed in
		// order by acc.flush.
		acc.record(k, c, s)
	}
}

// householder represents a Householder reflection v and its beta.
type householder struct {
	v    []float64
	beta float64
}

// householderVectorInto builds the standard reflect-to-first-coordinate
// Householder vector for x, where H = I - beta v vᵀ with H x = sign(x₀)‖x‖ e₁,
// writing the reflector into dst (which must have room for len(x)
// entries and may alias x) so sweep callers can reuse one buffer
// instead of allocating per reflection. Squared magnitudes are summed
// relative to the largest entry so a vector with entries near 1e154
// does not overflow on the way to its norm.
//
// beta = 2/(vᵀv) needs the square of v's largest entry to stay in the
// normal range. Past about 1.3e154 that square is +Inf and beta comes
// out as +0, which callers read as "no reflection at all", and below
// about 1e-154 it is subnormal or 0 and beta comes out as +Inf, whose
// product with a zero dot is NaN. Either way the reflector is silently
// lost. When the largest entry falls outside that range the reflector is
// re-expressed in units of the exact power of two that brings it into
// the safe window: v is scaled and beta carries the compensatory square,
// so H is the same reflection with every intermediate O(1).
func householderVectorInto(dst, x []float64) householder {
	maxAbs := 0.0
	for _, xi := range x {
		if v := math.Abs(xi); v > maxAbs {
			maxAbs = v
		}
	}
	if maxAbs == 0 {
		return householder{v: x, beta: 0}
	}
	scaled2 := 0.0
	for _, xi := range x {
		xi /= maxAbs
		scaled2 += xi * xi
	}
	norm := maxAbs * math.Sqrt(scaled2)
	dst[0] = x[0] + signOrNonZero(x[0])*norm
	copy(dst[1:], x[1:])
	vMax := 0.0
	for _, vi := range dst {
		if v := math.Abs(vi); v > vMax {
			vMax = v
		}
	}
	if vMax == 0 {
		return householder{v: dst, beta: 0}
	}
	v2 := 0.0
	for _, vi := range dst {
		vi /= vMax
		v2 += vi * vi
	}
	if beta := 2 / (vMax * vMax * v2); beta > 0 && beta < math.MaxFloat64 {
		return householder{v: dst, beta: beta}
	}
	ws := windowScale(vMax)
	scaleFloats(dst, ws)
	// w keeps max|v| inside the safe window, so w*w is a normal finite
	// number and 2/(w*w*v2) is finite and non-zero: the same reflection,
	// still built so that H x = sign(x₀)‖x‖ e₁.
	w := vMax * ws
	return householder{v: dst, beta: 2 / (w * w * v2)}
}

// signOrNonZero returns the sign of x, with 0 mapping to +1 (Householder
// convention).
func signOrNonZero(x float64) float64 {
	if x >= 0 {
		return 1
	}
	return -1
}

// isSymmetric reports whether the matrix is symmetric within a purely
// relative tolerance: a floor would admit an asymmetry that is a large
// fraction of a small-scale matrix, so the guard would call a matrix
// symmetric that is not, whatever the absolute size of its entries.
func isSymmetric(mat []float64, n int) bool {
	scale := 0.0
	for i := range n {
		for j := range n {
			a := math.Abs(mat[i*n+j])
			if a > scale {
				scale = a
			}
		}
	}
	tol := 1e-12 * scale
	for i := range n {
		for j := i + 1; j < n; j++ {
			if math.Abs(mat[i*n+j]-mat[j*n+i]) > tol {
				return false
			}
		}
	}
	return true
}

func countAboveThreshold(s []float64, eps float64, m, n int) int {
	if eps <= 0 {
		dim := max(n, m)
		biggest := 0.0
		for _, v := range s {
			if v > biggest {
				biggest = v
			}
		}
		eps = float64(dim) * biggest * base.EpsF
	}
	count := 0
	for _, v := range s {
		if v > eps {
			count++
		}
	}
	return count
}

func condValue(s []float64, eps float64, m, n int) float64 {
	if len(s) == 0 {
		return 0
	}
	biggest := 0.0
	smallest := math.Inf(1)
	for _, v := range s {
		if v > biggest {
			biggest = v
		}
		if v < smallest {
			smallest = v
		}
	}
	if biggest == 0 {
		// The zero matrix has no usable inverse direction; the condition
		// number is infinite, matching the doc contract.
		return math.Inf(1)
	}
	if smallest <= 0 {
		return math.Inf(1)
	}
	if eps > 0 && smallest <= eps {
		return math.Inf(1)
	}
	return biggest / smallest
}

func permuteCols(a []float64, rows, cols int, idx []int) {
	out := make([]float64, rows*cols)
	for j := range cols {
		for i := range rows {
			out[i*cols+j] = a[i*cols+idx[j]]
		}
	}
	copy(a, out)
}

// sortDescIndices orders the indices so the values descend. The sort is
// stable, matching the insertion sort it replaced bit for bit: ties keep
// their original order, so the returned permutation is identical.
func sortDescIndices(s []float64) []int {
	idx := make([]int, len(s))
	for i := range idx {
		idx[i] = i
	}
	slices.SortStableFunc(idx, func(a, b int) int {
		switch {
		case s[a] > s[b]:
			return -1
		case s[a] < s[b]:
			return 1
		default:
			return 0
		}
	})
	return idx
}

// sortAscIndices orders the indices so the values ascend, stable like
// sortDescIndices.
func sortAscIndices(s []float64) []int {
	idx := make([]int, len(s))
	for i := range idx {
		idx[i] = i
	}
	slices.SortStableFunc(idx, func(a, b int) int {
		switch {
		case s[a] < s[b]:
			return -1
		case s[a] > s[b]:
			return 1
		default:
			return 0
		}
	})
	return idx
}

// isIdentityPerm reports whether the permutation is the identity.
func isIdentityPerm(p []int) bool {
	for i, v := range p {
		if v != i {
			return false
		}
	}
	return true
}