Karim Belabas on Fri, 28 Nov 1997 05:28:09 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
patch9 |
* addresses a bug reported by Bill Dany (directly to us, patch by Henri Cohen). > (05:19) gp > psi(1+I) > %1 = 0.09465031687201272145140085631 + 1.076674036798763527702729101*I whereas the correct value would be: 0.09465032062247697727187848267 + 1.076674047468581174134050793*I (not that far, but still off...). This slight inaccuracy could affect the functions psi and lngamma for non-real arguments. It was already here in version 1.39. Cheers, Karim. ============================ patch9 (2.0.alpha) ============================ *** src/basemath/trans2.c.orig Fri Nov 14 04:53:32 1997 --- src/basemath/trans2.c Thu Nov 27 19:28:35 1997 *************** *** 1251,1261 **** { for (i=1; i<=n; i++) { ! addsrz(-1,p2,p2); p7 = (i==1)? p2: mulrr(p7,p2); } f=signe(p7); if (f<0) setsigne(p7,1); subrrz(p4,mplog(p7),p4); } if (u) { setlg(pitemp,l+1); p1 = mulrr(pitemp,x); --- 1251,1262 ---- { for (i=1; i<=n; i++) { ! addsrz(-1,p2,p2); p7 = (i==1)? rcopy(p2): mulrr(p7,p2); } f=signe(p7); if (f<0) setsigne(p7,1); subrrz(p4,mplog(p7),p4); } + else f=1; if (u) { setlg(pitemp,l+1); p1 = mulrr(pitemp,x); *************** *** 1362,1368 **** for (i=1; i<=n; i++) { addsrz(-1,(GEN)p2[1], (GEN)p2[1]); ! if (i==1) p7[1] = (long)p1; else p7 = gmul(p7,p2); } gsubz(p4,glog(p7,l+1), p4); } --- 1363,1369 ---- for (i=1; i<=n; i++) { addsrz(-1,(GEN)p2[1], (GEN)p2[1]); ! if (i==1) p7[1] = lrcopy((GEN)p2[1]); else p7 = gmul(p7,p2); } gsubz(p4,glog(p7,l+1), p4); } *************** *** 1532,1538 **** cxpsi(GEN z, long prec) /* version p=2 */ { long l,n,k,x,xx,av,av1,tetpil; ! GEN zk,u,v,a,b; l = precision(z); if (!l) l = prec; av=avma; x=(long)(1 + (bit_accuracy(l)>>1)*LOG2 + 1.58*rtodbl(gabs(z,DEFAULTPREC))); --- 1533,1539 ---- cxpsi(GEN z, long prec) /* version p=2 */ { long l,n,k,x,xx,av,av1,tetpil; ! GEN zk,u,v,a,b,p1; l = precision(z); if (!l) l = prec; av=avma; x=(long)(1 + (bit_accuracy(l)>>1)*LOG2 + 1.58*rtodbl(gabs(z,DEFAULTPREC))); *************** *** 1541,1547 **** b=cgetg(3,t_COMPLEX); b[1]=lgetr(l); b[2]=lgetr(l); gaffsg(1,b); u=cgetg(3,t_COMPLEX); u[1]=lgetr(l); u[2]=lgetr(l); v=cgetg(3,t_COMPLEX); v[1]=lgetr(l); v[2]=lgetr(l); gaffsg(1,v); ! a=glog(a,l); gaffect(a,u); av1=avma; for (k=1; k<=n; k++) { zk=(k>1) ? gaddsg(k-1,z) : z; --- 1542,1548 ---- b=cgetg(3,t_COMPLEX); b[1]=lgetr(l); b[2]=lgetr(l); gaffsg(1,b); u=cgetg(3,t_COMPLEX); u[1]=lgetr(l); u[2]=lgetr(l); v=cgetg(3,t_COMPLEX); v[1]=lgetr(l); v[2]=lgetr(l); gaffsg(1,v); ! p1=glog(a,l); gaffect(p1,a); gaffect(p1,u); av1=avma; for (k=1; k<=n; k++) { zk=(k>1) ? gaddsg(k-1,z) : z; -- Karim Belabas e-mail: Max-Planck-Institut fuer Mathematik karim@mpim-bonn.mpg.de Gottfried-Claren-Str. 26 tel: 53225 Bonn (Germany) (00 49 228) 402-245