301 lines
11 KiB
Go
301 lines
11 KiB
Go
// 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
|
||
}
|