135 lines
4.4 KiB
Go
135 lines
4.4 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
package core
|
||
|
||
import (
|
||
"math"
|
||
"testing"
|
||
)
|
||
|
||
// The third-kind integral and the Carlson forms behind it. Reference
|
||
// values are mpmath at 30 to 35 digits (elliprf, elliprj, elliprc, and
|
||
// Pi through the identity below, cross-checked against mpmath's
|
||
// quadrature); the library must reach double precision at the singular
|
||
// ends of the parameter square, where the product rule it replaced
|
||
// lost seven digits and worse.
|
||
|
||
// TestEllipticPiAccuracy pins Π against mpmath. The covered corners
|
||
// are n and m approaching 1 together (where Π reaches 1e6), n
|
||
// approaching 1, m approaching 1, negative parameters, and the
|
||
// degenerate m = 0.
|
||
func TestEllipticPiAccuracy(t *testing.T) {
|
||
cases := []struct {
|
||
n, m, want float64
|
||
}{
|
||
{0, 0.5, 1.8540746773013719184},
|
||
{0.5, 0.8, 3.4166601403243870137},
|
||
{-1, 0.3, 1.1936018953043909136},
|
||
{0.9, 0.9, 11.047747327040735532},
|
||
{0.99, 0.999, 193.39638212991131552},
|
||
{0.999, 0.5, 69.434652042115458236},
|
||
{0.5, 0, 2.2214414690791831235},
|
||
{-0.5, 0.5, 1.4878469926687983853},
|
||
{0.9, -0.5, 4.2505802986876906623},
|
||
{-2, -3, 0.68305896638359993325},
|
||
{0.999999, 0.999999, 1000003.8969974163894},
|
||
{0, 0.999999999, 11.747927296421043878},
|
||
{0.3, 0.999999999, 16.301444430032469214},
|
||
{0.99, 0.999999, 531.60692473385781712},
|
||
}
|
||
for _, tc := range cases {
|
||
got := ellipticPiScalar(tc.n, tc.m)
|
||
if rel := math.Abs(got-tc.want) / math.Abs(tc.want); rel > 1e-14 {
|
||
t.Errorf("Π(%g, %g) = %.17g, mpmath says %.17g (relative %.2g)",
|
||
tc.n, tc.m, got, tc.want, rel)
|
||
}
|
||
}
|
||
// The identities: Π(0, m) = K(m), and the divergence contract.
|
||
n := mustFloats(t, []float64{0}, 1)
|
||
m := mustFloats(t, []float64{0.999999}, 1)
|
||
pi, err := EllipticPi(n, m)
|
||
if err != nil {
|
||
t.Fatalf("EllipticPi: %v", err)
|
||
}
|
||
if rel := math.Abs(pi.FloatAt(0)-EllipticKScalar(0.999999)) / EllipticKScalar(0.999999); rel > 1e-15 {
|
||
t.Errorf("Π(0, 0.999999) = %.17g, K = %.17g", pi.FloatAt(0), EllipticKScalar(0.999999))
|
||
}
|
||
edgeN := mustFloats(t, []float64{1, 1.5}, 2)
|
||
edgeM := mustFloats(t, []float64{0.5, 0.5}, 2)
|
||
edge, err := EllipticPi(edgeN, edgeM)
|
||
if err != nil {
|
||
t.Fatalf("EllipticPi: %v", err)
|
||
}
|
||
if !math.IsInf(edge.FloatAt(0), 1) {
|
||
t.Errorf("Π(1, 0.5) = %v, want +Inf", edge.FloatAt(0))
|
||
}
|
||
if !math.IsNaN(edge.FloatAt(1)) {
|
||
t.Errorf("Π(1.5, 0.5) = %v, want NaN", edge.FloatAt(1))
|
||
}
|
||
}
|
||
|
||
// TestCarlsonRJAccuracy pins R_J, including the arguments the third-kind
|
||
// identity hands it: a zero first argument and a small fourth one.
|
||
func TestCarlsonRJAccuracy(t *testing.T) {
|
||
cases := []struct {
|
||
x, y, z, p, want float64
|
||
}{
|
||
{1, 1, 1, 1, 1.0},
|
||
{0, 0.5, 1, 0.5, 5.0832785087638745196},
|
||
{0, 1e-12, 1, 1e-06, 22802707.378631573755},
|
||
{1e-20, 1, 2, 3, 0.77688623771511264203},
|
||
{0, 1, 1, 1e-16, 471238893.32608005743},
|
||
{0.1, 0.2, 0.3, 0.9, 4.2398211507913140837},
|
||
}
|
||
for _, tc := range cases {
|
||
got := carlsonRJ(tc.x, tc.y, tc.z, tc.p)
|
||
if rel := math.Abs(got-tc.want) / math.Abs(tc.want); rel > 1e-13 {
|
||
t.Errorf("R_J(%g, %g, %g, %g) = %.17g, mpmath says %.17g (relative %.2g)",
|
||
tc.x, tc.y, tc.z, tc.p, got, tc.want, rel)
|
||
}
|
||
}
|
||
// A negative fourth argument needs the principal value, which this
|
||
// implementation refuses rather than approximates.
|
||
if got := carlsonRJ(1, 1, 1, -1); !math.IsNaN(got) {
|
||
t.Errorf("R_J(1, 1, 1, −1) = %v, want NaN", got)
|
||
}
|
||
}
|
||
|
||
// TestCarlsonRCAccuracy pins R_C, both argument orders and the Cauchy
|
||
// principal value for a negative second argument.
|
||
func TestCarlsonRCAccuracy(t *testing.T) {
|
||
cases := []struct {
|
||
x, y, want float64
|
||
}{
|
||
{1, 1, 1.0},
|
||
{0, 1, 1.5707963267948966192},
|
||
{1, 2, 0.78539816339744830962},
|
||
{2, 1, 0.881373587019543025},
|
||
{1.5, 0.5, 1.14621583478058884},
|
||
{0.5, 1.5, 0.955316618124509278},
|
||
{0.1, 10, 0.467396548061524389},
|
||
{10, 0.1, 0.951308668352240085},
|
||
{1, 0.25, 1.5206919926018927},
|
||
{0.25, 1, 1.20919957615614523},
|
||
{1e-30, 1, 1.5707963267948956192},
|
||
{1, 1e-12, 14.508657738531223752},
|
||
{1, -0.5, 0.93588131010357011049},
|
||
{1, -3, 0.27465307216702742285},
|
||
{0, -1, 0},
|
||
}
|
||
for _, tc := range cases {
|
||
got := carlsonRC(tc.x, tc.y)
|
||
if tc.want == 0 {
|
||
if got != 0 {
|
||
t.Errorf("R_C(%g, %g) = %v, want 0", tc.x, tc.y, got)
|
||
}
|
||
continue
|
||
}
|
||
if rel := math.Abs(got-tc.want) / math.Abs(tc.want); rel > 1e-14 {
|
||
t.Errorf("R_C(%g, %g) = %.17g, mpmath says %.17g (relative %.2g)",
|
||
tc.x, tc.y, got, tc.want, rel)
|
||
}
|
||
}
|
||
}
|