\\ This code is designed to compute Belyi Lattes maps.

\\ let's first increase precision

\p 1000

\\ the following function takes a vector specifying generators of the lattice of an elliptic curve, another vector and an integer expected to be 3,4 or 6 and returns the polynomial that admits elements in the vector as its roots
{
polynomial_from_roots(n,v) = 
my(w,pol,normalization);
pol = 1;
if(n == 3,w = ellperiods([1,exp(Pi*I/3)]);print("using p prime"); normalization=ellwp(w,1/3*(1 + exp(2*Pi*I/6)),1)[2]; for(i = 1,length(v), pol = pol*(x - ( (normalization - ellwp(w,v[i],1)[2])/(2*normalization) ) ););,if(n == 6, w = ellperiods([1,exp(Pi*I/3)]); print("using p cube"); for(i = 1,length(v),pol = pol*(x - (ellwp(w,v[i])^3 / ellwp(w,1/2)^3)););,if( n == 4, w = ellperiods([1,I]); print("using p square");for(i = 1,length(v),pol = pol*(x - (ellwp(w,v[i])/ellwp(w,1/2))^2););, print("unknown n");break; );  ););
pol
}

\\ the following two are functions to calculate the set over which the ramification points are calculated. In the second function, we need the extra vecextract command as zero is a pole of \wp'.

{
calculate_vector_of_zeros_3 (n) = 
my(v,d);
v = [];
d = (1+exp(2*Pi*I/6));
for(i = 0,n-1, for(j = 0,n-1, v = concat(v, 1/n* ( i*1 + j*exp(2*Pi*I/6) +  1/3*d  )  );  );  );
v
}


{
calculate_vector_of_poles_3 (n) = 
my(v);
v = [];
for(i = 0,n-1, for(j = 0,n-1, v = concat(v, 1/n* ( i*1 + j*exp(2*Pi*I/6) )  );  );  );
vecextract( v, concat("2..",Str(n^2)))
}

{
	calculate_vector_of_zeros_4 (n) = 
	my(v,d);
	v = [];
	d = (1+I)/2;
	for(i = 0,n-1, for(j = 0,n-1, v = concat(v, 1/n* ( i*1 + j*I +  d  )  );  );  );
	v
}

{
	calculate_vector_of_poles_4 (n) = 
	my(v);
	v = [];
	for(i = 0,n-1, for(j = 0,n-1, v = concat(v, 1/n* ( i*1 + j*I )  );  );  );
	vecextract( v, concat("2..",Str(n^2)))
}

{
	calculate_vector_of_zeros_6 (n) = 
	my(v,d);
	v = [];
	d = (1+exp(2*Pi*I/6))/3;
	for(i = 0,n-1, for(j = 0,n-1, v = concat(v, 1/n* ( i*1 + j*exp(2*Pi*I/6) +  d  )  );  );  );
	\\for(i = 0,n-1, for(j = 0,n-1, v = concat(v, 1/n* ( i*1 + j*exp(2*Pi*I/6) +  2*d  )  );  );  );
	v
}

{
	calculate_vector_of_poles_6 (n) = 
	my(v);
	v = [];
	for(i = 0,n-1, for(j = 0,n-1, v = concat(v, 1/n* ( i*1 + j*exp(Pi*I/3) )  );  );  );
	vecextract( v, concat("2..",Str(n^2)))
}

\\ the following function produces belyi lattes maps 

{
	belyi_lattes(s,n) = 
	if(s == 3, purify_belyi(polynomial_from_roots(3,calculate_vector_of_zeros_3(n)))/purify_belyi(polynomial_from_roots(3,calculate_vector_of_poles_3(n))),if(s == 4, purify_belyi(polynomial_from_roots(4,calculate_vector_of_zeros_4(n)))/purify_belyi(polynomial_from_roots(4,calculate_vector_of_poles_4(n))), if(s==6, purify_belyi(polynomial_from_roots(6,calculate_vector_of_zeros_6(n)))/purify_belyi(polynomial_from_roots(6,calculate_vector_of_poles_6(n))), print("unknown symmetry!") )  )  )
}


\\ the following function converts a given float to a rational number. In the end, we chose to use the existing bestappr() function.
{
convert(x)=
my(n=0);
while(x-floor(x*10^n)/10^n!=precision(0., 400),n++);
floor(x*10^n)/10^n
}

\\ for(i = 1,length(coef),if(imag(coef[i])>10^(-120), print("the "i"th coefficient has non-zero imaginary part!");break;, coef[i] = convert(real(coef[i])*precision(1., 200)) ) ); this is to check whether there is a value with non-zero imaginary part.

{
purify_belyi(pol) = 
my(p = 0,coef = []);
for(i = 0,poldegree(pol), coef = concat(coef,polcoef(pol,i)));
for(i = 1,length(coef), if(abs(real(coef[i]))<10^(-200), coef[i] = 0, coef[i] = bestappr(real(coef[i])); ); );
for(i=1,length(coef),p = p + coef[i]*x^(i-1););
p
}

{
output(n) = 
my(zeros = [],poles = []);
for(i = 3,n,  poles = concat(poles,purify_belyi(polynomial_from_roots(3,calculate_vector_of_poles_3(i)))); zeros = concat(zeros,purify_belyi(polynomial_from_roots(3,calculate_vector_of_zeros_3(i)))));
for(i = 1,n-2,printtex(zeros[i]/poles[i]);)
}

{
print_table3 (n,m) = 
my(a , b  ) ;
for(i = n,m, a = factor( purify_belyi(polynomial_from_roots(3,calculate_vector_of_zeros_3(i))) ); b = factor( purify_belyi(polynomial_from_roots(3,calculate_vector_of_poles_3(i))) ); writetex(table3,"\\phi_{3," ); writetex(table3,i ); writetex( table3,"}(z) &= "  ); writetex(table3,"\\frac{");  for(j = 1,matsize(a)[1], writetex(table3,"(" ); writetex(table3,a[j,1]); writetex(table3,")^{"); writetex(table3,a[j,2]); writetex(table3,"}");  ); writetex(table3, "}{"); for(j = 1,matsize(b)[1],writetex(table3,"("); writetex(table3,b[j,1]); writetex(table3,")^{"); writetex(table3,b[j,2]); writetex(table3,"}");  ); writetex(table3,"}\\\\");  );
}

{
	print_table4 (n,m) = 
	my(a , b  ) ;
	for(i = n,m, a = factor( purify_belyi(polynomial_from_roots(4,calculate_vector_of_zeros_4(i))) ); b = factor( purify_belyi(polynomial_from_roots(4,calculate_vector_of_poles_4(i))) ); writetex(table4,"\\phi_{4," ); writetex(table4,i ); writetex( table4,"}(z) &= "  ); writetex(table4,"\\frac{");  for(j = 1,matsize(a)[1], writetex(table4,"(" ); writetex(table4,a[j,1]); writetex(table4,")^{"); writetex(table4,a[j,2]); writetex(table4,"}");  ); writetex(table4, "}{"); for(j = 1,matsize(b)[1],writetex(table4,"("); writetex(table4,b[j,1]); writetex(table4,")^{"); writetex(table4,b[j,2]); writetex(table4,"}");  ); writetex(table4,"}\\\\");  );
}

{
	print_table6 (n,m) = 
	my(a , b  ) ;
	for(i = n,m, a = factor( purify_belyi(polynomial_from_roots(6,calculate_vector_of_zeros_6(i))) ); b = factor( purify_belyi(polynomial_from_roots(6,calculate_vector_of_poles_6(i))) ); writetex(table6,"\\phi_{6," ); writetex(table6,i ); writetex( table6,"}(z) &= "  ); writetex(table6,"\\frac{");  for(j = 1,matsize(a)[1], writetex(table6,"(" ); writetex(table6,a[j,1]); writetex(table6,")^{"); writetex(table6,a[j,2]); writetex(table6,"}");  ); writetex(table6, "}{"); for(j = 1,matsize(b)[1],writetex(table6,"("); writetex(table6,b[j,1]); writetex(table6,")^{"); writetex(table6,b[j,2]); writetex(table6,"}");  ); writetex(table6,"}\\\\");  );
}

{
print_to_file_3 (n) = 
my(zeros = [],poles = [], str = "");
for(i = 3,n,  poles = concat(poles,purify_belyi(polynomial_from_roots(3,calculate_vector_of_poles_3(i)))); zeros = concat(zeros,purify_belyi(polynomial_from_roots(3,calculate_vector_of_zeros_3(i)))));
for(i = 1,n-2, writetex(belyi31,i+2); writetex(belyi31," & "); writetex(belyi31,zeros[i]/poles[i]); writetex(belyi31," & "); writetex(belyi31,factor(zeros[i])); writetex(belyi31," & ");  writetex(belyi31,factor(poles[i])); writetex(belyi31,"\\ \\"); )
}


{
periodic_rational(word) = 
my(check = 0,string = "",n=2 );
word = Vec(Str(word));
while(check == 0, for(i = 1,n,if(word[n+i]==word[i],check =1;,check = 0;n +=1; break;)));
for(i=1,n,string = concat(string,word[i]););
eval(string)/(10^n-1)
}

{
get_rational_coeff_3(p) = 
my(polynomial);
for(i = 0, poldegree(p), polynomial = polynomial + (bestappr(real(polcoef(p,i)))+ bestappr(imag(polcoef(p,i))/sqrt(3)*I))*x^i);
polynomial
}

{
get_rational_coeff_1(p) = 
my(polynomial);
for(i = 0, poldegree(p), polynomial = polynomial + (bestappr(real(polcoef(p,i)))+ bestappr(imag(polcoef(p,i))*I))*x^i);
polynomial
}


{
get_rational_coeff_for_factorization(p) = 
my(polynomial);
for(i = 0, poldegree(p), polynomial = polynomial + (bestappr(real(polcoef(p,i)))+ bestappr(imag(polcoef(p,i))/sqrt(3)*I*abc))*x^i);
polynomial
}