Comparison of quad_qag with bessel functions using integer or float orders

On Nov. 5, 2011, I wrote:
------------------
>We compare using integer order with using 
>float order for the realpart (bessel_u (n,x)) and
>the imagpart (bessel_u (n,x)) for u  =  j, i, k, and  y,
> when using quag_qag.
--------------------------
The motivation is to know ahead of time that
we can use either integer or float orders for
all our bessel functions. At present, we
cannot, depending on the particular bessel
function and order, and arg.

The test reported Nov. 5 was based on the
Windows binary 5.25.1gcl, without adding
the improved bessel.lisp which I downloaded
from git Oct. 29.

That Nov. 4 test showed that bessel_y and bessel_k
didn't care about how the integral order value
was supplied, but both bessel_j and bessel_i
cared, in some cases, when the arg was %i*x.

Loading in the new bessel.lisp cures the two
cases found for bessel_j, but does not cure
the seven cases found for bessel_i.

Let's use the simplified notation:

  RP(b_u( n,z)) means realpart(bessel_u(n,z))

  IP (b_u (n,z)) means imagpart (bessel_u (n,z))

  I x  means %i*x
===========================
We summarize the status first for the bessel_j
cases.

1.  IP (b_j (1,I x)) is cured by new bessel.lisp

  (%i50) quad_qag(imagpart(bessel_j(1,%i*x)),x,0,3,3,limit=700);
(%o50) [3.880792585865024,4.3085452826219503E-14,31,0]

(%i51) float(imagpart(bessel_j(1,%i*2)));
(%o51) 1.590636854637329

2. IP (b_j (3, I x)) is cured by new bessel.lisp

 (%i53) quad_qag(imagpart(bessel_j(3,%i*x)),x,0,3,3,limit=700);
(%o53) [-0.60963229599488,6.7682781156861259E-15,31,0]

(%i54) float(imagpart(bessel_j(3,%i*2)));
(%o54) -0.21273995923985

===========================
We next summarize the status for the bessel_i cases.

3. RP (b_i (0, I x))  not cured

  (%i55) quad_qag(realpart(bessel_i(0,%i*x)),x,0,3,3,limit=700);
(%o55) quad_qag(bessel_i(0,%i*x),x,0,3,3,epsrel = 1.0E-8,epsabs = 0.0,
                limit = 700)

4.  IP (b_i (0.0, I x))  not cured and still hangs

(%i56) quad_qag(imagpart(bessel_i(0.0,%i*x)),x,0,3,3,limit=700);
(%o56) quad_qag('imagpart(bessel_i(0.0,%i*x)),x,0,3,3,epsrel = 1.0E-8,
                epsabs = 0.0,limit = 700)

5. RP (b_i (1.0, I x))  not cured

 (%i58) quad_qag(realpart(bessel_i(1.0,%i*x)),x,0,3,3,limit=700);
(%o58) quad_qag('realpart(bessel_i(1.0,%i*x)),x,0,3,3,epsrel = 1.0E-8,
                epsabs = 0.0,limit = 700)

6.  IP (b_i (1, I x))  not cured and shows incorrect float behavior

 (%i60) quad_qag(imagpart(bessel_i(1,%i*x)),x,0,3,3,limit=700);
(%o60) quad_qag(-%i*bessel_i(1,%i*x),x,0,3,3,epsrel = 1.0E-8,epsabs = 0.0,
                limit = 700)

???  (%i85) expand (float(imagpart(bessel_i(1,%i*2))));

???   (%o85) 0.57672480775687-3.5313043199302822E-17*%i


7. RP (b_i (2, I x))  not cured

(%i62) quad_qag(realpart(bessel_i(2,%i*x)),x,0,3,3,limit=700);
(%o62) quad_qag(bessel_i(2,%i*x),x,0,3,3,epsrel = 1.0E-8,epsabs = 0.0,
                limit = 700)


8. IP (b_i (3, I x)) not cured and shows incorrect float behavior

(%i64) quad_qag(imagpart(bessel_i(3,%i*x)),x,0,3,3,limit=700);
(%o64) quad_qag(-%i*bessel_i(3,%i*x),x,0,3,3,epsrel = 1.0E-8,epsabs = 0.0,
                limit = 700)

???   (%i71) expand(float(imagpart(bessel_i(3,%i*2))));
???   (%o71) 2.3685708388328505E-17*%i-0.1289432494744


9.  RP (b_i (4, I x))  not cured and shows incorrect float behavior

 (%i66) quad_qag(realpart(bessel_i(4,%i*x)),x,0,3,3,limit=700);
(%o66) quad_qag(bessel_i(4,%i*x),x,0,3,3,epsrel = 1.0E-8,epsabs = 0.0,
                limit = 700)

???   (%i68) float(realpart(bessel_i(4,%i*2)));
???   (%o68) 0.033995719807568-8.3262748958227084E-18*%i
-------------------------------------------------------------

Most of the problem cases of b_i return good values (ie., not a noun form)
if the float order is used.

However two of these problem b_i cases are anomalous, in that the
float order use returns a noun form, while the integer order use
returns a good result. These two cases are #4 and #5.

Ted Woollett