Re: besseli: Am I missing something?

[John W. Eaton]

On 12-Jan-2012, Robert T. Short wrote:

| Robert T. Short wrote:
| > Before I file a bug report, maybe I should check that I am not all wet.
| >
| > If I recall my Bessel function theory,
| >
| > besseli(n,x) = besseli(-n,x)
| >
| > if n is an integer.
| >
| > In octave 3.4.2 I get
| >
| > octave:3> besseli(1,1),besseli(-1,1)
| > ans =  0.565159103992485
| > ans =  0.565159103992485
| >
| > Looks good for n=1, but for n=10
| >
| > octave:4> besseli(10,1),besseli(-10,1)
| > ans =  2.75294803983687e-10
| > ans = -1.40610343427994e-07
| >
| > Not so good!
| >
| > And for n=100
| >
| > octave:5> besseli(100,1),besseli(-100,1)
| > ans =  8.47367400813812e-189
| > ans =  7.38028373423800e+170
| >
| > BTW, the numbers for positive n agree to about 9 or 10 decimal places 
| > with my tables.  Also, BTW, the relationships for besselj seem to work 
| > fine.
| >
| > Anybody?
| >
| > Bob
| >
| >
| Well, I will file the report.  I will also make some tests for this and 
| other cases.
| 
| I have already determined that the octave implementation is only good to 
| maybe 9 places in other situations.  Since Amos used asymptotic 
| expansions in some cases, I can't imagine how we would get more, but I 
| am not a numerical analysist.  Since I am messing with this stuff, I 
| will dig deeper and report back.  Thanks for the comments and pointers.

Does the following change fix the immediate problem of bad results
from besseli with negative integer orders?

diff --git a/liboctave/lo-specfun.cc b/liboctave/lo-specfun.cc
--- a/liboctave/lo-specfun.cc
+++ b/liboctave/lo-specfun.cc
@@ -852,6 +852,15 @@
 
       retval = bessel_return_value (Complex (yr, yi), ierr);
     }
+  else if (is_integer_value (alpha))
+    {
+      // zbesi can overflow as z->0, and cause troubles for generic case below
+      alpha = -alpha;
+      Complex tmp = zbesi (z, alpha, kode, ierr);
+      if ((static_cast <long> (alpha)) & 1)
+        tmp = - tmp;
+      retval = bessel_return_value (tmp, ierr);
+    }
   else
     {
       alpha = -alpha;

This is the same sort of special case already used for besselj.

There is a similar special case in bessely for negative half orders,
but I'm not sure that it is doing the right thing.

jwe