*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 will discuss in a future post, I am after `j1`

only, so… *my job here is done*.