JuliaMath / JuliaMath/Bessels.jl
FEAT request: compute bessel function across multiple n
- Dominant language
- Julia
- Stars
- 90
- Forks
- 11
- PR merge metrics
- No merged PRs in 30d
Description
Below is a poor-man hand crafted version to compute a series of Besselj efficiently. I am not sure if this repo want to target at this kind of feature.
```
function besselj_ladder!(out, M::Integer, z::T) where {T} # out[k+1] = J_k(z), k=0..M
# Below √eps(T) the one-term series J_k=(z/2)^k/k! is exact to O(z²)≤eps(T) rel,
# AND it is the safe path on denormal z, where the recurrence's 2n/z → Inf.
if abs(z) < sqrt(eps(T))
half = z / 2
acc = one(T)
@inbounds out[1] = acc
@inbounds for k in 1:M
acc *= half / k
out[k + 1] = acc
end
return out
end
# Miller downward recurrence: recurse `J_{n−1}=(2n/z)J_n−J_{n+1}` (stable downward for `n>z`)
# normalize by the Neumann identity `J_0 + 2(J_2+J_4+…) = 1`.
# rescale by inv(BIG) whenever it crosses BIG to stay clear of floatmax(T).
BIG = sqrt(floatmax(T))
s = inv(BIG)
base = max(M, ceil(Int, abs(z))) # seed must clear BOTH the wanted M and the turning point n≈z
# extra steps decay the seed error in the n>z stable region; digits gained scale
# with step count, so widen the margin for higher-precision T (≥ the Float64 margin).
margin = (sqrt(40 * (base + 1)) + 15) * max(1, precision(T) / 53)
N = base + ceil(Int, margin)
fkp1 = zero(T)
fk = sqrt(floatmin(T)) # tiny seed leaves full exponent range to grow downward
nrm = zero(T)
@inbounds for n in N:-1:1
fkm1 = (2 * n / z) * fk - fkp1
(n - 1 <= M) && (out[n] = fkm1)
iseven(n - 1) && (nrm += fkm1)
fkp1, fk = fk, fkm1
if abs(fk) > BIG # preserve ratios (hence J)
fk *= s; fkp1 *= s; nrm *= s
for j in n:(M + 1)
out[j] *= s
end
end
end
invn = inv(2 * nrm - out[1]) # out[1]=J_0 is double-counted in 2·Σ_even
@inbounds for k in 0:M
out[k + 1] *= invn
end
return out
end
```
Contributor guide
No contributing guide indexed for this repository
Research direction
No files or tests are named. Start by reviewing the proposed `besselj_ladder!` implementation alongside the repository's existing Bessel-function entry points, then determine whether multi-order computation belongs in the public API. Done means an agreed integration path and validated results for multiple n values, including the small-z and recurrence cases described.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- julia
- Domain
- backend
- Issue type
- Feature
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Activity status
- Quiet
- Clarity
- Mostly clear
- Newbie friendliness
- 45/100