| 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