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

// Behaviour pins for the package's entry points: every test here
// holds a corrected defect to the answer it must keep, among them the
// NaN quantile hang, the rolling dtype panics, the Quantile NaN
// panic, the Histogram2D complex panic and unbounded bin counts, the
// lower-half NormalQuantile bracket, the saturated GLM likelihood,
// the weighted regression intercept statistics, the missing
// finiteness and symmetry refusals, the BinomialDraws count guard,
// the Kolmogorov series' truncation, the NaN-parameter and
// complex-input sweeps, the KolmogorovSmirnovTest NaN loop and the
// degenerate-response F statistic. Each one fails with the old
// behaviour restored.

package stats

import (
	"math"
	"math/big"
	"testing"
	"time"

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

// pinVector builds a rank-1 float array.
func pinVector(t *testing.T, vals ...float64) *core.Array {
	t.Helper()
	a, err := core.FromFloats(vals, len(vals))
	if err != nil {
		t.Fatalf("FromFloats(%v): %v", vals, err)
	}
	return a
}

// pinWatchdog runs f in its own goroutine and fails the test if it
// has not returned within limit. It is what keeps the hang regression a
// failure rather than a stuck test binary.
func pinWatchdog(t *testing.T, limit time.Duration, what string, f func()) {
	t.Helper()
	done := make(chan struct{})
	go func() {
		defer close(done)
		f()
	}()
	select {
	case <-done:
	case <-time.After(limit):
		t.Fatalf("%s has not returned after %v: the loop does not terminate", what, limit)
	}
}

// TestDiscreteQuantileNaNQuantileDoesNotHang pins the P0: q = NaN passed
// the written-out range guard (`q < 0 || q > 1` is false for NaN), every
// bracket comparison below was then false too, and the doubling loop
// never terminated. Both discrete quantiles must refuse NaN instead.
func TestDiscreteQuantileNaNQuantileDoesNotHang(t *testing.T) {
	for _, tc := range []struct {
		name string
		call func() (float64, error)
	}{
		{"PoissonQuantile", func() (float64, error) { return PoissonQuantile(math.NaN(), 3) }},
		{"BinomialQuantile", func() (float64, error) { return BinomialQuantile(math.NaN(), 0.5, 5) }},
	} {
		t.Run(tc.name, func(t *testing.T) {
			var err error
			pinWatchdog(t, 10*time.Second, tc.name+"(NaN)", func() {
				_, err = tc.call()
			})
			if err == nil {
				t.Fatalf("%s(NaN) returned a value, want a refusal", tc.name)
			}
		})
	}
}

// bigSqrt2 is sqrt(2) at the working precision.
func bigSqrt2() *big.Float {
	const prec = 256
	return new(big.Float).SetPrec(prec).Sqrt(new(big.Float).SetPrec(prec).SetInt64(2))
}

// pinIntArray builds a rank-1 int array.
func pinIntArray(t *testing.T, vals ...int64) *core.Array {
	t.Helper()
	a, err := core.FromInts(vals, len(vals))
	if err != nil {
		t.Fatalf("FromInts(%v): %v", vals, err)
	}
	return a
}

// pinFloat32Array builds a rank-1 float32 array.
func pinFloat32Array(t *testing.T, vals ...float32) *core.Array {
	t.Helper()
	a, err := core.FromFloat32s(vals, len(vals))
	if err != nil {
		t.Fatalf("FromFloat32s(%v): %v", vals, err)
	}
	return a
}

// pinSameFloats pins two float slices to the bit for the dtype tests.
func pinSameFloats(t *testing.T, label string, got, want []float64) {
	t.Helper()
	if len(got) != len(want) {
		t.Fatalf("%s has length %d, want %d", label, len(got), len(want))
	}
	for i := range got {
		if got[i] != want[i] {
			t.Fatalf("%s[%d] = %.17g, want %.17g", label, i, got[i], want[i])
		}
	}
}

// TestRollingFamilyAcceptsIntAndFloat32 pins the P0: rollingCheck
// returned RawFloats()[:n], which is a nil slice with capacity 0 for an
// int or a float32 payload, so all four reductions panicked. The int,
// float32 and float routes must agree value for value.
func TestRollingFamilyAcceptsIntAndFloat32(t *testing.T) {
	vals := []float64{1, 4, 2, 8, 5, 7}
	floats := pinVector(t, vals...)
	ints := pinIntArray(t, 1, 4, 2, 8, 5, 7)
	f32 := pinFloat32Array(t, 1, 4, 2, 8, 5, 7)
	cases := []struct {
		name string
		call func(a *core.Array) (*core.Array, error)
	}{
		{"RollingMean", func(a *core.Array) (*core.Array, error) { return RollingMean(a, 3) }},
		{"RollingSum", func(a *core.Array) (*core.Array, error) { return RollingSum(a, 2) }},
		{"RollingMax", func(a *core.Array) (*core.Array, error) { return RollingMax(a, 3) }},
		{"RollingMin", func(a *core.Array) (*core.Array, error) { return RollingMin(a, 2) }},
	}
	for _, tc := range cases {
		t.Run(tc.name, func(t *testing.T) {
			want, err := tc.call(floats)
			if err != nil {
				t.Fatalf("float route: %v", err)
			}
			wantVals := make([]float64, want.Len())
			for i := range wantVals {
				wantVals[i] = want.FloatAt(i)
			}
			for _, dt := range []struct {
				name string
				a    *core.Array
			}{{"int", ints}, {"float32", f32}} {
				got, err := tc.call(dt.a)
				if err != nil {
					t.Fatalf("%s route: %v", dt.name, err)
				}
				gotVals := make([]float64, got.Len())
				for i := range gotVals {
					gotVals[i] = got.FloatAt(i)
				}
				pinSameFloats(t, dt.name, gotVals, wantVals)
			}
		})
	}
}

// TestRollingExtremesAllNaNWindow pins the P2: the window fold seeded
// its running best with ±Inf and only replaced it on a strict
// comparison, so an all-NaN window answered -Inf/+Inf where core.Min
// and core.Max answer NaN. A NaN that shares a window with a number
// must still lose.
func TestRollingExtremesAllNaNWindow(t *testing.T) {
	allNaN := pinVector(t, math.NaN(), math.NaN())
	max, err := RollingMax(allNaN, 2)
	if err != nil {
		t.Fatalf("RollingMax: %v", err)
	}
	if !math.IsNaN(max.FloatAt(0)) {
		t.Fatalf("RollingMax of an all-NaN window = %v, want NaN as core.Max gives", max.FloatAt(0))
	}
	min, err := RollingMin(allNaN, 2)
	if err != nil {
		t.Fatalf("RollingMin: %v", err)
	}
	if !math.IsNaN(min.FloatAt(0)) {
		t.Fatalf("RollingMin of an all-NaN window = %v, want NaN as core.Min gives", min.FloatAt(0))
	}
	// A NaN never wins against a number.
	mixed := pinVector(t, math.NaN(), 5, 2, math.NaN())
	wantMax := []float64{5, 5, 2}
	wantMin := []float64{5, 2, 2}
	gotMax, err := RollingMax(mixed, 2)
	if err != nil {
		t.Fatalf("RollingMax: %v", err)
	}
	gotMin, err := RollingMin(mixed, 2)
	if err != nil {
		t.Fatalf("RollingMin: %v", err)
	}
	for i := range wantMax {
		if gotMax.FloatAt(i) != wantMax[i] {
			t.Fatalf("RollingMax[%d] = %v, want %v (NaN never wins)", i, gotMax.FloatAt(i), wantMax[i])
		}
		if gotMin.FloatAt(i) != wantMin[i] {
			t.Fatalf("RollingMin[%d] = %v, want %v (NaN never wins)", i, gotMin.FloatAt(i), wantMin[i])
		}
	}
}

// TestQuantileRejectsNaNQuantile pins the P0: `q < 0 || q > 1` is
// false for NaN, so the quantile reached int(math.Floor(NaN)) and
// indexed the sorted values with it, panicking instead of refusing.
func TestQuantileRejectsNaNQuantile(t *testing.T) {
	a := pinVector(t, 1, 2, 3, 4)
	for _, qs := range [][]float64{{math.NaN()}, {0.5, math.NaN()}, {math.Inf(1)}, {math.Inf(-1)}} {
		out, err := Quantile(a, qs)
		if err == nil {
			t.Fatalf("Quantile(%v) = %v, want a refusal", qs, out)
		}
	}
	// The documented domain still works.
	out, err := Quantile(a, []float64{0.25, 0.5, 0.75})
	if err != nil {
		t.Fatalf("Quantile: %v", err)
	}
	pinSameFloats(t, "Quantile", []float64{out.FloatAt(0), out.FloatAt(1), out.FloatAt(2)},
		[]float64{1.75, 2.5, 3.25})
}

// TestHistogram2DRejectsComplexInput pins the P0: every other entry
// point of the package refuses complex input by name, while Histogram2D
// walked its float accessor into a nil int payload and panicked.
func TestHistogram2DRejectsComplexInput(t *testing.T) {
	real := pinVector(t, 1, 2, 3)
	cx, err := core.FromComplexes([]complex128{1, 2, 3}, 3)
	if err != nil {
		t.Fatalf("FromComplexes: %v", err)
	}
	if _, _, _, err := Histogram2D(cx, real, 2, 2); err == nil {
		t.Fatal("complex x accepted")
	}
	if _, _, _, err := Histogram2D(real, cx, 2, 2); err == nil {
		t.Fatal("complex y accepted")
	}
	// The real path is untouched.
	counts, _, _, err := Histogram2D(real, real, 2, 2)
	if err != nil {
		t.Fatalf("Histogram2D: %v", err)
	}
	if counts.Shape()[0] != 2 || counts.Shape()[1] != 2 {
		t.Fatalf("shape %v, want [2 2]", counts.Shape())
	}
}

// TestHistogramBinCountsAreBounded pins the P1: neither Histogram nor
// Histogram2D bounded its bin count, so a hostile request either asked
// the allocator for an unbounded buffer or, with a product that wraps,
// produced an empty count payload the first write then indexed. The
// guard must fire before anything is allocated.
func TestHistogramBinCountsAreBounded(t *testing.T) {
	a := pinVector(t, 0, 1)
	if _, _, err := Histogram(a, maxHistBins+1); err == nil {
		t.Fatal("a bin count past the limit was accepted")
	}
	if _, _, err := Histogram(a, 1<<21); err == nil {
		t.Fatal("a 2^21-bin request was accepted")
	}
	counts, edges, err := Histogram(a, 4)
	if err != nil {
		t.Fatalf("a modest histogram: %v", err)
	}
	if counts.Len() != 4 || edges.Len() != 5 {
		t.Fatalf("lengths %d and %d, want 4 and 5", counts.Len(), edges.Len())
	}
	if _, _, _, err := Histogram2D(a, a, maxHistBins+1, 1); err == nil {
		t.Fatal("a per-axis bin count past the limit was accepted")
	}
	// 1025 * 1024 = 1049600 is past the total, each axis alone is not.
	if _, _, _, err := Histogram2D(a, a, 1025, 1024); err == nil {
		t.Fatal("a bin product past the limit was accepted")
	}
	// A wide but legal request goes through: the guard refuses what is
	// past the limit, not what is merely large.
	if _, _, _, err := Histogram2D(a, a, 512, 512); err != nil {
		t.Fatalf("a 512 x 512 histogram: %v", err)
	}
}

// TestNormalQuantileLowerTail pins the P1: continuousQuantile brackets
// upwards from the seed only, and Phi(0) = 0.5 is the infimum it can
// reach, so every q < 0.5 failed to bracket. The mirrored tail must
// reproduce an independent high-precision inverse of Phi.
func TestNormalQuantileLowerTail(t *testing.T) {
	qs := []float64{1e-6, 1e-5, 1e-4, 1e-3, 0.01, 0.025, 0.05, 0.1, 0.25, 0.4, 0.4999}
	worst := 0.0
	for _, q := range qs {
		z, err := NormalQuantile(q)
		if err != nil {
			t.Fatalf("NormalQuantile(%g): %v", q, err)
		}
		if z >= 0 {
			t.Fatalf("NormalQuantile(%g) = %v, want a negative quantile", q, z)
		}
		ref := bigNormalQuantile(q)
		dev := math.Abs(z - ref)
		if dev > worst {
			worst = dev
		}
		if dev > 1e-9 {
			t.Fatalf("NormalQuantile(%g) = %.17g, the 256-bit bisection reference is %.17g (off by %.3g)",
				q, z, ref, dev)
		}
		// The symmetry is exact in the implementation and must hold to
		// the last bit: the mirrored call returns the negated value.
		mirror, err := NormalQuantile(1 - q)
		if err != nil {
			t.Fatalf("NormalQuantile(%g): %v", 1-q, err)
		}
		if z+mirror != 0 {
			t.Fatalf("NormalQuantile(%g) + NormalQuantile(%g) = %.17g, want exactly 0", q, 1-q, z+mirror)
		}
		// And the library's own CDF inverts the result.
		if dev := math.Abs(NormalCDF(z) - q); dev > 1e-12 {
			t.Fatalf("NormalCDF(NormalQuantile(%g)) is off by %.3g", q, dev)
		}
	}
	t.Logf("worst deviation from the 256-bit reference: %.3g", worst)
	// The upper half was already correct and must stay so.
	upper, err := NormalQuantile(0.975)
	if err != nil {
		t.Fatalf("NormalQuantile(0.975): %v", err)
	}
	if math.Abs(upper-1.9599639845400532) > 1e-12 {
		t.Fatalf("NormalQuantile(0.975) = %.17g, want 1.9599639845400532", upper)
	}
}

// TestNormalQuantileNaNRejected pins the P2: NaN passed the written-out
// range guard and every comparison of the bisection returned false, so
// the seed came back as the answer (v = 1 with a nil error).
func TestNormalQuantileNaNRejected(t *testing.T) {
	for _, q := range []float64{math.NaN(), math.Inf(1), math.Inf(-1), -0.5, 1.5, 0, 1} {
		v, err := NormalQuantile(q)
		if err == nil {
			t.Fatalf("NormalQuantile(%v) = %v, want a refusal", q, v)
		}
	}
	if v, err := StudentTQuantile(math.NaN(), 5); err == nil {
		t.Fatalf("StudentTQuantile(NaN, 5) = %v, want a refusal", v)
	}
}

// TestLogisticRegressionSaturatedLikelihoodFinite pins the P1: the
// convergence exit recomputed the sigmoid without the clamp the loop
// relies on, so one far-out covariate saturated Fitted to exactly 0/1
// and turned LogLikelihood into 0·log(0) = NaN while Converged was
// true. The reported values must be the ones the loop maximised.
func TestLogisticRegressionSaturatedLikelihoodFinite(t *testing.T) {
	design, err := core.FromFloats([]float64{
		1, 0, 1, 1, 1, 2, 1, 3, 1, 4, 1, 5, 1, 1e3,
	}, 7, 2)
	if err != nil {
		t.Fatalf("FromFloats: %v", err)
	}
	y := pinVector(t, 0, 1, 0, 1, 1, 0, 1)
	res, err := LogisticRegression(design, y)
	if err != nil {
		t.Fatalf("LogisticRegression: %v", err)
	}
	if !res.Converged {
		t.Fatal("the fit did not converge")
	}
	manual := 0.0
	for i, p := range res.Fitted {
		if p <= 0 || p >= 1 {
			t.Fatalf("Fitted[%d] = %v, want a probability strictly inside (0, 1)", i, p)
		}
		yv := y.FloatAt(i)
		manual += yv*math.Log(p) + (1-yv)*math.Log(1-p)
	}
	if math.IsNaN(res.LogLikelihood) || math.IsInf(res.LogLikelihood, 0) {
		t.Fatalf("LogLikelihood = %v, want the finite value the loop maximised", res.LogLikelihood)
	}
	if math.Abs(res.LogLikelihood-manual) > 1e-12 {
		t.Fatalf("LogLikelihood = %.17g, recomputing from Fitted gives %.17g",
			res.LogLikelihood, manual)
	}
	// The inference must stay finite as well.
	for j := range res.Coefficients {
		mustFinite := []float64{res.Coefficients[j], res.StandardErrors[j],
			res.ZStatistics[j], res.PValues[j]}
		for _, v := range mustFinite {
			if math.IsNaN(v) || math.IsInf(v, 0) {
				t.Fatalf("coefficient %d has a non-finite statistic: %v", j, mustFinite)
			}
		}
	}
}

// TestWeightedLinearRegressionInterceptStatistics pins the P1: the
// intercept column of the unweighted design stops being constant once
// the sqrt weights are applied, so the delegated fit reported the
// no-intercept conventions (uncentred TSS, DModel = p, an F test
// against the zero model) for a regression the caller described with an
// intercept. The statistics must be the weighted-theory ones, against
// the weighted mean.
func TestWeightedLinearRegressionInterceptStatistics(t *testing.T) {
	design, err := core.FromFloats([]float64{1, 0, 1, 1, 1, 2, 1, 3, 1, 4, 1, 5}, 6, 2)
	if err != nil {
		t.Fatalf("FromFloats: %v", err)
	}
	yVals := []float64{10, 10.1, 9.9, 10, 10.2, 9.8}
	y := pinVector(t, yVals...)
	plain, err := LinearRegression(design, y)
	if err != nil {
		t.Fatalf("LinearRegression: %v", err)
	}
	// Uniform weights are the ordinary fit, statistic for statistic.
	uniform, err := WeightedLinearRegression(design, y, pinVector(t, 1, 1, 1, 1, 1, 1))
	if err != nil {
		t.Fatalf("WeightedLinearRegression: %v", err)
	}
	if uniform.DModel != plain.DModel || uniform.RSquared != plain.RSquared ||
		uniform.FStatistic != plain.FStatistic || uniform.FPValue != plain.FPValue {
		t.Fatalf("uniform weights report DModel=%d R2=%.17g F=%.17g p=%.17g, the ordinary fit DModel=%d R2=%.17g F=%.17g p=%.17g",
			uniform.DModel, uniform.RSquared, uniform.FStatistic, uniform.FPValue,
			plain.DModel, plain.RSquared, plain.FStatistic, plain.FPValue)
	}
	for _, w6 := range []float64{1.0000001, 2} {
		w := pinVector(t, 1, 1, 1, 1, 1, w6)
		res, err := WeightedLinearRegression(design, y, w)
		if err != nil {
			t.Fatalf("WeightedLinearRegression(w6=%g): %v", w6, err)
		}
		// The reference from first principles, on the residuals the call
		// itself reports: sumW, the weighted mean, the weighted centred
		// total and the weighted residual sum.
		sumW, sumWY, rss := 0.0, 0.0, 0.0
		weights := []float64{1, 1, 1, 1, 1, w6}
		for i, wi := range weights {
			sumW += wi
			sumWY += wi * yVals[i]
			rss += wi * res.Residuals[i] * res.Residuals[i]
		}
		meanW := sumWY / sumW
		tss := 0.0
		for i, wi := range weights {
			d := yVals[i] - meanW
			tss += wi * d * d
		}
		wantR2 := 1 - rss/tss
		wantF := (tss - rss) / 1 / (rss / 4)
		if res.DModel != 1 {
			t.Fatalf("w6=%g: DModel = %d, want 1 (the design carries the intercept)", w6, res.DModel)
		}
		if math.Abs(res.RSquared-wantR2) > 1e-12 {
			t.Fatalf("w6=%g: R2 = %.17g, the weighted-centred reference is %.17g", w6, res.RSquared, wantR2)
		}
		if math.Abs(res.FStatistic-wantF) > 1e-9*wantF {
			t.Fatalf("w6=%g: F = %.17g, the weighted-centred reference is %.17g", w6, res.FStatistic, wantF)
		}
		// For a single slope F = t^2, so the model p-value is the
		// slope's own two-sided tail.
		if math.Abs(res.FPValue-res.PValues[1]) > 1e-9 {
			t.Fatalf("w6=%g: model p = %.17g, the slope's two-sided tail is %.17g",
				w6, res.FPValue, res.PValues[1])
		}
	}
	// The report's exact acceptance values on this data with w6 = 2, the
	// weighted-centred values (R2 0.172938829787234, F 0.836401640003216,
	// p 0.4121703772248073 from a numerical integration of F(1,4)):
	// the defect reported 0.9998404595, 12534.0045019697 and
	// 2.54531602909285e-08 instead.
	res, err := WeightedLinearRegression(design, y, pinVector(t, 1, 1, 1, 1, 1, 2))
	if err != nil {
		t.Fatalf("WeightedLinearRegression: %v", err)
	}
	if math.Abs(res.RSquared-0.172938829787234) > 1e-9 {
		t.Fatalf("R2 = %.17g, the weighted-centred value is 0.172938829787234", res.RSquared)
	}
	if math.Abs(res.FStatistic-0.836401640003216) > 1e-9 {
		t.Fatalf("F = %.17g, the weighted-centred value is 0.836401640003216", res.FStatistic)
	}
	if math.Abs(res.FPValue-0.4121703772248073) > 1e-9 {
		t.Fatalf("model p = %.17g, the weighted-centred value is 0.4121703772248073", res.FPValue)
	}
	if res.DModel != 1 || res.DResidual != 4 {
		t.Fatalf("degrees of freedom (%d, %d), want (1, 4)", res.DModel, res.DResidual)
	}
	// A single weight a hair away from uniform must not move the
	// statistics: the flip the defect produced was 0.05 to 0.9998.
	hair, err := WeightedLinearRegression(design, y, pinVector(t, 1, 1, 1, 1, 1, 1.0000001))
	if err != nil {
		t.Fatalf("WeightedLinearRegression: %v", err)
	}
	if math.Abs(hair.RSquared-uniform.RSquared) > 1e-6 {
		t.Fatalf("a 1e-7 weight change moved R2 from %.17g to %.17g", uniform.RSquared, hair.RSquared)
	}
	if math.Abs(hair.FPValue-uniform.FPValue) > 1e-6 {
		t.Fatalf("a 1e-7 weight change moved the model p from %.17g to %.17g",
			uniform.FPValue, hair.FPValue)
	}
}

// TestBinomialDrawsRejectsNonPositiveCount pins the P2: the only draw
// function in the file without the n < 1 guard returned a nil array
// with a nil error, which a caller checking only the error dereferences.
func TestBinomialDrawsRejectsNonPositiveCount(t *testing.T) {
	g := core.NewGenerator(1)
	for _, n := range []int{-1, 0} {
		out, err := BinomialDraws(g, n, 4, 0.5)
		if err == nil {
			t.Fatalf("BinomialDraws(n=%d) = %v, want a refusal", n, out)
		}
		if out != nil {
			t.Fatalf("BinomialDraws(n=%d) returned a non-nil array with the refusal", n)
		}
	}
	out, err := BinomialDraws(g, 3, 4, 0.5)
	if err != nil {
		t.Fatalf("BinomialDraws: %v", err)
	}
	if out.Len() != 3 {
		t.Fatalf("length %d, want 3", out.Len())
	}
}

// TestKolmogorovSmirnovRejectsNonFiniteSamples pins an unbounded loop
// found while verifying the package's finiteness policy: the merge walk
// in KolmogorovSmirnovTest advances only past values that compare true
// against the current support point, and a NaN never does, so the walk
// re-reads the same element for ever. A +Inf sample walks through but
// distorts the distance. The sibling tests (MannWhitneyU, ANOVAOneWay)
// refuse a non-finite sample by name and so must this one.
func TestKolmogorovSmirnovRejectsNonFiniteSamples(t *testing.T) {
	clean := pinVector(t, 1, 2, 3, 4, 5)
	nan := pinVector(t, 1, 2, math.NaN(), 4, 5)
	inf := pinVector(t, 1, 2, math.Inf(1), 4, 5)
	for _, tc := range []struct {
		name string
		a, b *core.Array
	}{
		{"NaN first sample", nan, clean},
		{"NaN second sample", clean, nan},
		{"+Inf first sample", inf, clean},
		{"+Inf second sample", clean, inf},
	} {
		t.Run(tc.name, func(t *testing.T) {
			var err error
			pinWatchdog(t, 5*time.Second, "KolmogorovSmirnovTest("+tc.name+")", func() {
				_, _, err = KolmogorovSmirnovTest(tc.a, tc.b)
			})
			if err == nil {
				t.Fatal("a non-finite sample was accepted")
			}
		})
	}
	// The finite path is untouched: identical samples are distance zero
	// with p = 1, and the report's {1,2,3} against {1.5,2.5,3.5} case
	// gives d = 1/3 exactly.
	d, p, err := KolmogorovSmirnovTest(clean, clean)
	if err != nil {
		t.Fatalf("KolmogorovSmirnovTest: %v", err)
	}
	if d != 0 || p != 1 {
		t.Fatalf("identical samples give d = %g, p = %g, want 0 and 1", d, p)
	}
	d, _, err = KolmogorovSmirnovTest(pinVector(t, 1, 2, 3), pinVector(t, 1.5, 2.5, 3.5))
	if err != nil {
		t.Fatalf("KolmogorovSmirnovTest: %v", err)
	}
	// The largest gap is the first support point's 1/3, computed from
	// integer counts: allow the couple of ulps that division costs.
	if math.Abs(d-1.0/3.0) > 1e-15 {
		t.Fatalf("d = %.17g, want 1/3", d)
	}
}

// TestKernelDensityRejectsNonFiniteInput pins the first half of the
// finiteness gap: KernelDensity never scanned the sample or the
// evaluation points, so a NaN came back as a NaN estimate with a nil
// error and a +Inf silently dropped that sample from the average (the
// estimate was then the one of the remaining n-1 points). The clean
// path is pinned against the hand-written KDE definition.
func TestKernelDensityRejectsNonFiniteInput(t *testing.T) {
	sample := pinVector(t, 0, 1, 2, 3)
	points := pinVector(t, 0.5, 2.5)
	for _, tc := range []struct {
		name string
		s, p *core.Array
	}{
		{"NaN sample", pinVector(t, 0, 1, math.NaN(), 3), points},
		{"+Inf sample", pinVector(t, 0, 1, math.Inf(1), 3), points},
		{"NaN point", sample, pinVector(t, 0.5, math.NaN())},
		{"-Inf point", sample, pinVector(t, 0, math.Inf(-1))},
	} {
		t.Run(tc.name, func(t *testing.T) {
			if out, err := KernelDensity(tc.s, 0.5, tc.p); err == nil {
				t.Fatalf("a non-finite input was accepted, out = %v", out)
			}
		})
	}
	// The definition, written out: 1/(n h sqrt(2 pi)) * sum_k
	// exp(-((x - sample[k])/h)^2 / 2) at x = 0.5 and x = 2.5, which are
	// mirror images of each other in this sample and so agree exactly.
	for _, x := range []float64{0.5, 2.5} {
		want := 0.0
		for _, s := range []float64{0, 1, 2, 3} {
			z := (x - s) / 0.5
			want += math.Exp(-0.5 * z * z)
		}
		want /= 4 * 0.5 * math.Sqrt(2*math.Pi)
		out, err := KernelDensity(sample, 0.5, pinVector(t, x))
		if err != nil {
			t.Fatalf("KernelDensity(%g): %v", x, err)
		}
		if got := out.FloatAt(0); got != want {
			t.Fatalf("KernelDensity at %g = %.17g, the hand-written KDE is %.17g", x, got, want)
		}
	}
}

// TestMultivariateNormalRejectsNonFiniteInput pins the second half of
// the finiteness gap: the covariance was scanned for non-finite
// entries, the mean and the point were not, so a NaN mean or point came
// back as a NaN log density and drew NaN vectors, all with a nil error.
func TestMultivariateNormalRejectsNonFiniteInput(t *testing.T) {
	mean := pinVector(t, 0, 0)
	cov, err := core.FromFloats([]float64{1, 0.5, 0.5, 1}, 2, 2)
	if err != nil {
		t.Fatalf("FromFloats: %v", err)
	}
	nanCov, err := core.FromFloats([]float64{1, math.NaN(), math.NaN(), 1}, 2, 2)
	if err != nil {
		t.Fatalf("FromFloats: %v", err)
	}
	x := pinVector(t, 1, 1)
	nanMean := pinVector(t, math.NaN(), 0)
	infMean := pinVector(t, math.Inf(1), 0)
	nanPoint := pinVector(t, 1, math.NaN())
	g := core.NewGenerator(7)
	for _, tc := range []struct {
		name string
		call func() error
	}{
		{"density NaN mean", func() error { _, err := MultivariateNormalLogDensity(nanMean, cov, x); return err }},
		{"density +Inf mean", func() error { _, err := MultivariateNormalLogDensity(infMean, cov, x); return err }},
		{"density NaN covariance", func() error { _, err := MultivariateNormalLogDensity(mean, nanCov, x); return err }},
		{"density NaN point", func() error { _, err := MultivariateNormalLogDensity(mean, cov, nanPoint); return err }},
		{"draws NaN mean", func() error { _, err := MultivariateNormalDraws(g, 2, nanMean, cov); return err }},
		{"draws NaN covariance", func() error { _, err := MultivariateNormalDraws(g, 2, mean, nanCov); return err }},
	} {
		t.Run(tc.name, func(t *testing.T) {
			if err := tc.call(); err == nil {
				t.Fatal("a non-finite input was accepted")
			}
		})
	}
	// The clean path is pinned by hand: with mean 0, cov [[1,.5],[.5,1]]
	// and x = (1,1), the log density is -ln(2 pi) - 0.5 ln(0.75) - 2/3.
	got, err := MultivariateNormalLogDensity(mean, cov, x)
	if err != nil {
		t.Fatalf("MultivariateNormalLogDensity: %v", err)
	}
	want := -math.Log(2*math.Pi) - 0.5*math.Log(0.75) - 2.0/3.0
	if math.Abs(got-want) > 1e-12 {
		t.Fatalf("log density = %.17g, hand value %.17g", got, want)
	}
	draws, err := MultivariateNormalDraws(g, 3, mean, cov)
	if err != nil {
		t.Fatalf("MultivariateNormalDraws: %v", err)
	}
	if draws.Len() != 6 {
		t.Fatalf("draws length %d, want 6", draws.Len())
	}
}

// TestMultivariateNormalRejectsAsymmetricCovariance pins the P2: the
// Cholesky factor read only the lower triangle, so an asymmetric matrix
// was silently symmetrised and the reported density described a
// different distribution from the one handed in. A mirror that differs
// beyond a relative tolerance is refused by name; a mirror that differs
// by rounding only is still accepted.
func TestMultivariateNormalRejectsAsymmetricCovariance(t *testing.T) {
	mean := pinVector(t, 0, 0)
	x := pinVector(t, 1, 1)
	asym, err := core.FromFloats([]float64{1, 100, 0.5, 1}, 2, 2)
	if err != nil {
		t.Fatalf("FromFloats: %v", err)
	}
	for _, tc := range []struct {
		name string
		call func() error
	}{
		{"density", func() error { _, err := MultivariateNormalLogDensity(mean, asym, x); return err }},
		{"draws", func() error { _, err := MultivariateNormalDraws(core.NewGenerator(3), 2, mean, asym); return err }},
	} {
		t.Run(tc.name, func(t *testing.T) {
			if err := tc.call(); err == nil {
				t.Fatal("the asymmetric covariance was accepted")
			}
		})
	}
	// Rounding-level asymmetry is not asymmetry: 0.5 + 1e-15 against 0.5
	// is a mirror pair an A*A^T assembly can produce.
	sym, err := core.FromFloats([]float64{1, 0.5 + 1e-15, 0.5, 1}, 2, 2)
	if err != nil {
		t.Fatalf("FromFloats: %v", err)
	}
	got, err := MultivariateNormalLogDensity(mean, sym, x)
	if err != nil {
		t.Fatalf("a rounding-level mirror difference was refused: %v", err)
	}
	want := -math.Log(2*math.Pi) - 0.5*math.Log(0.75) - 2.0/3.0
	if math.Abs(got-want) > 1e-12 {
		t.Fatalf("log density = %.17g, the symmetric hand value is %.17g", got, want)
	}
}

// TestKolmogorovSeriesTruncationRule pins the inert defect: the early
// exit asked whether the current term was 1e-14 times the next, which is
// exp(2(2k+1)l^2) > 1 for every l > 0, so it could never fire and the
// comment described a termination check that was not there. The
// truncation is now the first term below the rounding level; the count
// must be exactly that term, and the summed series must agree with a
// directly summed one to the last bit.
func TestKolmogorovSeriesTruncationRule(t *testing.T) {
	term := func(k int, lambda float64) float64 {
		return 2 * math.Exp(-2*float64(k)*float64(k)*lambda*lambda)
	}
	for _, lambda := range []float64{0.2, 0.3, 0.5, 0.75, 1, 1.5, 2, 3, 5} {
		n := kolmogorovTermCount(lambda)
		if n < 1 || n > maxKolmogorovTerms {
			t.Fatalf("lambda = %g: %d terms is outside [1, %d]", lambda, n, maxKolmogorovTerms)
		}
		if term(n, lambda) >= kolmogorovEps {
			t.Fatalf("lambda = %g: term %d is %g, still above the %g rounding level: the series stops too early",
				lambda, n, term(n, lambda), kolmogorovEps)
		}
		if n > 1 && term(n-1, lambda) < kolmogorovEps {
			t.Fatalf("lambda = %g: term %d is already below the rounding level: the series adds terms that cannot move the sum",
				lambda, n-1)
		}
		direct := 0.0
		for k := 1; k <= 400; k++ {
			if k%2 == 0 {
				direct -= term(k, lambda)
			} else {
				direct += term(k, lambda)
			}
		}
		want := min(1, max(0, direct))
		if got := kolmogorovTail(lambda); got != want {
			t.Fatalf("lambda = %g: tail = %.17g, a directly summed 400-term series gives %.17g", lambda, got, want)
		}
	}
	// The documented short circuit below lambda = 0.2 is untouched.
	if got := kolmogorovTail(0.1); got != 1 {
		t.Fatalf("kolmogorovTail(0.1) = %g, want 1", got)
	}
}

// bootstrapMean is the bootstrap statistic the sweeps below resample with.
func bootstrapMean(a *core.Array) (float64, error) {
	return core.Mean(a)
}

// TestNaNParametersRejected pins the class sweep over the package's
// float parameters. Each guard used to be written out (`x < 0 || x > 1`
// and friends), which a NaN walks through because every comparison
// against a NaN is false: the parameter reached the algorithm and the
// call answered a value, a NaN, or nothing at all. The NaN-rejecting
// form `!(x >= 0 && x <= 1)` must turn every one of them into the
// refusal the caller deserves. The watchdog turns a regression to the
// unbounded discrete-quantile loop into a failure rather than a stuck
// test binary.
func TestNaNParametersRejected(t *testing.T) {
	g := core.NewGenerator(1)
	clean := pinVector(t, 1, 2, 3, 4, 5)
	nanSample := pinVector(t, 1, 2, math.NaN(), 4, 5)
	cases := []struct {
		name string
		call func() error
	}{
		{"GammaLower(NaN shape)", func() error { _, err := GammaLower(math.NaN(), 1); return err }},
		{"GammaLower(NaN x)", func() error { _, err := GammaLower(1, math.NaN()); return err }},
		{"GammaUpper(NaN shape)", func() error { _, err := GammaUpper(math.NaN(), 1); return err }},
		{"GammaUpper(NaN x)", func() error { _, err := GammaUpper(1, math.NaN()); return err }},
		{"BetaIncomplete(NaN x)", func() error { _, err := BetaIncomplete(math.NaN(), 2, 2); return err }},
		{"BetaIncomplete(NaN a)", func() error { _, err := BetaIncomplete(0.5, math.NaN(), 2); return err }},
		{"BetaIncomplete(NaN b)", func() error { _, err := BetaIncomplete(0.5, 2, math.NaN()); return err }},
		{"ExponentialCDF(NaN x)", func() error { _, err := ExponentialCDF(math.NaN(), 2); return err }},
		{"ExponentialCDF(NaN rate)", func() error { _, err := ExponentialCDF(1, math.NaN()); return err }},
		{"GammaCDF(NaN x)", func() error { _, err := GammaCDF(math.NaN(), 2, 3); return err }},
		{"GammaCDF(NaN shape)", func() error { _, err := GammaCDF(1, math.NaN(), 3); return err }},
		{"GammaCDF(NaN rate)", func() error { _, err := GammaCDF(1, 2, math.NaN()); return err }},
		{"ChiSquareCDF(NaN x)", func() error { _, err := ChiSquareCDF(math.NaN(), 2); return err }},
		{"StudentTCDF(NaN t)", func() error { _, err := StudentTCDF(math.NaN(), 5); return err }},
		{"PoissonCDF(NaN lambda)", func() error { _, err := PoissonCDF(3, math.NaN()); return err }},
		{"BinomialCDF(NaN p)", func() error { _, err := BinomialCDF(3, 10, math.NaN()); return err }},
		{"ExponentialQuantile(NaN q)", func() error { _, err := ExponentialQuantile(math.NaN(), 2); return err }},
		{"ExponentialQuantile(NaN rate)", func() error { _, err := ExponentialQuantile(0.5, math.NaN()); return err }},
		{"GammaQuantile(NaN q)", func() error { _, err := GammaQuantile(math.NaN(), 2, 3); return err }},
		{"GammaQuantile(NaN shape)", func() error { _, err := GammaQuantile(0.5, math.NaN(), 3); return err }},
		{"GammaQuantile(NaN rate)", func() error { _, err := GammaQuantile(0.5, 2, math.NaN()); return err }},
		{"ChiSquareQuantile(NaN q)", func() error { _, err := ChiSquareQuantile(math.NaN(), 2); return err }},
		{"StudentTQuantile(NaN q)", func() error { _, err := StudentTQuantile(math.NaN(), 5); return err }},
		{"PoissonQuantile(NaN q)", func() error { _, err := PoissonQuantile(math.NaN(), 3); return err }},
		{"PoissonQuantile(NaN lambda)", func() error { _, err := PoissonQuantile(0.5, math.NaN()); return err }},
		{"BinomialQuantile(NaN q)", func() error { _, err := BinomialQuantile(math.NaN(), 0.5, 5); return err }},
		{"BinomialQuantile(NaN p)", func() error { _, err := BinomialQuantile(0.5, math.NaN(), 5); return err }},
		{"ExponentialDraws(NaN rate)", func() error { _, err := ExponentialDraws(g, 3, math.NaN()); return err }},
		{"GammaDraws(NaN shape)", func() error { _, err := GammaDraws(g, 3, math.NaN(), 2); return err }},
		{"GammaDraws(NaN rate)", func() error { _, err := GammaDraws(g, 3, 2, math.NaN()); return err }},
		{"PoissonDraws(NaN lambda)", func() error { _, err := PoissonDraws(g, 3, math.NaN()); return err }},
		{"BinomialDraws(NaN p)", func() error { _, err := BinomialDraws(g, 3, 4, math.NaN()); return err }},
		{"ChiSquareGoodnessOfFit(NaN observed)", func() error {
			_, _, _, err := ChiSquareGoodnessOfFit(nanSample, clean)
			return err
		}},
		{"ChiSquareGoodnessOfFit(NaN expected)", func() error {
			_, _, _, err := ChiSquareGoodnessOfFit(clean, nanSample)
			return err
		}},
		{"BootstrapCI(NaN level)", func() error { _, _, err := BootstrapCI(clean, bootstrapMean, math.NaN(), 20, 1); return err }},
		{"TrimmedMean(NaN fraction)", func() error { _, err := TrimmedMean(clean, math.NaN()); return err }},
		{"KernelDensity(NaN bandwidth)", func() error { _, err := KernelDensity(clean, math.NaN(), clean); return err }},
		{"WelchTTest(NaN sample)", func() error { _, _, _, err := WelchTTest(nanSample, clean); return err }},
	}
	for _, tc := range cases {
		t.Run(tc.name, func(t *testing.T) {
			var err error
			pinWatchdog(t, 5*time.Second, tc.name, func() { err = tc.call() })
			if err == nil {
				t.Fatalf("%s answered with no error for a NaN parameter", tc.name)
			}
		})
	}
	// No guard overreaches: the legal domain still answers, and the two
	// values below are hand-checkable. The incomplete beta is the
	// continued fraction's own rounding level, a few ulps at 0.5.
	if p, err := BetaIncomplete(0.5, 2, 2); err != nil || math.Abs(p-0.5) > 1e-15 {
		t.Fatalf("BetaIncomplete(0.5, 2, 2) = %v, %v, want 0.5 by symmetry", p, err)
	}
	if p, err := GammaLower(1, 1); err != nil || math.Abs(p-(1-math.Exp(-1))) > 1e-15 {
		t.Fatalf("GammaLower(1, 1) = %v, %v, want 1 - e^-1", p, err)
	}
	if q, err := ExponentialQuantile(0.5, 3); err != nil || math.Abs(q-math.Ln2/3) > 1e-15 {
		t.Fatalf("ExponentialQuantile(0.5, 3) = %v, %v, want ln 2 / 3", q, err)
	}
	if p, err := PoissonQuantile(0.5, 3); err != nil || p != 3 {
		t.Fatalf("PoissonQuantile(0.5, 3) = %v, %v, want 3", p, err)
	}
}

// TestComplexInputsRejected pins the dtype sweep: every entry point in
// the package refuses a complex array by name, and the ones exercised
// here (the three MVN inputs, BootstrapCI's data, KernelDensity's
// sample and points) walked the float accessor into a nil payload and
// panicked instead.
func TestComplexInputsRejected(t *testing.T) {
	cx2, err := core.FromComplexes([]complex128{1 + 1i, 2 - 1i}, 2)
	if err != nil {
		t.Fatalf("FromComplexes: %v", err)
	}
	cx5, err := core.FromComplexes([]complex128{1 + 1i, 2, 3, 4, 5}, 5)
	if err != nil {
		t.Fatalf("FromComplexes: %v", err)
	}
	cxMat, err := core.FromComplexes([]complex128{1, 0, 0, 1}, 2, 2)
	if err != nil {
		t.Fatalf("FromComplexes: %v", err)
	}
	mean := pinVector(t, 0, 0)
	x := pinVector(t, 1, 1)
	cov, err := core.FromFloats([]float64{1, 0.5, 0.5, 1}, 2, 2)
	if err != nil {
		t.Fatalf("FromFloats: %v", err)
	}
	g := core.NewGenerator(3)
	for _, tc := range []struct {
		name string
		call func() error
	}{
		{"MVN density complex mean", func() error { _, err := MultivariateNormalLogDensity(cx2, cov, x); return err }},
		{"MVN density complex covariance", func() error { _, err := MultivariateNormalLogDensity(mean, cxMat, x); return err }},
		{"MVN density complex point", func() error { _, err := MultivariateNormalLogDensity(mean, cov, cx2); return err }},
		{"MVN draws complex mean", func() error { _, err := MultivariateNormalDraws(g, 2, cx2, cov); return err }},
		{"MVN draws complex covariance", func() error { _, err := MultivariateNormalDraws(g, 2, mean, cxMat); return err }},
		{"BootstrapCI complex data", func() error { _, _, err := BootstrapCI(cx5, bootstrapMean, 0.95, 20, 1); return err }},
		{"KernelDensity complex sample", func() error { _, err := KernelDensity(cx5, 0.5, cx5); return err }},
		{"KernelDensity complex points", func() error { _, err := KernelDensity(pinVector(t, 1, 2, 3, 4, 5), 0.5, cx5); return err }},
		{"Histogram2D complex sample", func() error { _, _, _, err := Histogram2D(cx5, cx5, 2, 2); return err }},
	} {
		t.Run(tc.name, func(t *testing.T) {
			if err := tc.call(); err == nil {
				t.Fatal("a complex input was accepted")
			}
		})
	}
	// The real path is untouched: a constant sample has a degenerate
	// bootstrap interval at its own value.
	lo, hi, err := BootstrapCI(pinVector(t, 5, 5, 5, 5), bootstrapMean, 0.95, 50, 1)
	if err != nil {
		t.Fatalf("BootstrapCI: %v", err)
	}
	if lo != 5 || hi != 5 {
		t.Fatalf("the constant sample gives [%g, %g], want [5, 5]", lo, hi)
	}
}

// TestLinearRegressionConstantResponseStatistics pins the NaN guard a
// degenerate response needs: with no variation in y the fit reproduces
// it exactly, so the explained and residual sums are both zero and the
// F statistic is 0/0. A NaN there would spread through every consumer
// of the result, where the finite zero is the value the tail reads as
// p = 1: there is no evidence of a model.
func TestLinearRegressionConstantResponseStatistics(t *testing.T) {
	design, err := core.FromFloats([]float64{1, 0, 1, 1, 1, 2, 1, 3}, 4, 2)
	if err != nil {
		t.Fatalf("FromFloats: %v", err)
	}
	y := pinVector(t, 5, 5, 5, 5)
	for _, tc := range []struct {
		name string
		fit  func() (*LinearRegressionResult, error)
	}{
		{"ordinary", func() (*LinearRegressionResult, error) { return LinearRegression(design, y) }},
		{"weighted", func() (*LinearRegressionResult, error) {
			return WeightedLinearRegression(design, y, pinVector(t, 1, 1, 1, 2))
		}},
	} {
		t.Run(tc.name, func(t *testing.T) {
			res, err := tc.fit()
			if err != nil {
				t.Fatalf("the constant-response fit: %v", err)
			}
			if math.IsNaN(res.FStatistic) || res.FStatistic != 0 {
				t.Fatalf("FStatistic = %v, want the finite 0", res.FStatistic)
			}
			if res.FPValue != 1 {
				t.Fatalf("FPValue = %v, want 1", res.FPValue)
			}
			if res.DModel != 1 || res.DResidual != 2 {
				t.Fatalf("degrees of freedom (%d, %d), want (1, 2)", res.DModel, res.DResidual)
			}
		})
	}
}

// bigPhi returns Phi(x) in 256-bit arithmetic through the Taylor
// series of the error function at x/sqrt(2). Nothing here calls the
// library, so it is an independent reference for the quantile tests.
func bigPhi(x *big.Float) *big.Float {
	const prec = 256
	one := new(big.Float).SetPrec(prec).SetInt64(1)
	y := new(big.Float).SetPrec(prec).Quo(new(big.Float).SetPrec(prec).Set(x), bigSqrt2())
	neg := y.Sign() < 0
	ay := new(big.Float).SetPrec(prec).Abs(y)
	// erf(ay) = (2/sqrt(pi)) * sum_{n>=0} (-1)^n ay^(2n+1) / (n! (2n+1)).
	term := new(big.Float).SetPrec(prec).Set(ay) // ay^(2n+1)/n! at n = 0
	sum := new(big.Float).SetPrec(prec).Set(ay)
	ay2 := new(big.Float).SetPrec(prec).Mul(ay, ay)
	for n := 1; ; n++ {
		term.Mul(term, ay2)
		term.Quo(term, new(big.Float).SetPrec(prec).SetInt64(int64(n)))
		term.Neg(term)
		contrib := new(big.Float).SetPrec(prec).Quo(term,
			new(big.Float).SetPrec(prec).SetInt64(int64(2*n+1)))
		sum.Add(sum, contrib)
		if contrib.Sign() == 0 || contrib.MantExp(nil) < -400 {
			break
		}
	}
	// 2/sqrt(pi).
	twoOverSqrtPi := new(big.Float).SetPrec(prec).Quo(
		new(big.Float).SetPrec(prec).SetInt64(2),
		new(big.Float).SetPrec(prec).Sqrt(new(big.Float).SetPrec(prec).SetFloat64(math.Pi)))
	erf := new(big.Float).SetPrec(prec).Mul(sum, twoOverSqrtPi)
	if neg {
		erf.Neg(erf)
	}
	// Phi(x) = (1 + erf(x/sqrt(2)))/2.
	return new(big.Float).SetPrec(prec).Quo(
		new(big.Float).SetPrec(prec).Add(one, erf),
		new(big.Float).SetPrec(prec).SetInt64(2))
}

// bigNormalQuantile inverts bigPhi by bisection on [-12, 12], the
// independent reference the NormalQuantile test compares against.
func bigNormalQuantile(q float64) float64 {
	const prec = 256
	target := new(big.Float).SetPrec(prec).SetFloat64(q)
	lo := new(big.Float).SetPrec(prec).SetFloat64(-12)
	hi := new(big.Float).SetPrec(prec).SetFloat64(12)
	mid := new(big.Float).SetPrec(prec)
	for range 220 {
		mid.Add(lo, hi)
		mid.Quo(mid, new(big.Float).SetPrec(prec).SetInt64(2))
		if bigPhi(mid).Cmp(target) < 0 {
			lo.Set(mid)
		} else {
			hi.Set(mid)
		}
	}
	mid.Add(lo, hi)
	mid.Quo(mid, new(big.Float).SetPrec(prec).SetInt64(2))
	out, _ := mid.Float64()
	return out
}