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

package core

import (
	"math"
	"testing"
)

func mustFloats(t *testing.T, vals []float64, shape ...int) *Array {
	t.Helper()
	if len(shape) == 0 {
		shape = []int{len(vals)}
	}
	a, err := FromFloats(vals, shape...)
	if err != nil {
		t.Fatalf("FromFloats: %v", err)
	}
	return a
}

// TestGammaFamily checks Γ, ln|Γ| and B against exact values.
func TestGammaFamily(t *testing.T) {
	x := mustFloats(t, []float64{1, 5, 0.5, -0.5})
	g, err := Gamma(x)
	if err != nil {
		t.Fatalf("Gamma: %v", err)
	}
	// Γ(1)=1, Γ(5)=24, Γ(½)=√π, Γ(−½)=−2√π.
	want := []float64{1, 24, math.SqrtPi, -2 * math.SqrtPi}
	for i := range want {
		if math.Abs(g.FloatAt(i)-want[i]) > 1e-12*(1+math.Abs(want[i])) {
			t.Fatalf("Gamma(%v) = %v, want %v", x.FloatAt(i), g.FloatAt(i), want[i])
		}
	}

	lg, err := LnGamma(x)
	if err != nil {
		t.Fatalf("LnGamma: %v", err)
	}
	if math.Abs(lg.FloatAt(1)-math.Log(24)) > 1e-12 {
		t.Fatalf("LnGamma(5) = %v, want ln 24", lg.FloatAt(1))
	}

	bx := mustFloats(t, []float64{1, 2})
	by := mustFloats(t, []float64{1, 3})
	b, err := Beta(bx, by)
	if err != nil {
		t.Fatalf("Beta: %v", err)
	}
	// B(1,1) = 1 and B(2,3) = Γ(2)Γ(3)/Γ(5) = 2/24 = 1/12.
	if math.Abs(b.FloatAt(0)-1) > 1e-12 || math.Abs(b.FloatAt(1)-1.0/12) > 1e-12 {
		t.Fatalf("Beta = [%v, %v], want [1, 1/12]", b.FloatAt(0), b.FloatAt(1))
	}
	if _, err := Beta(mustFloats(t, []float64{1}), by); err == nil {
		t.Fatal("expected a shape-mismatch error for Beta")
	}
}

// TestErrorFunctions checks the error function family through the
// defining identity erf(x) + erfc(x) = 1 and known values.
func TestErrorFunctions(t *testing.T) {
	x := mustFloats(t, []float64{0, 0.5, 1, 3})
	e, err := Erf(x)
	if err != nil {
		t.Fatalf("Erf: %v", err)
	}
	c, err := Erfc(x)
	if err != nil {
		t.Fatalf("Erfc: %v", err)
	}
	for i := range x.Len() {
		if math.Abs(e.FloatAt(i)+c.FloatAt(i)-1) > 1e-12 {
			t.Fatalf("erf(%v)+erfc(%v) = %v, want 1",
				x.FloatAt(i), x.FloatAt(i), e.FloatAt(i)+c.FloatAt(i))
		}
	}
	if e.FloatAt(0) != 0 {
		t.Fatalf("erf(0) = %v, want 0", e.FloatAt(0))
	}
	if math.Abs(e.FloatAt(2)-0.8427007929497149) > 1e-12 {
		t.Fatalf("erf(1) = %v, want 0.8427007929497149", e.FloatAt(2))
	}
}

// TestSinc checks the normalised sinc: value 1 at the origin, zeros at
// non-zero integers, 2/π at a half.
func TestSinc(t *testing.T) {
	x := mustFloats(t, []float64{0, 0.5, 1, 2})
	s, err := Sinc(x)
	if err != nil {
		t.Fatalf("Sinc: %v", err)
	}
	if s.FloatAt(0) != 1 {
		t.Fatalf("Sinc(0) = %v, want 1", s.FloatAt(0))
	}
	if math.Abs(s.FloatAt(1)-2/math.Pi) > 1e-12 {
		t.Fatalf("Sinc(0.5) = %v, want 2/π", s.FloatAt(1))
	}
	for _, i := range []int{2, 3} {
		if math.Abs(s.FloatAt(i)) > 1e-15 {
			t.Fatalf("Sinc(%v) = %v, want 0", x.FloatAt(i), s.FloatAt(i))
		}
	}
}

// TestLegendrePolynomials checks P_l against closed forms and the
// Bonnet recurrence across degrees.
func TestLegendrePolynomials(t *testing.T) {
	x := mustFloats(t, []float64{-0.8, -0.3, 0, 0.42, 0.9})
	closed := []struct {
		l    int
		want func(v float64) float64
	}{
		{0, func(v float64) float64 { return 1 }},
		{1, func(v float64) float64 { return v }},
		{2, func(v float64) float64 { return 0.5 * (3*v*v - 1) }},
		{3, func(v float64) float64 { return 0.5 * (5*v*v*v - 3*v) }},
		// P_7's explicit expansion: (429v⁷ − 693v⁵ + 315v³ − 35v)/16.
		{7, func(v float64) float64 {
			return (429*math.Pow(v, 7) - 693*math.Pow(v, 5) + 315*v*v*v - 35*v) / 16
		}},
	}
	for _, tc := range closed {
		p, err := Legendre(tc.l, x)
		if err != nil {
			t.Fatalf("Legendre(%d): %v", tc.l, err)
		}
		for i := range x.Len() {
			v := x.FloatAt(i)
			want := tc.want(v)
			if math.Abs(p.FloatAt(i)-want) > 1e-12*(1+math.Abs(want)) {
				t.Fatalf("P_%d(%v) = %v, want %v", tc.l, v, p.FloatAt(i), want)
			}
		}
	}
	if _, err := Legendre(-1, x); err == nil {
		t.Fatal("expected an error for a negative degree")
	}
}

// TestLegendreAssociated checks the associated functions against
// closed forms, the negative-order relation and orthogonality.
func TestLegendreAssociated(t *testing.T) {
	x := mustFloats(t, []float64{-0.6, 0, 0.35, 0.77})
	// P_1^1 = −√(1−x²) with the Condon-Shortley phase.
	p11, err := LegendreAssociated(1, 1, x)
	if err != nil {
		t.Fatalf("LegendreAssociated(1,1): %v", err)
	}
	for i := range x.Len() {
		v := x.FloatAt(i)
		want := -math.Sqrt(1 - v*v)
		if math.Abs(p11.FloatAt(i)-want) > 1e-12 {
			t.Fatalf("P_1^1(%v) = %v, want %v", v, p11.FloatAt(i), want)
		}
	}
	// P_2^1 = −3x√(1−x²).
	p21, err := LegendreAssociated(2, 1, x)
	if err != nil {
		t.Fatalf("LegendreAssociated(2,1): %v", err)
	}
	for i := range x.Len() {
		v := x.FloatAt(i)
		want := -3 * v * math.Sqrt(1-v*v)
		if math.Abs(p21.FloatAt(i)-want) > 1e-12 {
			t.Fatalf("P_2^1(%v) = %v, want %v", v, p21.FloatAt(i), want)
		}
	}
	// P_1^0 must equal P_1.
	p10, err := LegendreAssociated(1, 0, x)
	if err != nil {
		t.Fatalf("LegendreAssociated(1,0): %v", err)
	}
	p1, _ := Legendre(1, x)
	for i := range x.Len() {
		if p10.FloatAt(i) != p1.FloatAt(i) {
			t.Fatalf("P_1^0(%v) = %v, want P_1 = %v", x.FloatAt(i), p10.FloatAt(i), p1.FloatAt(i))
		}
	}
	// Negative order: P_1^(−1) = −(0)!/(2)!·P_1^1 = −½ P_1^1.
	pm1, err := LegendreAssociated(1, -1, x)
	if err != nil {
		t.Fatalf("LegendreAssociated(1,−1): %v", err)
	}
	for i := range x.Len() {
		if math.Abs(pm1.FloatAt(i)+0.5*p11.FloatAt(i)) > 1e-12 {
			t.Fatalf("P_1^−1(%v) = %v, want %v", x.FloatAt(i), pm1.FloatAt(i), -0.5*p11.FloatAt(i))
		}
	}
	// Orthogonality on a coarse Simpson grid: ∫₋₁¹ P_2 P_3 dx = 0.
	const n = 2001
	grid := make([]float64, n)
	for i := range n {
		grid[i] = -1 + 2*float64(i)/float64(n-1)
	}
	g := mustFloats(t, grid, n)
	p2, _ := Legendre(2, g)
	p3, _ := Legendre(3, g)
	sum := 0.0
	for i := range n {
		w := 1.0
		if i == 0 || i == n-1 {
			w = 1
		} else if i%2 == 1 {
			w = 4
		} else {
			w = 2
		}
		sum += w * p2.FloatAt(i) * p3.FloatAt(i)
	}
	sum *= 2.0 / 3.0 / float64(n-1)
	if math.Abs(sum) > 1e-6 {
		t.Fatalf("∫P_2 P_3 = %v, want 0", sum)
	}
	if _, err := LegendreAssociated(2, 3, x); err == nil {
		t.Fatal("expected an error for |m| > l")
	}
}

// TestSphericalHarmonics checks Y_lm against closed forms and the
// conjugation symmetry Y_{l,−m} = (−1)^m conj(Y_{l,m}).
func TestSphericalHarmonics(t *testing.T) {
	theta := mustFloats(t, []float64{0.3, 1.1, 2.4})
	phi := mustFloats(t, []float64{0, 0.9, -1.3})

	// Y_0^0 = 1/√(4π).
	y00, err := SphericalHarmonic(0, 0, theta, phi)
	if err != nil {
		t.Fatalf("SphericalHarmonic(0,0): %v", err)
	}
	want := 1 / math.Sqrt(4*math.Pi)
	for i := range theta.Len() {
		if d := cmplxAbsDiff(y00.ComplexAt(i), complex(want, 0)); d > 1e-12 {
			t.Fatalf("Y_0^0 = %v, want %v", y00.ComplexAt(i), want)
		}
	}

	// Y_1^0 = √(3/4π) cosθ.
	y10, err := SphericalHarmonic(1, 0, theta, phi)
	if err != nil {
		t.Fatalf("SphericalHarmonic(1,0): %v", err)
	}
	for i := range theta.Len() {
		w := math.Sqrt(3/(4*math.Pi)) * math.Cos(theta.FloatAt(i))
		if d := cmplxAbsDiff(y10.ComplexAt(i), complex(w, 0)); d > 1e-12 {
			t.Fatalf("Y_1^0 = %v, want %v", y10.ComplexAt(i), w)
		}
	}

	// Y_1^1 = −√(3/8π) sinθ e^{iφ}.
	y11, err := SphericalHarmonic(1, 1, theta, phi)
	if err != nil {
		t.Fatalf("SphericalHarmonic(1,1): %v", err)
	}
	for i := range theta.Len() {
		mod := math.Sqrt(3/(8*math.Pi)) * math.Sin(theta.FloatAt(i))
		w := -complex(mod, 0) * cmplxExpStd(complex(0, phi.FloatAt(i)))
		if d := cmplxAbsDiff(y11.ComplexAt(i), w); d > 1e-12 {
			t.Fatalf("Y_1^1 = %v, want %v", y11.ComplexAt(i), w)
		}
	}

	// Conjugation symmetry.
	ym1, err := SphericalHarmonic(1, -1, theta, phi)
	if err != nil {
		t.Fatalf("SphericalHarmonic(1,−1): %v", err)
	}
	for i := range theta.Len() {
		if d := cmplxAbsDiff(ym1.ComplexAt(i), -cmplxConjStd(y11.ComplexAt(i))); d > 1e-12 {
			t.Fatalf("Y_1^−1 symmetry broken: %v vs %v",
				ym1.ComplexAt(i), -cmplxConjStd(y11.ComplexAt(i)))
		}
	}

	if _, err := SphericalHarmonic(1, 2, theta, phi); err == nil {
		t.Fatal("expected an error for |m| > l")
	}
	if _, err := SphericalHarmonic(0, 0, mustFloats(t, []float64{1}), phi); err == nil {
		t.Fatal("expected a shape-mismatch error")
	}
}

// cmplxAbsDiff returns |a − b| for two complex values.
func cmplxAbsDiff(a, b complex128) float64 {
	d := real(a) - real(b)
	e := imag(a) - imag(b)
	return math.Sqrt(d*d + e*e)
}

// cmplxExpStd is e^{i·z}.
func cmplxExpStd(z complex128) complex128 {
	return complex(math.Cos(imag(z)), math.Sin(imag(z)))
}

// cmplxConjStd conjugates z.
func cmplxConjStd(z complex128) complex128 {
	return complex(real(z), -imag(z))
}

// TestSphericalBessel checks j_l and y_l against their closed forms
// for l ≤ 2, the recurrence linking three consecutive degrees, and
// the small-x limits of j.
func TestSphericalBessel(t *testing.T) {
	x := mustFloats(t, []float64{0.05, 0.5, 1, 3.7, 20})

	j0, err := SphericalBesselJ(0, x)
	if err != nil {
		t.Fatalf("SphericalBesselJ(0): %v", err)
	}
	y0, err := SphericalBesselY(0, x)
	if err != nil {
		t.Fatalf("SphericalBesselY(0): %v", err)
	}
	for i := range x.Len() {
		xv := x.FloatAt(i)
		if math.Abs(j0.FloatAt(i)-math.Sin(xv)/xv) > 1e-12 {
			t.Fatalf("j_0(%v) = %v, want %v", xv, j0.FloatAt(i), math.Sin(xv)/xv)
		}
		if math.Abs(y0.FloatAt(i)-(-math.Cos(xv)/xv)) > 1e-12 {
			t.Fatalf("y_0(%v) = %v, want %v", xv, y0.FloatAt(i), -math.Cos(xv)/xv)
		}
	}

	j1, _ := SphericalBesselJ(1, x)
	j2, _ := SphericalBesselJ(2, x)
	y1, _ := SphericalBesselY(1, x)
	y2, _ := SphericalBesselY(2, x)
	for i := range x.Len() {
		xv := x.FloatAt(i)
		// j_1 = sin/x² − cos/x and y_1 = −cos/x² − sin/x.
		j1w := math.Sin(xv)/(xv*xv) - math.Cos(xv)/xv
		y1w := -math.Cos(xv)/(xv*xv) - math.Sin(xv)/xv
		if math.Abs(j1.FloatAt(i)-j1w) > 1e-11*(1+math.Abs(j1w)) {
			t.Fatalf("j_1(%v) = %v, want %v", xv, j1.FloatAt(i), j1w)
		}
		if math.Abs(y1.FloatAt(i)-y1w) > 1e-11*(1+math.Abs(y1w)) {
			t.Fatalf("y_1(%v) = %v, want %v", xv, y1.FloatAt(i), y1w)
		}
		// Three-term recurrence j_{l+1} = (2l+1)/x·j_l − j_{l−1} at l=1.
		if d := j2.FloatAt(i) - (3/xv*j1.FloatAt(i) - j0.FloatAt(i)); math.Abs(d) > 1e-10*(1+math.Abs(j2.FloatAt(i))) {
			t.Fatalf("j recurrence at %v broken: %v", xv, d)
		}
		if d := y2.FloatAt(i) - (3/xv*y1.FloatAt(i) - y0.FloatAt(i)); math.Abs(d) > 1e-10*(1+math.Abs(y2.FloatAt(i))) {
			t.Fatalf("y recurrence at %v broken: %v", xv, d)
		}
	}

	// Small-x limits: j_l(0) = δ_{l0}.
	atZero := mustFloats(t, []float64{0})
	for l := range 4 {
		jl, err := SphericalBesselJ(l, atZero)
		if err != nil {
			t.Fatalf("SphericalBesselJ(%d, 0): %v", l, err)
		}
		want := 0.0
		if l == 0 {
			want = 1
		}
		if jl.FloatAt(0) != want {
			t.Fatalf("j_%d(0) = %v, want %v", l, jl.FloatAt(0), want)
		}
	}
	// j_l below its degree is tiny compared to j_0 at the same point:
	// the downward recurrence must keep the small roots, not flood
	// them with round-off. Reference: the small-x asymptote
	// j_l(x) tends to x^l/(2l+1)!!; at x = 0.01 the first correction
	// term x²/(2(2l+3)) is 6.4e-7 relative, well inside the tolerance.
	small := mustFloats(t, []float64{0.01})
	j5, err := SphericalBesselJ(5, small)
	if err != nil {
		t.Fatalf("SphericalBesselJ(5): %v", err)
	}
	asymptote := math.Pow(0.01, 5) / 10395 // 11!! = 10395
	if d := math.Abs(j5.FloatAt(0) - asymptote); d > 1e-5*math.Abs(asymptote) {
		t.Fatalf("j_5(0.01) = %v, want ≈ %v (downward recurrence lost the root)",
			j5.FloatAt(0), asymptote)
	}
	if _, err := SphericalBesselJ(-1, small); err == nil {
		t.Fatal("expected an error for a negative degree")
	}
}