\\ ---------------  GP code  ---------------------------------------
\\
\\ Time-stamp: <2001-01-04 14:25:12 cdoche>
\\
\\ Description: Routines for computing mahler measure.
\\
\\
\\ Original author:	  Christophe Doche
\\			  cdoche@math.u-bordeaux.fr
\\			  Université Bordeaux I
\\			  Laboratoire A2X	
\\
\\ Created:               Sun Apr  01  2001
\\


\\======================================================================
\\Given a polynomial, computes its mahler measure by Graeffe's method. nb is the number of iteration to be
\\choosen between 10 and 15

{polgraeffe(f,nb)=
	local(pol,M1,M2,nbtour,rf,sf,i,L):
	pol=f:
	M1=0:
	M2=1:
	nbtour=0:
	while(nbtour<nb,rf=0:
	      sf=0:
	      L=0:
      	      forstep(i=0,poldegree(pol),2,rf=rf+polcoeff(pol,i)*x^(i/2)):
	      forstep(i=1,poldegree(pol),2,sf=sf+polcoeff(pol,i)*x^((i-1)/2)):
              nbtour=nbtour++:
              pol=rf^2-x*sf^2:
              for(i=0,poldegree(pol),L=L+abs(polcoeff(pol,i))):
              M2=M1:
              M1=L^(1/2^nbtour)):
	      return(M1)}


\\======================================================================
\\Given a polynomial, computes its mahler measure

{polmahler(f)=
	local(M,racines,i):
	M=abs(polcoeff(f,poldegree(f))):
	racines=vecsort(abs(polroots(f))):
	for(i=1,length(racines),M=M*max(1,abs(racines[i]))):
	return(M)}


\\======================================================================
\\Given a polynomial, computes its mahler measure by Riemann sum. N is the number of step. 
\\Usually take N=10000. This routine is useful for sparse polynomials of large degree.

{unity(u)=
	exp(I*Pi*u)}

{polmahlerint(f,N)=
	local(M,racines,i,S):
	M=1:S=0:
	for(i=1,N,S+=log(abs(subst(f,variable(f),unity(i/N))))):
	return(exp(S/N))}



\\======================================================================
\\Given a reciprocal polynomial P, returns Q such that P(z)=Q(z+1/z)z^deg(Q)

{polcondens(f)=
	local(tmp,pol,i):
	tmp=f:
	pol=0:
	for(i=0,poldegree(f)/2,pol=pol+polcoeff(tmp,poldegree(f)-2*i)*x^(poldegree(f)/2-i):
		tmp=(tmp-polcoeff(tmp,poldegree(f)-2*i)*(x^2+1)^(poldegree(f)/2-i))/x):
	if(numerator(subst(pol,x,x+1/x))==f,return(pol))}


\\======================================================================
\\Given a reciprocal polynomial, computes its mahler measure

{polmahlerrecip(f)=
	local(tmp,racines,i):
	g=polcondens(f):
	M=abs(polcoeff(f,poldegree(f))):
	racines=polroots(g):
	for(i=1,matsize(racines)[1],tmp=polroots(x^2+1-racines[i]*x):
        	M=M*max(1,abs(tmp[1]))*max(1,abs(tmp[2]))):
	        return(M)}


\\======================================================================
\\Given a reciprocal polynomial, computes the number of roots outside the unit circle

{polnurecip(f)=
return(poldegree(f)/2-polsturm(polcondens(f))+polsturm(f)/2)}

\\======================================================================
\\Given a two variable polynomial, exchanges x and y

{polswap(p)=
	q=subst(subst(p,x,axa),y,bxb):
	return(subst(subst(q,bxb,x),axa,y))}


\\======================================================================
\\Given a two variable Laurent polynomial, returns a true polynomial

{pollaurent2pol(p)=
	return(p/x^valuation(p,x)/y^valuation(p,y))}


\\======================================================================
\\Given a two variable polynomial P, tests if P is reciprocal

{polisreciproc(p)=
	local(u,t,q):
	q=subst(subst(p,x,1/x),y,1/y):
	t=p/q:u=x^valuation(t,x)*y^valuation(t,y):
	if((t-u)*(t+u)==0,return(1),return(0))}




\\======================================================================
\\Given a two variable polynomial P, returns a polynomial Q with same mahler measure 
\\such that deg_y Q> deg_x Q and P(x,y)=Q(x^a,y^b) if any.

{polnormalize(p)=
	local(i,lx,ly,bmax,bmay,g,q,r):
	if(p==0,return(0),p=pollaurent2pol(p):
	lx=[]:
	for(i=0,poldegree(p,x),if(polcoeff(p,i,x)<>0,lx=concat(lx,[i]))):
	g=0:
	for(i=1,length(lx),g=gcd(g,lx[i])):
	if(g<>0,bmax=poldegree(p,x)/g,bmax=0):
	q=0:
	for(i=0,bmax,q=q+polcoeff(p,g*i,x)*x^(i)):
	ly=[]:
	for(i=0,poldegree(q,y),if(polcoeff(q,i,y)<>0,ly=concat(ly,[i]))):
	g=0:
	for(i=1,length(ly),g=gcd(g,ly[i])):
	r=0:
	if(g<>0,bmay=poldegree(p,y)/g,bmay=0):
	for(i=0,bmay,r=r+polcoeff(q,g*i,y)*y^(i)):
	p=r:
	if(poldegree(p,y)<poldegree(p,x),q=polswap(p),return(p)))}


{deg2(pol)=
	local(a,b,c,d):
	a=polcoeff(pol,2):
	b=polcoeff(pol,1):
	c=polcoeff(pol,0):
	if(a==0,if(b<>0,return([-c/b]),return([0])),d=sqrt(b^2-4*a*c):
	return([(-b-d)/2/a,(-b+d)/2/a]))}

{deg2f(pol)=
	local(a,b,c,d):
	a=polcoeff(pol,2):
	b=polcoeff(pol,1):
	c=polcoeff(pol,0):
	if(a==0,if(b<>0,return(log(max(1,abs(c/b)))),return(0)),d=sqrt(b^2-4*a*c):
	return(log(max(1,abs((-b-d)/2/a)))+log(max(1,abs((-b+d)/2/a)))))}

{deg3(pol)=
	local(a,b,c,p,q,v,U,V):
	pol=pol/polcoeff(pol,3):
	a=polcoeff(pol,2):
	b=polcoeff(pol,1):
	c=polcoeff(pol,0):
	p=-a^2/3+b:
	q=2*a^3/27-b*a/3+c:
	v=deg2(x^2+q*x-p^3/27):
	U=v[1]^(1/3):
	if(U<>0,V=-p/3/U,V=0):
	return([(U+V)-a/3,(U*exp(2*I*Pi/3)+V*exp(4*I*Pi/3))-a/3,(U*exp(4*I*Pi/3)+V*exp(2*I*Pi/3))-a/3])}

{deg3f(pol)=
	local(a,b,c,p,q,v,U,V):
	pol=pol/polcoeff(pol,3):
	a=polcoeff(pol,2):
	b=polcoeff(pol,1):
	c=polcoeff(pol,0):
	p=-a^2/3+b:
	q=2*a^3/27-b*a/3+c:
	v=deg2(x^2+q*x-p^3/27):
	U=v[1]^(1/3):
	if(U<>0,V=-p/3/U,V=0):
	return(log(max(1,abs((U+V)-a/3)))+log(max(1,abs((U*exp(2*I*Pi/3)+V*exp(4*I*Pi/3))-a/3)))+log(max(1,abs((U*exp(4*I*Pi/3)+V*exp(2*I*Pi/3))-a/3))))}


{deg4f(pol)=
	local(a,b,c,d,r,t,y,w,j):
	pol=pol/polcoeff(pol,4):
	a=polcoeff(pol,3):
	b=polcoeff(pol,2):
	c=polcoeff(pol,1):
	d=polcoeff(pol,0):
	r=deg3(x^3-b*x^2+(a*c-4*d)*x+d*(4*b-a^2)-c^2):
	t=[a^2/4-b+r[1],a^2/4-b+r[2],a^2/4-b+r[3]]:
	j=1:
	if(abs(t[2])>abs(t[1]),j=2):
	if(abs(t[3])>abs(t[j]),j=3):
	w=sqrt(t[j]):
	y=r[j]:
	if(w==0,return(4*log(max(1,abs(a/4)))),return(deg2f(x^2+(w+a/2)*x+a*y/4/w-c/2/w+y/2)+deg2f(-x^2+(w-a/2)*x+a*y/4/w-c/2/w-y/2)))}


\\======================================================================
\\Given a two variable polynomial, computes its mahler measure. N is the number of step. Take N=10000.
\\ If flag=0 the calculus should be a priori faster but less acurate

{polmahler2var(f,N,flag)=
	local(i,S,a,b,c,d,e,l):
	if(flag==0,resolution(pol)=if(poldegree(pol)>4,return(log(polgraeffe(pol,10))-log(abs(pollead(pol)))),
			   	   if(poldegree(pol)==4,deg4f(pol),
				   if(poldegree(pol)==3,deg3f(pol),
				   if(poldegree(pol)<=2,deg2f(pol))))),
		   resolution(pol)=if(poldegree(pol)>4,return(log(polmahler(pol))-log(abs(pollead(pol)))),
				   if(poldegree(pol)==4,deg4f(pol),
				   if(poldegree(pol)==3,deg3f(pol),
				   if(poldegree(pol)<=2,deg2f(pol)))))):
	S=0:
	l=pollead(f,x):
	for(i=1,N,
	g=subst(f,y,unity(i/N)):
	S+=resolution(g)):
	return(polmahler(l)*exp(S/N))}



{helpmahler()=
print("-------------------------------------------------------------------");
print("      		Computational number Theory                        "); 
print("-------------------------------------------------------------------");
print("");
print("Description: Routines for computing mahler measure of polynomials");
print("in 1 and 2  variables");
print("");
print("Original Author:     Christophe Doche");
print("                     cdoche@math.u-bordeaux.fr");
print("                     Université Bordeaux I");
print("                     Laboratoire A2X");
print("");
print("-------------------------------------------------------------------");
print("");
print("Help functions");
print("");
print("helppolgraeffe       helppolmahler      helppolmahlerint      helppolcondens");
print("helppolmahlerrecip   helppolnurecip     helppollaurent2pol    helppolnormalize");
print("helppolisreciprocal  helppolmahler2var");
print("");}

{helppolgraeffe()=
print("");
print("polgraeffe(p,nb) Given a polynomial p, computes its mahler measure by Graeffe's
method. nb is the number of iteration to be choosen between 10 and 15");
print("");}

{helppolmahler()=
print("");
print("Given a polynomial, computes its mahler measure");
print("");}

{helppolmahlerint()=
print("");
print("polgraeffe(p,N) Given a polynomial p, computes its mahler measure by Riemann sum. 
N is the number of step. Usually take N=10000. 
This routine is useful for sparse polynomials of large degree");
print("");}

{helppolcondens()=
print("");
print("Given a reciprocal polynomial P, returns Q such that P(z)=Q(z+1/z)z^deg(Q)");
print("");}

{helppolmahlerrecip()=
print("");
print("Given a reciprocal polynomial, computes its mahler measure");
print("");}

{helppolnurecip()=
print("");
print("Given a reciprocal polynomial, returns the number of roots outside the unit circle");
print("");}

{helppollaurent2pol()=
print("");
print("Given a Laurent polynomial in two variables, returns a true polynomial");
print("");}

{helppolnormalize()=
print("");
print("Given a two variable polynomial P, returns a polynomial Q with same 
mahler measure such that deg_y Q> deg_x Q and P(x,y)=Q(x^a,y^b) if any.");
print("");}

{helppolisreciprocal()=
print("");
print("Given a polynomial in two variables P, return 1 if P is reciprocal, 0 otherwise");
print("");}

{helppolmahler2var()=
print("");
print("polmahler2var(p,N,flag) Given a polynomial in two variables p, computes 
its mahler measure. N is the number of step. Take N=10000. 
If flag=0 the calculus should be a priori faster but less acurate");
print("");}
print("");
print("Type helpmahler for getting the list of help functions of")
print("the file mahler.gp")
print("")