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