It seems to me that this has been requested from time to time.
Here's a version that should work..
push(c,l)::=buildq([c,l],l:cons(c,l))$ /*utility*/
/* Compute Bessel J[n](x) for non-negative integer n and float or
bigfloat real x*/
bjn(n,x):=block([lam:0, r:[1,0],xi:1/x,li:0],
if ?floatp(x) then bessel_j(n,x) else
(for k:n+findn(x,fpprec)-2 step -1 thru 0 do
(push(2*(k+1)*xi*first(r)-second(r), r),
if evenp(k) then lam:lam+first(r)),
r[n+1]/(2*lam-first(r))))$
/* findn return a recommended value for iteration count in the bj routine,
given the argument x and d=number of decimal digits */
findn(x,d):=block([r:1, y:[1,0],lim:2*10^d],
while abs(first(y))<lim do
(push((2*(r-1)/x*first(y)-second(y)), y), r:r+1),
r-1)
If someone cares to translate this into lisp or bigfloat package (Ray?)
it should be fine. This routine should be good for about fpprec digits
except for arguments very close to a zero of Bessel J. In this rare
situation, a simple fix is to request twice as many digits.
RJF