bessel_i quadpack bug 5.26.0



On Feb. 29, 2012, Raymond Toy wrote

>Still need to make bessel_i(n,%i*float)
>return a purely imaginary result, though.  
>bessel_j apparently already does this.

I use the blunt force tool fchop in my work:

-----------------------------------
(%i1) cfloat(p):= expand(float(rectform(p)))$
(%i2) _small% : 1.0e-14$

(%i3) _chop%(ex) := 
block([eex],
  if numberp(ex) then
      (eex:float(ex),
       if ( abs(eex) < _small% ) then 0.0
       else eex)
   else ex)$

(%i4) fchop(expr) := if mapatom(expr) then _chop%(expr)
            else scanmap(_chop%,expr)$

(%i5) for n thru 10 do
    print(n,"   ",cfloat(bessel_i(n,%i*2)))$

1     0.57672480775687*%i+3.5313043199302822E-17 
2     4.3208279328745645E-17*%i-0.35283402861564 
3     -0.1289432494744*%i-2.3685708388328505E-17 
4     0.033995719807568-8.3262748958227084E-18*%i 
5     0.0070396297558717*%i+2.1551938318993842E-18 
6     4.4175064766224005E-19*%i-0.00120242897179 
7     -1.7494407486827413E-4*%i-7.4983168976543158E-20 
8     2.2179552287925885E-5-1.0864488262691872E-20*%i 
9     2.4923434351330654E-6*%i+1.3734628229571489E-21 
10     1.5401790120662821E-22*%i-2.5153862827167347E-7 

(%i6) for n thru 10 do
    print(n,"   ",fchop(cfloat(bessel_i(n,%i*2))))$
1     0.57672480775687*%i 
2     -0.3528340286156 
3     -0.1289432494744*%i 
4     0.033995719807568 
5     0.0070396297558717*%i 
6     -0.00120242897179 
7     -1.7494407486827413E-4*%i 
8     2.2179552287925885E-5 
9     2.4923434351330654E-6*%i 
10     -2.515386282716735E-7 
------------------------------------------------------
Ted