\\we redefine the bernoulli numbers to simplify our formulas.

{newbern(n)=bernreal(n+1)/(n+1)!}
{k(D,n)=kronecker(D,n)}

\\we redefine the bernoulli numbers associated to a character to simplify our formulas.

{newbernchar(D,n)=sum(a1=1,abs(D)-1,k(D,a1)*sum(m=0,n+1,(abs(D)-a1)^m/m!*abs(D)^(n-m)*newbern(n-m)))}

\\g is a constant which facilitates the use of Crandalls method for multiple polylogarithms at roots of unity.

{g(D)=Pi/abs(D)}

\\q(D) is a Dth root of unity.

{q(D)=exp(2*Pi*I/abs(D))}

\\we define bernoulli-type numbers corresponding to roots of unity.

{a(z,D,n)=sum(j=1,D,z^j*sum(m=0,n+1,(D-j)^m/m!*newbern(n-m)*D^(n-m)))}



\\the following two commands compute riemann zeta and l-series values at positive integer values.

{rzeta(s,N1,N2)=(sum(n=-1,N1,newbern(n)/(n+s))+sum(n=1,N2,sum(j=0,s-1,(s-1)!/j!*exp(-n)/n^(s-j))))/(s-1)!}

{lser(s,D,N1,N2)=(sum(n=0,N1,newbernchar(D,n)*g(D)^(n+s)/(n+s))+sum(n=1,N2,k(D,n)*sum(j=0,s-1,(s-1)!/j!*g(D)^j*exp(-g(D)*n)/n^(s-j))))/(s-1)!}



\\the following three commands compute single, double, and triple polylogarithms at roots of unity by finding partial sums of Dirichlet series using Bailey acceleration.  These converge SLOWLY.

{sumsingle(x,A,N)=sum(n=1,N,x^n/n^A,0.)}

{sumdouble(x,y,A,B,N,sum1,sum2)=sum1=0.;sum2=0.;for(q1=1,N,sum2=sum2+x^q1/q1^A;sum1=sum1+y^(q1+1)/(q1+1)^B*sum2);sum1}

{sumtriple(x,y,z,A,B,C,N,sum1,sum2,sum3)=sum1=0.;sum2=0.;sum3=0.;for(q1=1,N,sum3=sum3+x^q1/q1^A;sum2=sum2+y^(q1+1)/(q1+1)^B*sum3;sum1=sum1+z^(q1+2)/(q1+2)^C*sum2);sum1}



\\dbllseries computes the double l-value L(\chi,1;s1,s2) for s2>1.  Bailey acceleration has been incorporated.

{dbllseries(s1,s2,D,N1,N2,m,n,j,k1,k2,G)=
	G=sum(k2=1,abs(D)-1,k(D,k2)*q(D)^k2);
	v=vector(N1+2,x,newbern(x-2));
	vc=vector(N1+1,x,newbernchar(D,x-1));
	(sum(j=0,s1-1,(-1)^j*(s1-1)!/j!/(s1-j-1)!*
(sum(m=0,N1,
sum(n=-1,N1,vc[m+1]*v[n+2]*g(D)^(s1+s2+m+n)/(s2+j+n)/(s1+s2+m+n)))
	+sum(n=-1,N1,v[n+2]*g(D)^(s2+n+j)/(s2+n+j))*
	sum(j1=0,s1-j-1,(s1-j-1)!/j1!*g(D)^j1*
	sum(m=1,N2,k(D,m)*exp(-g(D)*m)/m^(s1-j-j1)))
	+sum(j1=0,s2+j-1,(s2+j-1)!/j1!*
	(sum(j2=0,s1-j-1,(s1-j-1)!/j2!*g(D)^(j1+j2)*
	sum(m=1,N2,k(D,m)*exp(-g(D)*m)/m^(s1-j-j2))*
	sum(n=1,N2,exp(-g(D)*n)/n^(s2+j-j1)))-
	sum(j3=0,s1-j+j1-1,(s1-j+j1-1)!/j3!*g(D)^j3*
	sum(k1=1,abs(D)-1,k(D,k1)*
	sumdouble(q(D)^(-k1),q(D)^k1*exp(-g(D)),s2+j-j1,s1-j+j1-j3,N2))/G))))))
/(s1-1)!/(s2-1)!}



\\tpllseries computes the triple l-value L(\chi,1,1;s1,s2,s3) for s3>1.  Bailey acceleration has been incorporated.

{tpllseries(s1,s2,s3,D,N1,N2,Gausssum)=
	Gausssum=sum(k2=1,abs(D)-1,k(D,k2)*q(D)^k2);	
	v=vector(N1+2,x,newbern(x-2));
	vc=vector(N1+1,x,newbernchar(D,x-1));
(sum(j=0,s1-1,sum(k1=0,s2-1,(-1)^(j+k1)*(s1-1)!/j!/(s1-j-1)!*(s2-1)!/k1!/(s2-k1-1)!*
(sum(m=0,N1,sum(n1=-1,N1,sum(n2=-1,N1,vc[m+1]*v[n1+2]*v[n2+2]*g(D)^(s1+s2+s3+m+n1+n2)/(s1+s2+s3+m+n1+n2)/(s2+s3+n1+n2+j)/(s3+n2+k1))))+
sum(n1=-1,N1,sum(n2=-1,N1,v[n1+2]*v[n2+2]*g(D)^(s2+s3+n1+n2+j)/(s2+s3+n1+n2+j)/(s3+n2+k1)))*sum(j1=0,s1-j-1,(s1-j-1)!/j1!*g(D)^j1*sum(m=1,N2,k(D,m)*exp(-g(D)*m)/m^(s1-j-j1)))+
sum(m=-1,N1,v[m+2]*g(D)^(s3+m+k1)/(s3+m+k1))*sum(j1=0,s2-k1+j-1,(s2-k1+j-1)!/j1!*(sum(j2=0,s1-j-1,(s1-j-1)!/j2!*g(D)^(j1+j2)*sum(n1=1,N2,k(D,n1)*exp(-g(D)*n1)/n1^(s1-j-j2))*sum(n2=1,N2,exp(-g(D)*n2)/n2^(s2-k1+j-j1)))-
sum(j3=0,s1-j+j1-1,(s1-j+j1-1)!/j3!*g(D)^j3*sum(k2=1,abs(D)-1,k(D,k2)*sumdouble(q(D)^(-k2),q(D)^k2*exp(-g(D)),s2-k1+j-j1,s1-j+j1-j3,N2))/Gausssum)))+
sum(j1=0,s3+k1-1,(s3+k1-1)!/j1!*(sum(j2=0,s2-k1+j-1,(s2-k1+j-1)!/j2!*(sum(j4=0,s1-j-1,(s1-j-1)!/j4!*g(D)^(j1+j2+j4)*sum(n1=1,N2,k(D,n1)*exp(-g(D)*n1)/n1^(s1-j-j4))*sum(n2=1,N2,exp(-g(D)*n2)/n2^(s2-k1+j-j2))*sum(n3=1,N2,exp(-g(D)*n3)/n3^(s3+k1-j1)))-
sum(j5=0,s1-j+j2-1,(s1-j+j2-1)!/j5!*g(D)^(j1+j5)*sum(k2=1,abs(D)-1,k(D,k2)*sumdouble(q(D)^(-k2),q(D)^k2*exp(-g(D)),s2-k1+j-j2,s1-j+j2-j5,N2))/Gausssum*sum(n3=1,N2,exp(-g(D)*n3)/n3^(s3+k1-j1)))))-
sum(j3=0,s2-k1+j+j1-1,(s2-k1+j+j1-1)!/j3!*(sum(j6=0,s1-j-1,(s1-j-1)!/j6!*g(D)^(j3+j6)*sum(n1=1,N2,k(D,n1)*exp(-g(D)*n1)/n1^(s1-j-j6))*sumdouble(1,exp(-g(D)),s3+k1-j1,s2-k1+j+j1-j3,N2))-
sum(j7=0,s1-j+j3-1,(s1-j+j3-1)!/j7!*g(D)^j7*sum(k2=1,abs(D)-1,k(D,k2)*sumtriple(1,q(D)^(-k2),q(D)^k2*exp(-g(D)),s3+k1-j1,s2-k1+j+j1-j3,s1-j+j3-j7,N2))/Gausssum)))))))))
/(s1-1)!/(s2-1)!/(s3-1)!}



\\Bailey acceleration has not been incorporated into the remaining series.


\\dzeta computes the double zeta value zeta(s1,s2) for s2>1.

{dzeta(s1,s2,N1,N2)=
       v=vector(N1+2,x,newbern(x-2));
(sum(k1=0,s1-1,(-1)^k1*(s1-1)!/k1!/(s1-k1-1)!*(sum(m=-1,N1,sum(n=-1,N1,v[m+2]*v[n+2]/(s2+k1+n)/(s1+s2+m+n)))+sum(m=1,N2,sum(n=-1,N1,v[n+2]*exp(-m)/(s2+k1+n)*sum(j=0,s1-k1-1,(s1-k1-1)!/j!/m^(s1-k1-j))))+sum(m=1,N2,sum(n=1,N2,exp(-m-n)*sum(j1=0,s2+k1-1,(s2+k1-1)!/j1!/n^(s2+k1-j1)*(sum(j2=0,s1-k1-1,(s1-k1-1)!/j2!/m^(s1-k1-j2))-sum(j3=0,s1-k1+j1-1,(s1-k1+j1-1)!/j3!/(m+n)^(s1-k1+j1-j3)))))))))/(s1-1)!/(s2-1)!}

\\dp11 computes the double polylog Li_{1,1}(z1*z2^{-1},z2) for z1, z2 roots of unity.  Here z2 has order E and z1*z2^{-1} has order D.

{dp11(z1,D,z2,E,N1,N2)=
       vz1=vector(N1+2,x,a(z1,D,x-2));
       vz2=vector(N1+1,x,a(z2,E,x-1));
       g1=min(g(D),g(E));
sum(m=0,N1,sum(n=-1,N1,vz2[m+1]*vz1[n+2]/(m+1)/(m+n+2)*g1^(m+n+2)))+sum(m=0,N1,vz2[m+1]*g1^(m+1)/(m+1))*sum(n=1,N2,z1^n*exp(-n*g1)/n)+sum(m=1,N2,sum(n=1,N2,z1^m*z2^n*exp(-(m+n)*g1)/m/(m+n)))}

\\li computes the value of the sth polylogarithm at a root of unity z of order D.

{li(z,D,s,N1,N2)=(sum(n=0,N1,a(z,D,n)*g(D)^(n+s)/(n+s))+sum(n=1,N2,z^n*sum(j=0,s-1,(s-1)!/j!*g(D)^j*exp(-g(D)*n)/n^(s-j))))/(s-1)!}

\\dps1 computes Li_{s,1}(z1*z2^{-1},z2) where z1 and z2 are roots of unity.  Here z2 has order E and z1*z2^{-1} has order D.

{dps1(z1,D,z2,E,s,N1,N2)=
        vz1=vector(N1+2,x,a(z1,D,x-2));
        vz2=vector(N1+1,x,a(z2,E,x-1));
        g1=min(g(D),g(E));
sum(k1=0,s-1,(-1)^k1*(s-1)!/k1!/(s-k1-1)!*(sum(m=-1,N1,sum(n=0,N1,vz1[m+2]*vz2[n+1]/(s+m+n+1)/(k1+n+1)*g1^(s+m+n+1)))+sum(m=1,N2,z1^m*sum(j=0,s-k1-1,(s-k1-1)!/j!*exp(-m*g1)*g1^j/m^(s-k1-j)))*sum(n=0,N1,vz2[n+1]/(k1+n+1)*g1^(k1+n+1))+sum(m=1,N2,sum(n=1,N2,z1^m*z2^n*exp(-(m+n)*g1)*sum(j1=0,k1,k1!/j1!/n^(k1+1-j1)*(sum(j2=0,s-k1-1,(s-k1-1)!/j2!*g1^(j1+j2)/m^(s-k1-j2))-sum(j3=0,s-k1+j1-1,(s-k1+j1-1)!/j3!*g1^j3/(m+n)^(s-k1+j1-j3))))))))/(s-1)!}

\\dp computes a double polylog Li_{s1,s2}(z1*z2^{-1},z2), where s2>1.  Here z2 has order E and z1*z2^{-1} has order D.

{dp(z1,D,z2,E,s1,s2,N1,N2)=
        vz1=vector(N1+2,x,a(z1,D,x-2));
        vz2=vector(N1+2,x,a(z2,E,x-2));
        g1=min(g(D),g(E));
sum(k1=0,s1-1,(-1)^k1*(s1-1)!/k1!/(s1-k1-1)!*(sum(m=-1,N1,sum(n=-1,N1,vz1[m+2]*vz2[n+2]*g1^(s1+s2+m+n)/(s1+s2+m+n)/(s2+k1+n)))+sum(n=-1,N1,vz2[n+2]*g1^(s2+k1+n)/(s2+k1+n))*sum(m=1,N2,z1^m*sum(j=0,s1-k1-1,(s1-k1-1)!/j!*g1^j*exp(-m*g1)/m^(s1-k1-j)))+sum(m=1,N2,sum(n=1,N2,z1^m*z2^n*exp(-(m+n)*g1)*sum(j1=0,s2+k1-1,(s2+k1-1)!/j1!/n^(s2+k1-j1)*(sum(j2=0,s1-k1-1,(s1-k1-1)!/j2!*g1^(j1+j2)/m^(s1-k1-j2))-sum(j3=0,s1-k1+j1-1,(s1-k1+j1-1)!/j3!*g1^j3/(m+n)^(s1-k1+j1-j3))))))))/(s1-1)!/(s2-1)!}