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

package core

import "math"

// Modified Bessel functions of the first and second kind, the workhorses
// of heat conduction, waveguides and filtered noise. I₀ and I₁ use the
// polynomial fits below 3.75 (Abramowitz & Stegun 9.8.1 and 9.8.2) and
// their integral representations beyond, K₀ and K₁ integrate their
// kernels directly, and integer orders follow by the stable direction
// of the recurrence: upward for I, downward for K.
// All are element-wise over real arrays; ints and float32 promote.

// BesselI0 returns the modified Bessel function of the first kind of
// order zero, I₀(x), of each element. I₀ is even with I₀(0) = 1.
func BesselI0(a *Array) (*Array, error) {
	return a.realFunc("BesselI0", besselI0)
}

// BesselI1 returns the modified Bessel function of the first kind of
// order one, I₁(x), of each element. I₁ is odd with I₁(0) = 0.
func BesselI1(a *Array) (*Array, error) {
	return a.realFunc("BesselI1", besselI1)
}

// BesselK0 returns the modified Bessel function of the second kind of
// order zero, K₀(x), of each element, defined for x > 0. K₀ diverges
// logarithmically at 0 and returns +Inf there.
func BesselK0(a *Array) (*Array, error) {
	return a.realFunc("BesselK0", besselK0)
}

// BesselK1 returns the modified Bessel function of the second kind of
// order one, K₁(x), of each element, defined for x > 0. K₁ diverges
// like 1/x at 0 and returns +Inf there.
func BesselK1(a *Array) (*Array, error) {
	return a.realFunc("BesselK1", besselK1)
}

// BesselIn returns the modified Bessel function of the first kind of
// integer order n, Iₙ(x), of each element. The recurrence runs in its
// stable direction: orders well below |x| climb upward from the
// quadrature I₀ and I₁ at O(n) cost, and orders comparable to or above
// it run Miller's downward algorithm, which is where the upward climb
// would amplify the parasitic K component the seeds carry. The
// downward walk starts above the turning point at order |x|, so it
// costs O(n + |x|) steps: for an order far below a large argument that
// is the difference between a microsecond and minutes, which is why
// the climb exists.
func BesselIn(n int, a *Array) (*Array, error) {
	if n < 0 {
		return nil, errf("BesselIn: order must be ≥ 0, got %d", n)
	}
	return a.realFunc("BesselIn", func(x float64) float64 {
		if n == 0 {
			return besselI0(x)
		}
		if n == 1 {
			return besselI1(x)
		}
		// The relative contamination the climb suffers grows roughly
		// like e^{n²/x}, and the seeds it climbs from are only as good
		// as besselI0 and besselI1 are there: below 3.75 those are
		// polynomial fits worth about 2.5e-8, so a short climb over a
		// small argument multiplies a weak seed by the condition number
		// and lands coarser than the downward walk, which renormalises
		// through the same seed once. The climb is therefore taken only
		// where its seeds are strong, from besselIUpwardMinX up, and the
		// split follows n² ≤ 4x inside that: above it the downward walk
		// is both accurate and affordable, since the order is then
		// comparable to the argument and O(n + |x|) steps is the
		// answer's own price.
		ax := math.Abs(x)
		if ax >= besselIUpwardMinX && float64(n)*float64(n) <= 4*ax {
			return besselIUpward(n, x)
		}
		return besselInMiller(n, x)
	})
}

// besselIUpwardMinX is the argument from which the upward climb beats
// Miller's downward walk on accuracy as well as on cost: the asymptotic
// I₀ and I₁ it climbs from are near full precision there, while the
// polynomial pair below 3.75 is not. Miller's own cost is O(n + |x|), so
// staying with it under this bound costs nothing that matters.
const besselIUpwardMinX = 12

// besselIUpward climbs I₀ and I₁ to order n by the upward recurrence
// Iₖ₊₁ = Iₖ₋₁ − (2k/x)·Iₖ (the modified kind carries the minus sign
// where J carries a plus), the stable direction while the order stays
// well below |x|: the parasitic K component the seeds carry decays
// relative to I there, whereas near and above the argument the same
// recurrence amplifies it, which is what besselInMiller's downward walk
// avoids. The quadrature seeds carry the absolute scale and the parity,
// so a negative argument needs no extra handling. The caller guarantees
// x ≠ 0, 2 ≤ n, n² ≤ 4|x| and |x| ≥ besselIUpwardMinX.
func besselIUpward(n int, x float64) float64 {
	i0, i1 := besselI0(x), besselI1(x)
	// The quadrature overflows to Inf once |x| passes about 714, where
	// the true value overflows too: report that instead of letting the
	// recurrence fold Inf - Inf into NaN. The sign follows the parity
	// Iₙ(−x) = (−1)ⁿ Iₙ(x).
	if math.IsInf(i0, 0) || math.IsInf(i1, 0) {
		if x < 0 && n%2 == 1 {
			return math.Inf(-1)
		}
		return math.Inf(1)
	}
	for k := 1; k < n; k++ {
		i0, i1 = i1, i0-2*float64(k)/x*i1
	}
	return i1
}

// besselInMiller evaluates I_n by the downward Miller recurrence. The
// upward recurrence for I is unstable once the order passes the
// argument (the terms subtract nearly equal neighbours), while the
// downward direction amplifies no rounding. The arbitrary seed scale
// is removed against the exact I_0 at the end.
func besselInMiller(n int, x float64) float64 {
	if x == 0 {
		if n == 0 {
			return 1
		}
		return 0
	}
	// |x| keeps the start above n for negative arguments too: a start
	// at or below n would never pass the order on the way down and
	// renormalise against a garbage seed.
	start := n + int(math.Abs(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
		}
		// Downward recurrence: j_{k−1} = j_{k+1} + 2k/x·j_k.
		jp, j = j, jp+2*float64(k)/x*j
		if aj := math.Abs(j); aj > 1e200 {
			// I grows monotonically down the walk, and the unscaled
			// seed overflows for tiny x or high orders whose true
			// value is representable; a common rescale cancels in the
			// final ratio.
			jp /= aj
			j /= aj
			jn /= aj
		}
	}
	// j now holds the unscaled I_0. The quotient is formed before the
	// multiplication: above about x = 300 both walk values pass 1e100,
	// and the product I₀·jₙ overflows even though the answer, which is
	// the quotient times the seed scale, is an ordinary number. Taking
	// the product first turned every such evaluation into +Inf.
	return besselI0(x) * (jn / j)
}

// BesselKn returns the modified Bessel function of the second kind of
// integer order n, Kₙ(x), of each element, by the upward recurrence
// Kₙ₊₁ = Kₙ₋₁ + 2n/x·Kₙ from K₀ and K₁, which is the stable direction
// for the second-kind functions. K is defined for x > 0 at every
// order: outside that domain the answer is not order-dependent but
// undefined, so any element that is NaN or non-positive is an error
// naming the element, the same contract BesselY enforces at one
// point, never a silent NaN or a spurious +Inf.
func BesselKn(n int, a *Array) (*Array, error) {
	if n < 0 {
		return nil, errf("BesselKn: order must be ≥ 0, got %d", n)
	}
	if a.dt == Complex {
		return nil, errf("BesselKn: complex arrays are not supported")
	}
	for i := range a.Len() {
		if x := a.floatAt(i); math.IsNaN(x) || x <= 0 {
			return nil, errf("BesselKn: the argument must be positive, got %g at element %d", x, i)
		}
	}
	return a.realFunc("BesselKn", func(x float64) float64 {
		k0, k1 := besselK0(x), besselK1(x)
		if n == 0 {
			return k0
		}
		if n == 1 {
			return k1
		}
		km := k0
		k := k1
		for j := 1; j < n; j++ {
			km, k = k, km+2*float64(j)/x*k
		}
		return k
	})
}

// besselI0 evaluates I₀: the standard 3.75-piecewise polynomial fit
// below, and the exact integral representation I₀ = (1/π)∫₀^π
// e^(x·cos θ) dθ above, where the polynomial's accuracy degrades.
func besselI0(x float64) float64 {
	ax := math.Abs(x)
	if ax < 3.75 {
		t := x / 3.75
		t2 := t * t
		return 1 + t2*(3.5156229+t2*(3.0899424+t2*(1.2067492+t2*(0.2659732+t2*(0.0360768+t2*0.0045813)))))
	}
	return quadratureBesselI(x, false)
}

// besselI1 evaluates I₁: the polynomial fit below 3.75, the integral
// representation I₁ = (1/π)∫₀^π e^(x·cos θ)·cos θ dθ above. The odd
// parity comes out of the quadrature by itself; the small branch
// applies the sign of x explicitly.
func besselI1(x float64) float64 {
	ax := math.Abs(x)
	var r float64
	if ax < 3.75 {
		t := x / 3.75
		t2 := t * t
		r = ax * (0.5 + t2*(0.87890594+t2*(0.51498869+t2*(0.15084934+t2*(0.02658733+t2*(0.00301532+t2*0.00032411))))))
		if x < 0 {
			r = -r
		}
		return r
	}
	return quadratureBesselI(x, true)
}

// quadratureBesselI evaluates I₀ (weight false) or I₁ (weight true)
// through the exact integral representations over θ ∈ [0, π] by
// composite Simpson with 2048 intervals, which resolves the peak at
// θ = 0 for any practical |x|.
func quadratureBesselI(x float64, weightCos bool) float64 {
	const intervals = 2048
	h := math.Pi / float64(intervals)
	kern := func(th float64) float64 {
		v := math.Exp(x * math.Cos(th))
		if weightCos {
			v *= math.Cos(th)
		}
		return v
	}
	sum := kern(0) + kern(math.Pi)
	for i := 1; i < intervals; i++ {
		w := 2.0
		if i%2 == 1 {
			w = 4.0
		}
		sum += w * kern(float64(i)*h)
	}
	return sum * h / (3 * math.Pi)
}

// besselK0 evaluates K₀(x) = ∫₀^∞ e^(−x·cosh t) dt by composite
// Simpson quadrature over a truncated domain. The tail beyond the
// cutoff falls below 1e-100 for every x > 0, and the integrand is
// smooth, so the quadrature carries roughly ten significant digits.
// That is the deliberate trade: a provably correct entry-level K
// without hand-typed approximation coefficients.
func besselK0(x float64) float64 {
	if x <= 0 {
		return math.Inf(1)
	}
	return integrateCoshKernel(x, func(_ float64, coshU float64) float64 {
		return math.Exp(-x * coshU)
	})
}

// besselK1 evaluates K₁(x) = ∫₀^∞ e^(−x·cosh t)·cosh t dt, the
// derivative partner of K₀, by the same quadrature.
func besselK1(x float64) float64 {
	if x <= 0 {
		return math.Inf(1)
	}
	return integrateCoshKernel(x, func(_ float64, coshU float64) float64 {
		return math.Exp(-x*coshU) * coshU
	})
}

// integrateCoshKernel integrates ∫₀^∞ kernel(u; x) du where the
// kernel carries the factor e^(−x·cosh u). The cutoff follows from
// requiring the tail e^(−x·e^U/2) below 1e-100; for large x the
// result legitimately underflows to zero like the true value.
func integrateCoshKernel(x float64, kernel func(float64, float64) float64) float64 {
	const intervals = 2001
	u := math.Log(200/x) + 2
	if u < 2 {
		u = 2
	}
	h := u / float64(intervals)
	sum := kernel(0, 1) + kernel(u, math.Cosh(u))
	for i := 1; i < intervals; i++ {
		tu := float64(i) * h
		w := 2.0
		if i%2 == 1 {
			w = 4.0
		}
		sum += w * kernel(tu, math.Cosh(tu))
	}
	return sum * h / 3
}