What should bessel_j(-1.0,1) return? Currently, it gives an error
because it doesn't know that bessel_j(-1,x) = -bessel_j(1,x).
This is easy to fix.
However, what should we do when given a floating-point order of -1.0?
Treat it as if we got an integer order of -1? I think this is
reasonable, but I can see arguments where we should just leave it as
bessel_j(-1.0,1). (We can't currently numerically evaluate Bessel J
for arbitrary negative orders, and get a SLATEC error when given such
orders.)
Ray