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

package core

import "math"

// Ordinary Bessel functions of the first and second kind at one real
// point, the scalar companions of the modified pair in besselmod.go.
// J grows out of a convergent power series below the crossover and a
// downward Miller recurrence above it, the direction that amplifies
// no rounding once the order passes the argument. Y climbs the
// upward recurrence from the seeds Y₀ and Y₁, the stable direction
// for the second kind, with the seeds themselves from the Frobenius
// series below the crossover and from the asymptotic expansion above
// it, where the ascending series would start paying for cancellation.

// besselCrossover splits the power-series and recurrence regimes. At
// the crossover both sides still carry twelve significant digits, so
// the exact split point is a matter of taste rather than accuracy.
const besselCrossover = 15.0

// BesselJ returns the Bessel function of the first kind of integer
// order n at the real point x, Jₙ(x). Arguments with |x| below the
// crossover are served by the convergent power series
// Σ (−1)^k (x/2)^{2k+n}/(k!·Γ(k+n+1)); above it the recurrence runs in
// its stable direction, orders at or below the argument climbing
// upward from the large-argument asymptotic J₀ and J₁ at O(n) cost and
// higher orders running the downward Miller walk anchored on the same
// asymptotic J₀, whose start sits above the turning point at order
// n + |x|. A negative argument follows the parity law
// Jₙ(−x) = (−1)ⁿ Jₙ(x) and a negative order the law
// J₋ₙ(x) = (−1)ⁿ Jₙ(x); Jₙ is finite for every real x, so nothing here
// can fail.
func BesselJ(n int, x float64) float64 {
	if x == 0 {
		if n == 0 {
			return 1
		}
		return 0
	}
	// Fold both parity laws into one sign: (−1)ⁿ only cares about the
	// order's parity, which |n| preserves.
	order, ax, flip := n, math.Abs(x), 1.0
	if order < 0 {
		if order%2 != 0 {
			flip = -1
		}
		order = -order
	}
	if x < 0 && order%2 != 0 {
		flip = -flip
	}
	if ax < besselCrossover {
		return flip * besselJSeries(order, ax)
	}
	if float64(order) <= ax {
		return flip * besselJUpward(order, ax)
	}
	return flip * besselJMiller(order, ax)
}

// BesselY returns the Bessel function of the second kind of integer
// order n at the real point x, Yₙ(x), defined for x > 0; Y diverges
// at the origin and a non-positive argument is an error, not a NaN.
// Y₀ and Y₁ are seeded below the crossover by the Frobenius series
// (Abramowitz & Stegun 9.1.11 with ψ(k+1) = H_k − γ written out
// through the harmonic numbers H_k) and above it by the
// large-argument asymptotic expansion; higher orders climb the upward
// recurrence Yₙ₊₁ = 2n/x·Yₙ − Yₙ₋₁, stable for the second kind
// because the parasitic J component the seeds carry decays relative
// to Y at every step. A negative order follows the parity law
// Y₋ₙ(x) = (−1)ⁿ Yₙ(x).
//
// Errors: x ≤ 0 or NaN.
func BesselY(n int, x float64) (float64, error) {
	if math.IsNaN(x) || x <= 0 {
		return 0, errf("BesselY: the argument must be positive, got %g", x)
	}
	flip, order := 1.0, n
	if order < 0 {
		if order%2 != 0 {
			flip = -1
		}
		order = -order
	}
	var v float64
	switch {
	case order == 0:
		v = besselY0(x)
	case order == 1:
		v = besselY1(x)
	default:
		ym, y := besselY0(x), besselY1(x)
		for k := 1; k < order; k++ {
			ym, y = y, 2*float64(k)/x*y-ym
		}
		v = y
	}
	return flip * v, nil
}

// besselJSeries evaluates Jₙ(x) by the convergent power series,
// stepping the summand along tₖ = tₖ₋₁·(−(x/2)²)/(k(k+n)). The first
// term (x/2)ⁿ/n! goes through logarithms so a large order never
// overflows the factorial on the way to a small answer.
func besselJSeries(n int, x float64) float64 {
	half := 0.5 * x
	// The k = 0 term (x/2)ⁿ/n! is positive; the (−1)^k alternation
	// enters through the recursion step below.
	t := math.Exp(float64(n)*math.Log(half) - LnFactorial(n))
	sum := t
	q := half * half
	for k := 1; k <= 400; k++ {
		t *= -q / (float64(k) * float64(k+n))
		sum += t
		if math.Abs(t) <= 1e-17*math.Abs(sum) {
			break
		}
	}
	return sum
}

// besselJUpward climbs J₀ and J₁ from the large-argument asymptotic to
// order n by the upward recurrence Jₖ₊₁ = (2k/x)·Jₖ − Jₖ₋₁. The
// asymptotic seeds carry the absolute scale, so nothing needs
// renormalising, and the climb is the stable direction while the order
// stays at or below the argument: there the parasitic Y component the
// seeds carry stays bounded relative to J, while the downward walk
// would start above the turning point at order x and pay O(x) steps
// for an answer of order n. The caller guarantees x ≥ besselCrossover
// and 0 ≤ n ≤ x.
func besselJUpward(n int, x float64) float64 {
	j0, _ := besselAsymptotic(0, x)
	if n == 0 {
		return j0
	}
	j1, _ := besselAsymptotic(1, x)
	for k := 1; k < n; k++ {
		j0, j1 = j1, 2*float64(k)/x*j1-j0
	}
	return j1
}

// besselJMiller evaluates Jₙ(x) for |x| above the crossover by the
// downward Miller recurrence, the same stable scheme
// besselInMiller uses for the modified kind: start well above n with
// an arbitrary scale, recurse down through J_{k−1} = (2k/x)·J_k −
// J_{k+1}, then renormalise the arbitrary seed scale against an
// independent J₀, here the asymptotic value the way besselInMiller
// leans on its quadrature I₀. The caller guarantees x ≠ 0.
func besselJMiller(n int, x float64) float64 {
	start := n + int(x) + 40
	jp, j := 0.0, 1.0 // J_{k+1}, J_k, seeded at k = start
	jn := 0.0
	for k := start; k >= 1; k-- {
		if k == n {
			jn = j
		}
		jp, j = j, 2*float64(k)/x*j-jp
		if aj := math.Abs(j); aj > 1e200 {
			// The unscaled seed grows like k!/x^k on the way down and
			// would overflow for high orders whose true value is
			// representable; the renormalisation cancels any common
			// factor, so rescaling the running pair (and the captured
			// order-n value) is exact up to rounding.
			jp /= aj
			j /= aj
			jn /= aj
		}
	}
	if n == 0 {
		jn = j
	}
	anchor, _ := besselAsymptotic(0, x)
	return anchor * jn / j
}

// besselY0 evaluates Y₀(x) for x > 0. Below the crossover the
// Frobenius series in the form the digamma reduction gives,
//
//	Y₀ = (2/π)·[(ln(x/2) + γ)·J₀(x) + Σ (−1)^{k+1} H_k (x/2)^{2k}/(k!)²],
//
// above it the asymptotic expansion, whose terms still fall fast
// enough at the crossover to keep twelve significant digits.
func besselY0(x float64) float64 {
	if x >= besselCrossover {
		_, y := besselAsymptotic(0, x)
		return y
	}
	half := 0.5 * x
	q := half * half
	u := q   // u_k = (x/2)^{2k}/(k!)², starting at k = 1
	h := 1.0 // H_k
	sum := u // the k = 1 term carries the + sign
	for k := 2; k <= 200; k++ {
		u *= q / float64(k*k)
		h += 1 / float64(k)
		term := h * u
		if k%2 == 1 {
			sum += term
		} else {
			sum -= term
		}
		if math.Abs(term) <= 1e-17*math.Abs(sum) {
			break
		}
	}
	return (2 / math.Pi) * ((math.Log(half)+eulerGamma)*besselJSeries(0, x) + sum)
}

// besselY1 evaluates Y₁(x) for x > 0, the order-one twin of
// besselY0's series branch:
//
//	Y₁ = (2/π)(ln(x/2) + γ)·J₁(x) − 2/(πx)
//	     − (1/π)·Σ (−1)^k (H_k + H_{k+1})·(x/2)^{2k+1}/(k!(k+1)!),
//
// where the 2/(πx) term is the one-entry finite sum of A&S 9.1.11
// and the ascending series converges for every x, paying only the
// cancellation that caps its usable range at the crossover.
func besselY1(x float64) float64 {
	if x >= besselCrossover {
		_, y := besselAsymptotic(1, x)
		return y
	}
	half := 0.5 * x
	q := half * half
	u := half // u_k = (x/2)^{2k+1}/(k!(k+1)!), starting at k = 0
	h := 1.0  // H_{k+1}, starting at H_1 = 1 (H_0 = 0)
	sum := u  // the k = 0 term: (H_0 + H_1)·u₀ with the + sign
	for k := 1; k <= 200; k++ {
		hk := h // H_k, before the update below
		u *= q / float64(k*(k+1))
		h += 1 / float64(k+1)
		term := (hk + h) * u
		if k%2 == 1 {
			sum -= term
		} else {
			sum += term
		}
		if math.Abs(term) <= 1e-17*math.Abs(sum) {
			break
		}
	}
	return (2/math.Pi)*((math.Log(half)+eulerGamma)*besselJSeries(1, x)) -
		2/(math.Pi*x) - sum/math.Pi
}

// besselAsymptotic evaluates the pair J_ν, Y_ν above the crossover
// from the large-argument expansion (Abramowitz & Stegun 9.2.5 through
// 9.2.6):
//
//	J_ν ~ sqrt(2/πx)·[cos ω·Σ(−1)^k a_{2k} − sin ω·Σ(−1)^k a_{2k+1}]
//	Y_ν ~ sqrt(2/πx)·[sin ω·Σ(−1)^k a_{2k} + cos ω·Σ(−1)^k a_{2k+1}],
//
// with ω = x − νπ/2 − π/4 and aₘ = ∏(μ − (2s−1)²)/(m!·(8x)^m),
// μ = 4ν². The order is a float64: the integer callers pass exact
// integers, whose float64 arithmetic reproduces the int path bit for
// bit, and the real-order BesselJRealOrder passes the fractional seed
// orders. Both kinds share the one coefficient sweep, and the sum
// stops at the optimal truncation where the terms turn around and
// start growing again, which at the crossover still leaves about
// twelve significant digits.
func besselAsymptotic(nu float64, x float64) (j, y float64) {
	mu := 4 * nu * nu
	omega := besselPhase(nu, x)
	eighth := 8 * x
	a, prev := 1.0, math.Inf(1)
	even, odd := 1.0, 0.0 // Σ(−1)^k a_{2k}, Σ(−1)^k a_{2k+1}
	for m := 1; m <= 200; m++ {
		a *= (mu - float64((2*m-1)*(2*m-1))) / (float64(m) * eighth)
		// m = 2k and m = 2k+1 share the integer k = m/2 and the sign
		// (−1)^k of their sum's k-th term.
		if (m/2)%2 == 0 {
			if m%2 == 0 {
				even += a
			} else {
				odd += a
			}
		} else {
			if m%2 == 0 {
				even -= a
			} else {
				odd -= a
			}
		}
		if math.Abs(a) > prev {
			break
		}
		prev = math.Abs(a)
	}
	factor := math.Sqrt(2 / (math.Pi * x))
	sin, cos := math.Sin(omega), math.Cos(omega)
	j = factor * (cos*even - sin*odd)
	y = factor * (sin*even + cos*odd)
	return j, y
}

// 2π split into two float64 parts, the sum of which is 2π to about
// 1e-32: the phase reduction below needs the low part to hold the
// fraction of a large argument.
const (
	twoPiHi = 6.283185307179586      // fl(2π)
	twoPiLo = 2.4492935982947064e-16 // 2π − fl(2π)
)

// besselPhase returns ω = x − νπ/2 − π/4 reduced modulo 2π for the
// large-argument expansion. The reduction is what keeps the phase of a
// large argument meaningful: the plain float64 difference rounds the
// fraction of x away, at x = 1e9 to about 1e-7 absolute, which the
// amplitude √(2/πx) turns into a relative error of 1e-8, and at
// x = 1e12 into 1e-4. The remainder is taken with the two-part 2π:
// x − hi is exact by Sterbenz's lemma, math.FMA gives the exact
// residual of n·2π, and the rest is a handful of flops on quantities
// below π, so the phase keeps its own last ulp while |ω| < 2^53·2π
// (beyond that the float64 grid of x is coarser than a radian and the
// plain difference is as good as anything). The caller passes |x|,
// which is at least the crossover.
func besselPhase(nu float64, x float64) float64 {
	nuPi2 := nu * math.Pi / 2
	quarterPi := math.Pi / 4
	omega := x - nuPi2 - quarterPi
	n := math.Round(omega / twoPiHi)
	if math.Abs(n) >= 1<<53 {
		return omega
	}
	hi := n * twoPiHi
	lo := math.FMA(n, twoPiHi, -hi)
	return x - hi - lo - n*twoPiLo - nuPi2 - quarterPi
}

// BesselJRealOrder returns the Bessel function of the first kind of
// real order ν at the real point x, J_ν(x). The order must be ≥ 0 and
// the argument positive; J diverges at the origin for ν > 0 and a
// non-positive argument or a NaN order is an error, not a NaN.
//
// The regimes follow the integer BesselJ's, with the fractional part
// of the order taking the role of the seeds: below the crossover the
// ascending Frobenius series carries the answer, above it the order
// resolves into its integer and fractional parts, the fractional pair
// (ν₀, ν₀+1) is seeded from the large-argument expansion and the
// three-term recurrence runs in its stable direction, upward while the
// order stays at or below the argument and by the downward Miller walk
// anchored on the ν₀ seed when the order passes it. Orders within 1e-8
// of a non-negative integer are served by the integer algorithm
// itself, which is the continuous limit there: the general route's
// normalisation would cancel against sin(πν) and lose every digit.
func BesselJRealOrder(nu, x float64) (float64, error) {
	if math.IsNaN(nu) || math.IsNaN(x) {
		return 0, errf("BesselJRealOrder: the order and the argument must be finite, got %g and %g", nu, x)
	}
	if nu < 0 {
		return 0, errf("BesselJRealOrder: the order must be zero or greater, got %g", nu)
	}
	if x <= 0 {
		return 0, errf("BesselJRealOrder: the argument must be positive, got %g", x)
	}
	if n := math.Round(nu); math.Abs(nu-n) <= 1e-8 {
		return BesselJ(int(n), x), nil
	}
	if x < besselRealCrossover {
		return besselJSeriesReal(nu, x), nil
	}
	whole := math.Floor(nu)
	frac := nu - whole
	if nu <= x {
		// The climb: seed the fractional pair from the expansion and
		// run the upward recurrence, the stable direction while the
		// order stays at or below the argument. whole = 0 means the
		// answer is the first seed itself.
		jm, _ := besselAsymptotic(frac, x)
		if whole == 0 {
			return jm, nil
		}
		j, _ := besselAsymptotic(frac+1, x)
		for k := 1; k < int(whole); k++ {
			jm, j = j, 2*(frac+float64(k))/x*j-jm
		}
		return j, nil
	}
	// The downward Miller walk through one fractional residue class,
	// anchored on the fractional seed the same way besselJMiller
	// anchors its integer walk on J₀. The walk counts integer steps
	// from the fractional floor instead of testing k against nu and
	// frac: start descends in floats, and a capture missed by one ulp
	// of the running k would answer from an unseeded slot.
	steps := int(math.Floor(x)) + int(whole) + 40
	jp, j := 0.0, 1.0 // J_{k+1}, J_k, seeded at k = frac + steps
	uNu := 0.0
	uFrac := 0.0
	for i := steps; i >= 0; i-- {
		k := frac + float64(i)
		if i == int(whole) {
			uNu = j
		}
		if i == 0 {
			uFrac = j
		}
		jp, j = j, 2*k/x*j-jp
		if aj := math.Abs(j); aj > 1e200 {
			jp /= aj
			j /= aj
			uNu /= aj
			uFrac /= aj
		}
	}
	anchor, _ := besselAsymptotic(frac, x)
	return anchor * uNu / uFrac, nil
}

// lnGammaReal returns lnΓ(z) for z > 0 at one point, the scalar
// companion the real-order series needs; the sign of Γ is positive on
// that domain, so only the logarithm comes back.
func lnGammaReal(z float64) float64 {
	lg, _ := math.Lgamma(z)
	return lg
}

// besselRealCrossover splits the real-order series and recurrence
// regimes, lower than the integer one because the two sides' quality
// decides differently here: the series pays the same cancellation
// (about x·ln10/2 digits at x = 12, still leaving better than nine),
// while the expansion the seeds come from truncates to full precision
// for the half-integer orders and to ten-plus digits for the small
// fractional seeds the climb and the Miller walk lean on.
const besselRealCrossover = 12.0

// besselJSeriesReal evaluates J_ν(x) by the convergent Frobenius
// series for real ν, stepping the summand along
// tₖ = tₖ₋₁·(−(x/2)²)/(k(k+ν)). It is the real-order companion of
// besselJSeries, kept separate so the integer series keeps its
// exact LnFactorial opening and its recorded bits.
func besselJSeriesReal(nu, x float64) float64 {
	half := 0.5 * x
	// The opening term never overflows below the crossover: with
	// half < 7.5 the exponent ν·log(half) − lnΓ(ν+1) turns downward
	// past ν ≈ 20 and stays negative.
	t := math.Exp(nu*math.Log(half) - lnGammaReal(nu+1))
	sum := t
	q := half * half
	for k := 1; k <= 400; k++ {
		t *= -q / (float64(k) * (float64(k) + nu))
		sum += t
		if math.Abs(t) <= 1e-17*math.Abs(sum) {
			break
		}
	}
	return sum
}