294 lines
9.6 KiB
Go
294 lines
9.6 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
package core
|
||
|
||
import (
|
||
"math"
|
||
"testing"
|
||
)
|
||
|
||
// simpsonRef integrates f over [0, phi] with a dense composite
|
||
// Simpson rule: the independent oracle for the elliptic values.
|
||
func simpsonRef(f func(float64) float64, phi float64) float64 {
|
||
const n = 200000
|
||
h := phi / float64(n)
|
||
s := f(0) + f(phi)
|
||
for i := 1; i < n; i++ {
|
||
if i%2 == 1 {
|
||
s += 4 * f(float64(i)*h)
|
||
} else {
|
||
s += 2 * f(float64(i)*h)
|
||
}
|
||
}
|
||
return s * h / 3
|
||
}
|
||
|
||
// TestEllipticK pins K against the closed form at m = 1/2, the AGM
|
||
// degenerate values, and Simpson.
|
||
func TestEllipticK(t *testing.T) {
|
||
k, err := EllipticK(mustFloats(t, []float64{0, 0.5, 0.9, -1.5}, 4))
|
||
if err != nil {
|
||
t.Fatalf("EllipticK: %v", err)
|
||
}
|
||
if math.Abs(k.FloatAt(0)-math.Pi/2) > 1e-15 {
|
||
t.Fatalf("K(0) = %.15f, want π/2", k.FloatAt(0))
|
||
}
|
||
// K(1/2) = Γ(1/4)² / (4√π).
|
||
want := math.Gamma(0.25) * math.Gamma(0.25) / (4 * math.Sqrt(math.Pi))
|
||
if math.Abs(k.FloatAt(1)-want) > 1e-14 {
|
||
t.Fatalf("K(0.5) = %.15f, want %.15f", k.FloatAt(1), want)
|
||
}
|
||
// m = 0.9 against Simpson on the integrand.
|
||
f := func(th float64) float64 { return 1 / math.Sqrt(1-0.9*math.Sin(th)*math.Sin(th)) }
|
||
if math.Abs(k.FloatAt(2)-simpsonRef(f, math.Pi/2)) > 1e-11 {
|
||
t.Fatalf("K(0.9) = %.14f, Simpson says %.14f", k.FloatAt(2), simpsonRef(f, math.Pi/2))
|
||
}
|
||
// Negative parameter against Simpson through its own integrand.
|
||
fn := func(th float64) float64 { return 1 / math.Sqrt(1+1.5*math.Sin(th)*math.Sin(th)) }
|
||
if math.Abs(k.FloatAt(3)-simpsonRef(fn, math.Pi/2)) > 1e-11 {
|
||
t.Fatalf("K(-1.5) = %.14f, Simpson says %.14f", k.FloatAt(3), simpsonRef(fn, math.Pi/2))
|
||
}
|
||
inf, err := EllipticK(mustFloats(t, []float64{1}, 1))
|
||
if err != nil || !math.IsInf(inf.FloatAt(0), 1) {
|
||
t.Fatalf("K(1) must be +Inf, got %v err %v", inf.FloatAt(0), err)
|
||
}
|
||
}
|
||
|
||
// TestEllipticE pins E against Simpson and the AGM-series degenerate
|
||
// values, plus the negative-parameter transform.
|
||
func TestEllipticE(t *testing.T) {
|
||
e, err := EllipticE(mustFloats(t, []float64{0, 0.5, 0.99, -2.0, 1.0}, 5))
|
||
if err != nil {
|
||
t.Fatalf("EllipticE: %v", err)
|
||
}
|
||
if math.Abs(e.FloatAt(0)-math.Pi/2) > 1e-15 {
|
||
t.Fatalf("E(0) = %.15f, want π/2", e.FloatAt(0))
|
||
}
|
||
for i, m := range []float64{0.5, 0.99} {
|
||
f := func(th float64) float64 { return math.Sqrt(1 - m*math.Sin(th)*math.Sin(th)) }
|
||
if math.Abs(e.FloatAt(1+i)-simpsonRef(f, math.Pi/2)) > 1e-11 {
|
||
t.Fatalf("E(%g) = %.14f, Simpson says %.14f", m, e.FloatAt(1+i), simpsonRef(f, math.Pi/2))
|
||
}
|
||
}
|
||
fn := func(th float64) float64 { return math.Sqrt(1 + 2.0*math.Sin(th)*math.Sin(th)) }
|
||
if math.Abs(e.FloatAt(3)-simpsonRef(fn, math.Pi/2)) > 1e-11 {
|
||
t.Fatalf("E(-2) = %.14f, Simpson says %.14f", e.FloatAt(3), simpsonRef(fn, math.Pi/2))
|
||
}
|
||
if math.Abs(e.FloatAt(4)-1) > 1e-15 {
|
||
t.Fatalf("E(1) = %.15f, want 1", e.FloatAt(4))
|
||
}
|
||
}
|
||
|
||
// TestEllipticPi pins Pi against Simpson, including the Pi(0, m) = K
|
||
// identity.
|
||
func TestEllipticPi(t *testing.T) {
|
||
n := mustFloats(t, []float64{0, 0.5, -1.0}, 3)
|
||
m := mustFloats(t, []float64{0.5, 0.8, 0.3}, 3)
|
||
pi, err := EllipticPi(n, m)
|
||
if err != nil {
|
||
t.Fatalf("EllipticPi: %v", err)
|
||
}
|
||
k05 := EllipticKScalar(0.5)
|
||
if math.Abs(pi.FloatAt(0)-k05) > 1e-13 {
|
||
t.Fatalf("Π(0, 0.5) = %.14f, K(0.5) = %.14f", pi.FloatAt(0), k05)
|
||
}
|
||
for i, pair := range [][2]float64{{0.5, 0.8}, {-1.0, 0.3}} {
|
||
nn, mm := pair[0], pair[1]
|
||
f := func(th float64) float64 {
|
||
sq := math.Sin(th)
|
||
return 1 / ((1 - nn*sq*sq) * math.Sqrt(1-mm*sq*sq))
|
||
}
|
||
if math.Abs(pi.FloatAt(1+i)-simpsonRef(f, math.Pi/2)) > 1e-11 {
|
||
t.Fatalf("Π(%g, %g) = %.14f, Simpson says %.14f", nn, mm, pi.FloatAt(1+i), simpsonRef(f, math.Pi/2))
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestJacobiIdentities pins sn, cn, dn against their defining
|
||
// identities, degenerate parameters, the derivative relation and the
|
||
// inversion that defines them (F(am(u)) = u).
|
||
func TestJacobiIdentities(t *testing.T) {
|
||
us := mustFloats(t, []float64{0.3, 1.2, 2.7, -0.9, 5.5}, 5)
|
||
for _, m := range []float64{0.0, 0.37, 0.96} {
|
||
sn, err := JacobiSN(us, m)
|
||
if err != nil {
|
||
t.Fatalf("JacobiSN: %v", err)
|
||
}
|
||
cn, err := JacobiCN(us, m)
|
||
if err != nil {
|
||
t.Fatalf("JacobiCN: %v", err)
|
||
}
|
||
dn, err := JacobiDN(us, m)
|
||
if err != nil {
|
||
t.Fatalf("JacobiDN: %v", err)
|
||
}
|
||
for i := range 5 {
|
||
u := us.FloatAt(i)
|
||
s, c, d := sn.FloatAt(i), cn.FloatAt(i), dn.FloatAt(i)
|
||
if math.Abs(s*s+c*c-1) > 1e-12 {
|
||
t.Fatalf("m=%g u=%g: sn²+cn² = %.14f", m, u, s*s+c*c)
|
||
}
|
||
if math.Abs(m*s*s+d*d-1) > 1e-12 {
|
||
t.Fatalf("m=%g u=%g: m·sn²+dn² = %.14f", m, u, m*s*s+d*d)
|
||
}
|
||
if u < 0 && math.Abs(s+sn0(us, m, -u)) > 1e-13 {
|
||
t.Fatalf("m=%g: sn is not odd", m)
|
||
}
|
||
if u < 0 && math.Abs(d-dn0(us, m, -u)) > 1e-13 {
|
||
t.Fatalf("m=%g: dn is not even", m)
|
||
}
|
||
// d/du sn = cn·dn by central differences.
|
||
const h = 1e-6
|
||
ups := mustFloats(t, []float64{u + h, u - h}, 2)
|
||
sp, _ := JacobiSN(ups, m)
|
||
der := (sp.FloatAt(0) - sp.FloatAt(1)) / (2 * h)
|
||
if math.Abs(der-c*d) > 1e-6 {
|
||
t.Fatalf("m=%g u=%g: sn' = %.8f, cn·dn = %.8f", m, u, der, c*d)
|
||
}
|
||
}
|
||
}
|
||
// Degenerate parameters: m = 0 circular, m = 1 hyperbolic.
|
||
us2 := mustFloats(t, []float64{0.4, 1.3}, 2)
|
||
sn, _ := JacobiSN(us2, 0)
|
||
cn, _ := JacobiCN(us2, 0)
|
||
if math.Abs(sn.FloatAt(0)-math.Sin(0.4)) > 1e-14 || math.Abs(cn.FloatAt(1)-math.Cos(1.3)) > 1e-14 {
|
||
t.Fatal("m = 0 must degenerate to sin/cos")
|
||
}
|
||
// The hyperbolic limit m tending to 1: sn approaches tanh from
|
||
// inside a boundary layer of width ~sqrt(1−m), hence the loose
|
||
// tolerance.
|
||
sn1, err := JacobiSN(us2, 1-1e-12)
|
||
if err != nil {
|
||
t.Fatalf("JacobiSN near m=1: %v", err)
|
||
}
|
||
if math.Abs(sn1.FloatAt(1)-math.Tanh(1.3)) > 1e-4 {
|
||
t.Fatalf("m tending to 1: sn(1.3) = %.12f, tanh = %.12f", sn1.FloatAt(1), math.Tanh(1.3))
|
||
}
|
||
if _, err := JacobiSN(us2, 1); err == nil {
|
||
t.Fatal("m = 1 must be rejected (K diverges there)")
|
||
}
|
||
// The defining inversion: F(asin(sn(u, m)), m) = u.
|
||
// Valid below K(m) ~ 1.995: past it the amplitude passes π/2 and
|
||
// asin folds the branch away.
|
||
for _, u := range []float64{0.7, 1.5} {
|
||
ua := mustFloats(t, []float64{u}, 1)
|
||
snu, _ := JacobiSN(ua, 0.6)
|
||
phi := math.Asin(snu.FloatAt(0))
|
||
if math.Abs(ellipticF(phi, 0.6)-u) > 1e-10 {
|
||
t.Fatalf("F(asin sn(%g)) = %.12f, want %g", u, ellipticF(phi, 0.6), u)
|
||
}
|
||
}
|
||
// Periodicity: sn(u + 4K) = sn(u).
|
||
k := EllipticKScalar(0.5)
|
||
uper := mustFloats(t, []float64{0.9 + 4*k}, 1)
|
||
uplain := mustFloats(t, []float64{0.9}, 1)
|
||
s1, _ := JacobiSN(uper, 0.5)
|
||
s2, _ := JacobiSN(uplain, 0.5)
|
||
if math.Abs(s1.FloatAt(0)-s2.FloatAt(0)) > 1e-10 {
|
||
t.Fatalf("periodicity broke: %.12f vs %.12f", s1.FloatAt(0), s2.FloatAt(0))
|
||
}
|
||
}
|
||
|
||
// sn0/dn0 are scalar convenience reads for the parity checks.
|
||
func sn0(us *Array, m, u float64) float64 {
|
||
one, _ := FromFloats([]float64{u}, 1)
|
||
s, err := JacobiSN(one, m)
|
||
if err != nil {
|
||
panic(err)
|
||
}
|
||
return s.FloatAt(0)
|
||
}
|
||
|
||
func dn0(us *Array, m, u float64) float64 {
|
||
one, _ := FromFloats([]float64{u}, 1)
|
||
d, err := JacobiDN(one, m)
|
||
if err != nil {
|
||
panic(err)
|
||
}
|
||
return d.FloatAt(0)
|
||
}
|
||
|
||
// TestHypergeometric2F1 pins 2F1 against its closed forms.
|
||
func TestHypergeometric2F1(t *testing.T) {
|
||
x := mustFloats(t, []float64{0, 0.3, 0.7, 0.97, -0.8, 1.0}, 6)
|
||
f, err := Hypergeometric2F1(1, 1, 2, x)
|
||
if err != nil {
|
||
t.Fatalf("Hypergeometric2F1: %v", err)
|
||
}
|
||
for i, v := range x.RawFloats() {
|
||
want := -math.Log(1-v) / v
|
||
if v == 0 {
|
||
want = 1
|
||
}
|
||
if math.Abs(f.FloatAt(i)-want) > 1e-12 {
|
||
t.Fatalf("2F1(1,1;2;%g) = %.14f, want %.14f", v, f.FloatAt(i), want)
|
||
}
|
||
}
|
||
// 2F1(a, b; b; x) = (1−x)^{−a}.
|
||
g, err := Hypergeometric2F1(0.7, 3.5, 3.5, x)
|
||
if err != nil {
|
||
t.Fatalf("Hypergeometric2F1: %v", err)
|
||
}
|
||
for i, v := range x.RawFloats() {
|
||
if v == 1 {
|
||
continue
|
||
}
|
||
want := math.Pow(1-v, -0.7)
|
||
if math.Abs(g.FloatAt(i)-want) > 1e-12 {
|
||
t.Fatalf("2F1(a,b;b;%g) = %.14f, want %.14f", v, g.FloatAt(i), want)
|
||
}
|
||
}
|
||
// Terminating polynomial: 2F1(-3, 2; 1.5; x), built term by term below.
|
||
h, err := Hypergeometric2F1(-3, 2, 1.5, mustFloats(t, []float64{0.6}, 1))
|
||
if err != nil {
|
||
t.Fatalf("Hypergeometric2F1: %v", err)
|
||
}
|
||
// k=0: 1; k=1: (−3·2/1.5)·0.6 = −2.4; k=2: (−3·−2·2·3/(1.5·2.5·2))·0.36;
|
||
// k=3: (−3·−2·−1·2·3·4/(1.5·2.5·3.5·6))·0.216.
|
||
poly := 1.0
|
||
poly += (-3.0 * 2.0 / 1.5) * 0.6
|
||
poly += (-3.0 * -2.0 * 2.0 * 3.0 / (1.5 * 2.5 * 2.0)) * 0.36
|
||
poly += (-3.0 * -2.0 * -1.0 * 2.0 * 3.0 * 4.0 / (1.5 * 2.5 * 3.5 * 6.0)) * 0.216
|
||
if math.Abs(h.FloatAt(0)-poly) > 1e-14 {
|
||
t.Fatalf("2F1(-3,2;1.5;0.6) = %.14f, want %.14f", h.FloatAt(0), poly)
|
||
}
|
||
// c a non-positive integer is an error.
|
||
if _, err := Hypergeometric2F1(1, 1, -1, x); err == nil {
|
||
t.Fatal("c = −1 accepted")
|
||
}
|
||
}
|
||
|
||
// TestJacobiCDScalar pins the AGM evaluation against mpmath 3.16
|
||
// reference values (ellipfun cn/dn at 30 digits) over the parameter
|
||
// range the Cauer filter design walks: m near both ends and u across
|
||
// several periods.
|
||
func TestJacobiCDScalar(t *testing.T) {
|
||
cases := []struct {
|
||
u, m, want float64
|
||
}{
|
||
{0.3, 0.5, 9.77250336444249856e-01},
|
||
{1.0, 0.5, 7.24009721659370498e-01},
|
||
{1.7, 0.9, 7.12053939073435838e-01},
|
||
{2.5, 0.1, -7.69223750223807068e-01},
|
||
{3.2, 0.7, -8.39726459768202815e-01},
|
||
{5.0, 0.99, -8.64163207523496957e-01},
|
||
{7.5, 0.6, 9.81944916682632396e-01},
|
||
{12.0, 0.25, 1.98099136005263327e-01},
|
||
{0.0, 0.4, 1},
|
||
{1.0, 0.0, 5.40302305868139765e-01},
|
||
}
|
||
for _, c := range cases {
|
||
got := JacobiCDScalar(c.u, c.m)
|
||
if math.Abs(got-c.want) > 5e-16*math.Max(1, math.Abs(c.want)) {
|
||
t.Fatalf("cd(%.3g, %.3g) = %.16g, want %.16g", c.u, c.m, got, c.want)
|
||
}
|
||
}
|
||
if v := JacobiCDScalar(1, math.NaN()); !math.IsNaN(v) {
|
||
t.Fatal("NaN parameter accepted")
|
||
}
|
||
if v := JacobiCDScalar(1, 1); !math.IsNaN(v) {
|
||
t.Fatal("m = 1 accepted")
|
||
}
|
||
}
|