From Wikipedia
Iterative implementation
I came up with the above implementation, with which I plotted these.
J0-4 in [0..2pi]
2*abs(J1-4(x)/x) in [0..2pi]
Iterative methods are a terrible choice for many use cases because of error accumulation, jitter, the very presence of a loop, etc… But in this case this method also happens to be pretty weak and way too easy to break:
J0-4 in [0..8pi]
Deus ex machina
After some unsuccessful googling I was about to work on my own Taylor Series approximation or at least study the proposed implementations in Numerical Recipes. But I learnt to my surprise :o) that <math.h> provides j0/j1/jn/y0/y1/yn and even CUDA provides j0f/j1f/jnf/y0f/y1f/ynf.
I checked whether the C/CUDA implementations match my iterative implementation, and they do.
For reasons I may discuss in a future post, I am after j1 only, so… my job here is done.