JuliaMath / JuliaMath/Bessels.jl

FEAT request: compute bessel function across multiple n

Open
#135 1 comment 0 reactions 0 assignees View on GitHub
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

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.