Re: 2.0.8 bug in besselk?

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

Leonhard Moehring wrote:
> ? besselk(I,16)
> ----> (apparently) infinite loop goes here.
> 
> besselk seems to enter an infinite loop whenever the index is
> complex and the argument is roughly between 16 and 20.
> (while it works fine with GP/PARI 1.3.9)
> 
> Actually the lines:
> 
> >   /* this 10 should really be a 5, but then kbessel(15.99) enters oo loop */
> >  lbin = 10 - bit_accuracy(l); av1=avma;
> 
> in trans3.c seem to imply that the work in this routine isn't done yet. ;)

In fact the work _was_ done in kbessel, but not in older brother hyperu which
is eventually called when index is complex.

Problem: the result is less accurate, but round-off errors prevent "optimal"
accuracy in general (we'd need interval arithmetic, and I won't do that
before the beta...).

(this patch changes slightly the output of the "trans" bench).

Karim.

*** src/basemath/trans3.c.orig	Mon May  4 14:54:57 1998
--- src/basemath/trans3.c	Wed Jun  3 11:23:58 1998
***************
*** 148,154 ****
      d=cgetr(l1); e=cgetr(l1); f=cgetr(l1);
    }
    n=(long)(bit_accuracy(l)*LOG2 + PI*sqrt(gtodouble(gabs(gmul(a,a1),l1))));
!   lbin = 5-bit_accuracy(l); av1=avma;
    if (gcmpgs(x,n)<0)
    {
      gn=stoi(n); zf=gpui(gn,gneg(a),l1);
--- 148,154 ----
      d=cgetr(l1); e=cgetr(l1); f=cgetr(l1);
    }
    n=(long)(bit_accuracy(l)*LOG2 + PI*sqrt(gtodouble(gabs(gmul(a,a1),l1))));
!   lbin = 10-bit_accuracy(l); av1=avma;
    if (gcmpgs(x,n)<0)
    {
      gn=stoi(n); zf=gpui(gn,gneg(a),l1);
--
Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud             Tel: (33 1) 01 69 15 57 48
F-91405 Orsay (France)           Fax: (33 1) 01 69 15 60 19
--
PARI/GP Home Page: http://pari.home.ml.org