tensor/internal/core/besseljreal_test.go

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

package core

import (
	"math"
	"testing"
)

// jHalf evaluates the closed forms of the half-integer orders
// J₁/₂, J₃/₂ and J₅/₂, the exact referents no other Bessel test here
// enjoys.
func jHalf(nu, x float64) float64 {
	s := math.Sqrt(2 / (math.Pi * x))
	sin, cos := math.Sincos(x)
	switch nu {
	case 0.5:
		return s * sin
	case 1.5:
		return s * (sin/x - cos)
	case 2.5:
		return s * ((3/(x*x)-1)*sin - 3*cos/x)
	}
	panic("jHalf: unsupported order")
}

// climbHalf climbs the closed forms from orders 1/2 and 3/2 up to
// 0.5+m by the exact three-term recurrence, the referent for orders
// the closed forms themselves do not cover.
func climbHalf(m int, x float64) float64 {
	jm := jHalf(0.5, x)
	j := jHalf(1.5, x)
	for k := 1; k < m; k++ {
		jm, j = j, 2*(0.5+float64(k))/x*j-jm
	}
	return j
}

// TestBesselJRealOrderHalfIntegers holds every regime against the
// closed forms: the series below the crossover, the climb seeded from
// the expansion above it, and a phase reduction at an argument whose
// plain float64 phase would already be losing digits.
func TestBesselJRealOrderHalfIntegers(t *testing.T) {
	// The sweep starts at 0.5: below that the closed forms themselves
	// cancel catastrophically and stop being referents. The small-x
	// behaviour is pinned separately, against the series' own leading
	// terms.
	xs := []float64{0.5, 1, 2, 7, 12.1, 14.9, 15, 20, 40, 100, 1000}
	for _, nu := range []float64{0.5, 1.5, 2.5} {
		for _, x := range xs {
			got, err := BesselJRealOrder(nu, x)
			if err != nil {
				t.Fatalf("BesselJRealOrder(%g, %g): %v", nu, x, err)
			}
			want := jHalf(nu, x)
			if d := math.Abs(got-want) / math.Abs(want); d > 5e-12 {
				t.Fatalf("BesselJRealOrder(%g, %g) = %.17g, want %.17g (relative %.3g)", nu, x, got, want, d)
			}
		}
	}
	// At x = 0.05 the true J₅/₂ agrees with the series' first three
	// terms to five parts in 1e10, the third term being the first one
	// the referent omits.
	const x, nu = 0.05, 2.5
	got, err := BesselJRealOrder(nu, x)
	if err != nil {
		t.Fatal(err)
	}
	half := 0.5 * x
	t0 := math.Exp(nu*math.Log(half) - lnGammaReal(nu+1))
	q := half * half
	referent := t0 * (1 - q/(nu+1)*(1-q/(2*(nu+2))))
	if d := math.Abs(got-referent) / t0; d > 1e-10 {
		t.Fatalf("BesselJRealOrder(%g, %g) = %.17g, want %.17g (relative %.3g)", nu, x, got, referent, d)
	}
}

// TestBesselJRealOrderMiller pins the fractional Miller walk: an order
// past the argument at an argument past the crossover, against the
// exact closed forms climbed up by the recurrence.
func TestBesselJRealOrderMiller(t *testing.T) {
	const x = 15
	got, err := BesselJRealOrder(20.5, x)
	if err != nil {
		t.Fatalf("BesselJRealOrder: %v", err)
	}
	want := climbHalf(20, x)
	if d := math.Abs(got-want) / math.Abs(want); d > 1e-10 {
		t.Fatalf("BesselJRealOrder(20.5, 15) = %.17g, want %.17g (relative %.3g)", got, want, d)
	}
}

// TestBesselJRealOrderRecurrence checks the three-term recurrence
// across the map of regimes, the generic verifier no closed form can
// replace: every pair of orders the test touches lives on one walk.
func TestBesselJRealOrderRecurrence(t *testing.T) {
	cases := [][2]float64{
		{1.3, 15}, {2.3, 15}, {7.7, 40}, {12.3, 100},
		{25.5, 15}, {60.5, 40}, {1.7, 1000},
	}
	for _, c := range cases {
		nu, x := c[0], c[1]
		jm, err := BesselJRealOrder(nu-1, x)
		if err != nil {
			t.Fatalf("BesselJRealOrder(%g, %g): %v", nu-1, x, err)
		}
		j, err := BesselJRealOrder(nu, x)
		if err != nil {
			t.Fatalf("BesselJRealOrder(%g, %g): %v", nu, x, err)
		}
		jp, err := BesselJRealOrder(nu+1, x)
		if err != nil {
			t.Fatalf("BesselJRealOrder(%g, %g): %v", nu+1, x, err)
		}
		scale := math.Max(math.Abs(jm), math.Max(math.Abs(jp), math.Abs(2*nu/x*j)))
		residual := math.Abs(jm + jp - 2*nu/x*j)
		// Near the crossover the expansion's truncation caps the seeds
		// at about eight digits; further out it truncates past the
		// rounding floor and the residual follows it down.
		tol := 1e-11
		if x < 20 {
			tol = 5e-8
		}
		if residual > tol*scale {
			t.Fatalf("the recurrence residual at (ν = %g, x = %g) is %.3g against scale %.3g", nu, x, residual, scale)
		}
	}
}

// TestBesselJRealOrderIntegerDelegate pins the integer delegation: an
// exact integer order returns the integer algorithm's own bits, and an
// order a whisper away from it lands within a whisper of the same
// value through the general route.
func TestBesselJRealOrderIntegerDelegate(t *testing.T) {
	for _, x := range []float64{2, 20} {
		for _, n := range []int{0, 3, 17} {
			got, err := BesselJRealOrder(float64(n), x)
			if err != nil {
				t.Fatalf("BesselJRealOrder(%d, %g): %v", n, x, err)
			}
			if want := BesselJ(n, x); got != want {
				t.Fatalf("BesselJRealOrder(%d, %g) = %.17g, want the integer %.17g", n, x, got, want)
			}
		}
	}
	near, err := BesselJRealOrder(3+1e-12, 20)
	if err != nil {
		t.Fatal(err)
	}
	if d := math.Abs(near - BesselJ(3, 20)); d > 1e-11 {
		t.Fatalf("an order 1e-12 off the integer moved J by %g", d)
	}
}

// TestBesselJRealOrderSeriesVsClimb crosses the two routes against
// each other just above the crossover, where both still carry roughly
// nine digits: the series pushed past its regime and the climb from
// the expansion must agree to that shared quality.
func TestBesselJRealOrderSeriesVsClimb(t *testing.T) {
	const nu, x = 2.3, 12.1
	got, err := BesselJRealOrder(nu, x)
	if err != nil {
		t.Fatal(err)
	}
	want := besselJSeriesReal(nu, x)
	if d := math.Abs(got-want) / math.Abs(want); d > 1e-8 {
		t.Fatalf("the climb gives %.17g against the series %.17g (relative %.3g)", got, want, d)
	}
}

// TestBesselJRealOrderErrors pins the domain: a negative or NaN order
// and a non-positive or NaN argument are errors naming themselves.
func TestBesselJRealOrderErrors(t *testing.T) {
	for _, c := range [][2]float64{{-1, 2}, {math.NaN(), 2}, {1, 0}, {1, -2}, {math.NaN(), math.NaN()}} {
		if _, err := BesselJRealOrder(c[0], c[1]); err == nil {
			t.Fatalf("BesselJRealOrder(%g, %g): want an error", c[0], c[1])
		}
	}
}