Hello Robert,
I had a second look at the code. This part of the code is in my opinion
definitley the part which does the correct calculation for GCL. I don't see any
other possibility. I have tested it directly. My problem is that I am not able
to reproduce the wrong results.
In the code in the email there are missing two parenthesis (It is allways a bit
dangerous to reformat code in an email). This should be more correct:
(let* ((arg (float arg))
(order (float order))
(result (* (expt arg order)
(expt (complex 0.0 arg) (- order))
(bessel-j order (complex 0.0 arg)))))
;; Try to clean up result if we know the result is
;; purely real or purely imaginary.
(cond ((>= arg 0)
...
In the testsuite there are pairs of values which test this part of code:
bessel_i(-1, -2.0) and bessel_i(-1, +2.0)
bessel_i(-1.5,-1.5) and bessel_i(-1.5,+1.5)
bessel_i(-1.8,-1.5) and bessel_i(-1.8,+1.5)
bessel_i(-2 ,-1.5) and bessel_i(-2, +1.5)
I don't see why the examples above work but not:
bessel_i(-2.5,-1.5) and bessel_i(-2.5,+1.5)
I have the same problem with the third wrong result. Here also pairs are tested
with an order 1.5, 1.8, 2.0, 2.5 and an argument +1.0 and -1.0. Only the result
for bessel_i(2.5,-1.0) is wrong. What is special with the number "2.5".
Here are some more pairs calculated with GCL and checked with
functions.wolfram.com:
(%i6) bessel_i(-3.0,-1.5);
(%o6) - 0.080774113016092
(%i7) bessel_i(-3.0,+1.5);
(%o7) 0.080774113016092
(%i8) bessel_i(+3.0,-1.5);
(%o8) - 0.080774113016092
(%i9) bessel_i(+3.0,+1.5);
(%o9) 0.080774113016092
(%i10) bessel_i(-3.3,-1.5);
(%o10) 0.80626434106559 - 1.109727662163278 %i
(%i11) bessel_i(-3.3,+1.5);
(%o11) - 1.371698826945733
(%i12) bessel_i(+3.3,-1.5);
(%o12) - 0.040233131323362 %i - 0.029231080941247
(%i13) bessel_i(+3.3,+1.5);
(%o13) 0.049730885263352
(%i14) bessel_i(-3.5,-1.5);
(%o14) - 2.306410036086646 %i
(%i15) bessel_i(-3.5,+1.5);
(%o15) - 2.306410036086646
(%i16) bessel_i(+3.5,-1.5);
(%o16) - 0.035543108479936 %i
(%i17) bessel_i(+3.5,+1.5);
(%o17) 0.035543108479936
Dieter Kaiser
-----Urspr?ngliche Nachricht-----
Von: robert.dodier at gmail.com [mailto:robert.dodier at gmail.com]
Gesendet: Montag, 11. August 2008 03:02
An: Dieter Kaiser
Cc: maxima at math.utexas.edu
Betreff: Re: [Maxima] Maxima 5.16.1 release
On 8/10/08, Dieter Kaiser <drdieterkaiser at web.de> wrote:
> Perhaps the following code does it better:
>
> (let* ((arg (float arg)
> (order (float order)))
> (result (* (expt arg order)
> (expt (complex 0.0 arg) (- order))
> (bessel-j order (complex 0.0 arg))))
Unfortunately putting in this code (with a minor change --- I think
there is a parenthesis missing from the end of the first line)
doesn't seem to change the failed tests.
best
Robert Dodier