use openlibm jn and yn instead of AMOS

[stevengj]

For example, besselj(1,1.5) gives 0.5579365079100995, but bessely(1,1.5) gives a Complex result -0.4123086269739113 + 0.0im. There is no reason not to return a Real result for Real arguments.

[staticfloat]

bessely is Real for positive-real input arguments, but for negative real it turns Complex again:

julia>  bessely(3,-3)
0.5385416161050315 - 0.6181254445105032im

Therefore, I don't think we can use type-inference on the inputs to determine the types outputs.

[ivarne]

The sqrt() function solves this by throwing a DomainError if the argument is Real and negative. Thus you have to convert the argument to Complex if you want the complex result. I don't know the bessely function so this might not be intuitive for users, but I really like the approach for sqrt().

[StefanKarpinski]

I suspect that following the pattern of sqrt here is the way to go.

[stevengj]

I agree that following sqrt is the thing to do.

Note also that we are currently inconsistent: bessely0 and bessely1 give real results for positive real arguments, and NaN for negative real arguments.

(Also, it's not clear why we don't use openlibm's jn and yn functions for integer order and real arguments, rather than the slower complex AMOS routines, about 2× slower on my machine.)

[ivarne]

besseli apparently also need to be changed.