Re: dilog() bug

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

[Michael Somos:]
> gp> polylog(2,10)
> %1 = 0.5363012873578627365501597698 - 7.233784412415464812490046550*I
> gp> dilog(10) 
> %2 = 0.5363012873578627365501597699 - 7.233784412415464812490046550*I
> gp> \ps3
>    seriesprecision = 3 significant terms
> gp> polylog(2,x)
> %3 = x + 1/4*x^2 + O(x^3)
> gp> dilog(x)
>   ***   incorrect type in comparison.
> 
> It seems to me that polylog(2,x) and dilog(x) shuould give essentially
> the same answer for any value of x. In fact, I would define it as follows:
> 
> dilog(x)=polylog(2,x)
> 
> Any reason why this is not the case?

The functional equation (relating values at x and 1/x) can be written in a
(slightly) more efficient way when m = 2. But in fact, the way it was
written, it became (slightly) _less_ efficient (to save 3 multiplications,
one ordinary logarithm had to be recomputed...).

I trivialized dilog() the way you suggested and special cased m = 2 in
polylog(), rewriting the optimization so that it works as intended (skipping
the log computation). The patch is big, since I cleaned up the neighbouring
code.

   Karim.

Index: src/basemath/trans3.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/basemath/trans3.c,v
retrieving revision 1.2
diff -c -r1.2 trans3.c
*** src/basemath/trans3.c	1999/09/23 17:50:57	1.2
--- src/basemath/trans3.c	1999/11/12 19:05:00
***************
*** 809,885 ****
  GEN
  polylog(long m, GEN x, long prec)
  {
!   long av,av1,limpile,tetpil,l,e,sx,i;
!   GEN p1,p2,n,y,logx,logx2;
!   GEN *gptr[4];
  
    if (m<0) err(talker,"negative index in polylog");
    if (!m) return gneg(ghalf);
    if (gcmp0(x)) return gcopy(x);
!   av=avma;
    if (m==1)
!   {
!     p1=glog(gsub(gun,x),prec); tetpil=avma;
!     return gerepile(av,tetpil,gneg(p1));
!   }
!   l=precision(x);
    if (!l) { l=prec; x=gmul(x, realun(l)); }
!   e=gexpo(gnorm(x)); if (!e || e== -1) return cxpolylog(m,x,prec);
!   if (e>0)
    {
!     logx=glog(x,l); sx=gsigne(gimag(x));
      if (!sx)
      {
!       if (m&1)
!         sx = gsigne(greal(gsub(gun,x)));
!       else
!         sx = -gsigne(greal(x));
      }
!     x=ginv(x);
    }
    av1=avma; limpile=stack_lim(av1,1);
!   y=x; n=gun; p1=x;
!   do
    {
!     n=addsi(1,n); p1=gmul(x,p1); p2=gdiv(p1,gpuigs(n,m));
!     y=gadd(y,p2);
      if (low_stack(limpile, stack_lim(av1,1)))
!     {
        if(DEBUGMEM>1) err(warnmem,"polylog");
!       gptr[0]=&y; gptr[1]=&n; gptr[2]=&p1; gptr[3]=&p2;
!       gerepilemany(av1,gptr,4);
      }
    }
!   while (gexpo(p2) > -bit_accuracy(l));
!   tetpil=avma;
!   if (e<=0) return gerepile(av,tetpil,gcopy(y));
!   logx2=gsqr(logx); p1=gneg_i(ghalf);
!   for (i=m-2; i>=0; i-=2)
    {
!     p1 = gadd(izeta(m-i,l), gmul(p1,gdivgs(logx2,(i+1)*(i+2))));
    }
!   if (m&1) p1 = gmul(logx,p1);
!   p2=cgetg(3,t_COMPLEX); p2[1]=zero;
!   p2[2]=ldiv(gmulsg(sx,mppi(l)),mpfact(m-1));
!   p1=gadd(gmul2n(p1,1), gmul(p2,gpuigs(logx,m-1)));
!   if ((m&1) == 0) y = gneg_i(y);
!   tetpil=avma; return gerepile(av,tetpil,gadd(p1,y));
  }
  
  GEN
  dilog(GEN x, long prec)
  {
!   GEN p1,p2,p3,y;
!   long av=avma,tetpil;
! 
!   if (gcmpgs(gnorm(x),1)<1)
!   {
!     tetpil=avma; return gerepile(av,tetpil,polylog(2,x,prec));
!   }
!   y=gneg_i(polylog(2,ginv(x),prec)); p3=mppi(prec); p2=gdivgs(gsqr(p3),6);
!   p1=glog(gneg_i(x),prec); p1=gadd(p2,gmul2n(gsqr(p1),-1));
!   p1 = gneg_i(p1); tetpil=avma;
!   return gerepile(av,tetpil,gadd(y,p1));
  }
  
  GEN
--- 809,883 ----
  GEN
  polylog(long m, GEN x, long prec)
  {
!   long av,av1,limpile,l,e,i,G;
!   GEN z,p1,p2,n,y,logx;
  
    if (m<0) err(talker,"negative index in polylog");
    if (!m) return gneg(ghalf);
    if (gcmp0(x)) return gcopy(x);
!   av = avma;
    if (m==1)
!     return gerepileupto(av, gneg(glog(gsub(gun,x), prec)));
! 
!   l = precision(x);
    if (!l) { l=prec; x=gmul(x, realun(l)); }
!   e = gexpo(gnorm(x)); if (!e || e== -1) return cxpolylog(m,x,prec);
!   if (e > 0)
    {
!     long sx = gsigne(gimag(x));
!     logx = glog(x,l);
      if (!sx)
      {
!       if (m&1) sx = gsigne(gsub(gun,greal(x)));
!       else     sx = - gsigne(greal(x));
      }
!     z = cgetg(3,t_COMPLEX);
!     z[1] = zero;
!     z[2] = ldivri(mppi(l), mpfact(m-1));
!     if (sx < 0) z[2] = lnegr((GEN)z[2]);
!     x = ginv(x);
    }
+   G = -bit_accuracy(l);
+   n = icopy(gun);
    av1=avma; limpile=stack_lim(av1,1);
!   y=x; p1=x;
!   for (i=2; ; i++)
    {
!     n[2] = i; p1 = gmul(x,p1); p2 = gdiv(p1,gpuigs(n,m));
!     y = gadd(y,p2);
!     if (gexpo(p2) <= G) break;
! 
      if (low_stack(limpile, stack_lim(av1,1)))
!     { GEN *gptr[2]; gptr[0]=&y; gptr[1]=&p1;
        if(DEBUGMEM>1) err(warnmem,"polylog");
!       gerepilemany(av1,gptr,2);
      }
    }
!   if (e < 0) return gerepileupto(av, y);
! 
!   if (m == 2)
!   { /* same formula as below, written more efficiently */
!     y = gneg_i(y);
!     p1 = gmul2n(gsqr(gsub(logx, z)), -1); /* = (log(-x))^2 / 2 */
!     p1 = gadd(divrs(gsqr(mppi(l)), 6), p1);
!     p1 = gneg_i(p1);
!   }
!   else
    {
!     GEN logx2 = gsqr(logx); p1 = gneg_i(ghalf);
!     for (i=m-2; i>=0; i-=2)
!       p1 = gadd(izeta(m-i,l), gmul(p1,gdivgs(logx2,(i+1)*(i+2))));
!     if (m&1) p1 = gmul(logx,p1); else y = gneg_i(y);
!     p1 = gadd(gmul2n(p1,1), gmul(z,gpuigs(logx,m-1)));
    }
! 
!   return gerepileupto(av, gadd(y,p1));
  }
  
  GEN
  dilog(GEN x, long prec)
  {
!   return gpolylog(2, x, prec);
  }
  
  GEN
__
Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud             Tel: (00 33) 1 69 15 57 48
F-91405 Orsay (France)           Fax: (00 33) 1 69 15 60 19
--
PARI/GP Home Page: http://hasse.mathematik.tu-muenchen.de/ntsw/pari/