# This program computes degree 4 Belyi minus 2 map (which produces 5 singularities from 0,1,infty) from given "Singularity structure". The program sends the root of poly corresponding to exp. diff. 1/3 to infinity
# and calls FindF4.
# Branching pattern:  (1,3),(2,2),(1,1,1,1) above 0,1,infty  with exp diff. 1/3,1/2,0 respectively.

Find_F4:= proc(P,x::name,B::set)    # Here, P is a set of lists [monic irred poly in Q, exp. diff.]. 
  # option trace;          
     local A1,A2,P1,P2,P3,B1,a,i,RD;
       A1:= {};  A2:= {};
     # Following writes the field extension in terms of RootOf only: 
       RD:= radfield(indets(P,{radical,nonreal}));
       P1:= eval(P,RD[1]);
          for i in P1 do if type(i[2],integer) then A1 := A1 union {i[1]}; elif denom(i[2]) = 3 and type(numer(i[2]),integer) then A2 := A2 union {i[1]}; else return {} fi od;
          # We need four 0's and one 1/3 mod Z.
          if nops(A2) <> 1 or deg(A2[1],x) <> 1 or add(deg(i,x),i = A1) <> 4 then return {} fi;
             P2:= 1;
                for i in P1 do if i[1] <> infinity then P2:= sort(expand(P2*i[1]),x) fi; od;
                     if nargs=2 then  B1:= indets(P2, {RootOf}) fi;                 
                          if A2[1] = infinity then  return factor(eval(Find_F4a(P2,x,B1),RD[2]));
                          else
                             a:= solve(A2[1],x);
                             P3:= factor(eval(P2,x= x + a));
                             P3:= factor(eval(P3,x=1/x));
                                if degree(P2,x) = 4 then    return  factor(eval(Find_F4a(numer(P3)*x,x,B1),{x=1/(x-a),op(RD[2])}));
                                else  return  factor(eval(Find_F4a(numer(P3),x,B1),{x=1/(x-a),op(RD[2])}));
                                fi; 
                         fi;                        
 end:

# This program computes degree 4 Belyi minus 2 map for a given "singularity structure" assuming that the poly with exp. diff. 1/3 is at infinity.

Find_F4a := proc(P, x::name, B::set)         # The case when infty falls at zero of f.
   # option trace;
	local sx, f, a,b,c,d,p0,p1,p2, i, EQ9, res, EQa, j, b1v, EQb1, k, av, eqns, FB, so, F;
	if degree(P,x) <> 4 then
		return {}
	elif coeff(P,x,4) <> 1 then
		return procname( collect(P/lcoeff(P,x), x, evala), x,B )
	elif coeff(P,x,3) <> 0 then
		sx := coeff(P,x,3)/4;
		f := procname(collect(subs(x = x-sx, P), x, evala), x,B);
		return factor(subs(x=x+sx, f))
	fi;
	p0,p1,p2 := seq(coeff(P,x,i),i=0..2);

       # This is the equation in T=b1 we find using elimination method.
	EQ9 := T^9+24*p2*T^7-168*p1*T^6-78*p2^2*T^5+1080*p0*T^5+336*p1*p2*T^4+80*p2^3*T^3+1728*p0*p2*T^3-636*p1^2*T^3-168*p1*p2^2*T^2-864*p0*p1*T^2-27*p2^4*T-432*p0^2*T+216*p2^2*p0*T-120*p2*p1^2*T-8*p1^3;
	EQ9 := factors( evala(Primpart(EQ9,T)), B);
	res := NULL;
	    for i in EQ9[2] do if degree(i[1],T)=1 then
			b1v := evala( -coeff(i[1],T,0)/coeff(i[1],T,1) );			
	                EQb1 := -p1+b1*p2-b1^3-6*b1^2*a-6*b1*a^2;
			EQa := evala(Factors( evala(Primpart(eval(EQb1,b1=b1v),a)), B));
			  for k in EQa[2] do if degree(k[1],a)=1 then
				av := evala( -coeff(k[1],a,0)/coeff(k[1],a,1) );
				eqns := {b0^2-p0+2*b1*a^3, -p1+2*b0*b1-6*b1*a^2, 2*b0-p2+b1^2+6*b1*a};
				# FB := 2*b1*(x-a)^3/(x^2+b1*x+b0)^2;
				eqns := factor(eval(eqns, {a=av,b1=b1v}));
				FB := 2*b1*(x-a)^3/(x^2+b1*x+b0)^2;
				so := solve(eqns, {b0}); # I should just gcd them.
				if so = NULL then  next  fi;
				F := eval(eval(FB, so[1]), {a = av, b1=b1v});
				F := traperror( evala(F) );
				if F<>lasterror and degree(denom(F),x)=4 then
                                 # Now we need to adjust 0,1,infty.. Branching type: [1,3],[2,2],[1,1,1,1].
					res := res,sort(factor(F/(F-1)),x);
				fi
			fi od;
	    fi od;
	{res}
end:


deg:= proc(P,x)
 local P1;
   P1 := P;  
     if P1 = infinity then return 1;  else return degree(P1,x); fi;
end: