Karim BELABAS on Fri, 12 Nov 1999 20:09:40 +0100 (MET)

 Re: dilog() bug

```[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
{
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)
{
}
!   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));
!   if ((m&1) == 0) y = gneg_i(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 = gneg_i(p1); tetpil=avma;
}

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));
!     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)
!     if (m&1) p1 = gmul(logx,p1); else y = gneg_i(y);
}
!
}

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
--