// 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)
		}
	}
}