384 lines
11 KiB
Go
384 lines
11 KiB
Go
// 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")
|
||
}
|
||
}
|