wikipedia: bessel function

there are multiple kinds of bessel functions that are grouped together just because they are similar. they are called first kind, second kind, and so on.

x values for which the bessel function of the first kind crosses zero:

• order 0: 2.404825557695773, 5.520078110286311, 8.653727912911013
• order 1: 3.8317059702075125, 7.015586669815619, 10.173468135062722
• order 2: 5.135622301840682556301, 8.417244140399864857784, 11.61984117214905942709
• order 3: 6.380161895923983506237, 9.761023129981669678545, 13.01520072169843441983

• libc (standard) has jn(n, x) and more
• gsl has gsl_sf_bessel_Jn(n, x)
• implementation in gsl see specfunc/bessel_Jn.c
• implementation of jn in glibc see e_jnl.c in the directories sysdeps/ieee754/dbl-64 and similar (128 bit implementation in ldbl-128)
• j. bremer, 2017, an algorithm for the rapid numerical evaluation of bessel functions of real orders and arguments (code)
• leung, po kin; gammie, charles f.; noble, scott c. - bessel: fast bessel function jn(z) routine for large n,z (code)
• j. harrison, 2009, fast and accurate bessel function computation
• digital library of mathematical functions: methods of computation
• bessel programs in c/c++
• bessel.c (1998) approximations based on ideas from "Numerical Recipes, Press et. al."

bessel(order, x) = sum(n, 0, infinity,
  ((-1) ** n / (factorial(n) * gamma(n + order + 1)))
  * (x / 2) ** (2 * n + order))

• for integers, the gamma function can be replaced with the factorial
• the most direct translation to programming code quickly reaches the limits of numerical representation because of the exponential (2 * n + order)

this version is simple and compact. however, in software, it only gives accurate values for about x < 36 and then starts to rapidly move towards infinity.

bessel(order, x) = (x / 2) ** order
  * sum(n, 0, infinity, (-1 ** n * (x ** 2 / 4) ** n) / (n! * (n + order)!))

-1 * (x ** 2 / 4) / (n * (order + n))

(define (factorial n) (if (> 1 n) 1 (* n (factorial (- n 1)))))

(define (bessel order x term-count)
  (* (expt (/ x 2) order)
    (let loop ((n 1) (term (/ 1 (factorial order))))
      (if (< n term-count)
        (+ term (loop (+ 1 n) (* term -1 (/ (/ (* x x) 4) (* n (+ order n))))))
        term))))

factorial = (n) -> if 1 > n then 1 else n * factorial(n - 1)

bessel = (order, x, term_count) ->
  term = 1 / factorial order
  sum = term
  for n in [1...term_count]
    term = term * -1 * (x * x / 4) / (n * (order + n))
    sum += term
  sum * Math.pow(x / 2, order)

this version also works well for larger values. it was extracted and translated from SheetJs/bessel.

bessel_j = do ->
  M = Math
  W = 2 / M.PI

  # Coefficients for small and large x approximations for J0(x) and J1(x)
  coeff_j0_small_x_numerator = [57568490574.0, -13362590354.0, 651619640.7, -11214424.18, 77392.33017, -184.9052456].reverse()
  coeff_j0_small_x_denominator = [57568490411.0, 1029532985.0, 9494680.718, 59272.64853, 267.8532712, 1.0].reverse()
  coeff_j0_large_x_cosine = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse()
  coeff_j0_large_x_sine = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7].reverse()
  coeff_j1_small_x_numerator = [72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606].reverse()
  coeff_j1_small_x_denominator = [144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0].reverse()
  coeff_j1_large_x_cosine = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse()
  coeff_j1_large_x_sine = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse()

  # Horner's method for polynomial evaluation
  horner = (coefficients, variable) ->
    result = 0
    for coefficient in coefficients
      result = variable * result + coefficient
    result

  # Bessel function recurrence relation
  bessel_recurrence = (x, n, j0, j1, sign) ->
    return j0 if n == 0
    return j1 if n == 1
    two_over_x = 2 / x
    current_j = j1
    for order in [2..n]
      current_j = j1 * (order - 1) * two_over_x + sign * j0
      j0 = j1
      j1 = current_j
    current_j

  # Bessel function J0(x) computation
  bessel_j0 = (x) ->
    y = x * x
    if x < 8
      numerator = horner(coeff_j0_small_x_numerator, y)
      denominator = horner(coeff_j0_small_x_denominator, y)
      numerator / denominator
    else
      xx = x - 0.785398164  # x - pi/4
      y = 64 / y
      cosine_term = horner(coeff_j0_large_x_cosine, y)
      sine_term = horner(coeff_j0_large_x_sine, y)
      M.sqrt(W / x) * (M.cos(xx) * cosine_term - (M.sin(xx) * sine_term * 8 / x))

  # Bessel function J1(x) computation
  bessel_j1 = (x) ->
    y = x * x
    if M.abs(x) < 8
      numerator = x * horner(coeff_j1_small_x_numerator, y)
      denominator = horner(coeff_j1_small_x_denominator, y)
      numerator / denominator
    else
      xx = M.abs(x) - 2.356194491  # abs(x) - 3*pi/4
      y = 64 / y
      cosine_term = horner(coeff_j1_large_x_cosine, y)
      sine_term = horner(coeff_j1_large_x_sine, y)
      result = M.sqrt(W / M.abs(x)) * (M.cos(xx) * cosine_term - (M.sin(xx) * sine_term * 8 / M.abs(x)))
      result = -result if x < 0
      result

  # Main Bessel function J_n(x)
  (n, x) ->
    n = M.round(n)
    if !isFinite(x)
      return if isNaN(x) then x else 0
    if n < 0
      return (if n % 2 then -1 else 1) * bessel_j(x, -n)
    if x < 0
      return (if n % 2 then -1 else 1) * bessel_j(-x, n)
    return bessel_j0(x) if n == 0
    return bessel_j1(x) if n == 1
    return 0 if x == 0
    if x > n
      bessel_recurrence(x, n, bessel_j0(x), bessel_j1(x), -1)
    else
      m = 2 * M.floor((n + M.floor(M.sqrt(40 * n))) / 2)
      sum_sign = false
      previous_j = 0.0
      sum_j = 0.0
      current_j = 1.0
      previous_j_minus_1 = 0.0
      two_over_x = 2 / x
      for j in [m..1]
        previous_j_minus_1 = j * two_over_x * current_j - previous_j
        previous_j = current_j
        current_j = previous_j_minus_1
        if M.abs(current_j) > 1e10
          current_j *= 1e-10
          previous_j *= 1e-10
          sum_j *= 1e-10
        if sum_sign
          sum_j += current_j
        sum_sign = !sum_sign
        if j == n
          previous_j = previous_j_minus_1
      sum_j = 2.0 * sum_j - current_j
      previous_j / sum_j