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