Karim BELABAS on Wed, 3 Jun 1998 12:42:52 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: 2.0.8 bug in besselk? |
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