# This program computes a Belyi minus 2 map which produces 5 singularities from 0,1,infty.
# This program takes "singularity structure" for a Belyi minus two map of degree 6 and returns the map. All singularities come from exp diff 0. This program
# sends each linear factor to infinity and calls Find_F6a.
# Branching pattern:  (3,3),(2,2,2),(1,1,1,1,4)  above 0,1,infty with exp diff. 1/3,1/2,0 respectively.

Find_F6:= proc(P,x::name,B::set)   # Here, P is a set of lists [monic irred poly in Q, exp. diff.]. 
   local A1,P1,P2,P3,B1,a,i,Res,RD;
    # Following writes the field extension in terms of RootOf only: 
       RD:= radfield(indets(P,{radical,nonreal}));
       P1:= eval(P,RD[1]);
       A1:= {};
          # Each exp. diff. must be zero mod Z.
            for i in P1 do if type(i[2],integer) then A1 := A1 union {i[1]};  else return {} fi od;
                if add(deg(i,x),i=A1) <> 5 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;                 
                            Res:= {};
                              for i in A1 do 
                                 if deg(i,x) <> 1 then next;
                                 elif i = infinity then  Res:= Res union factor(eval(Find_F6a(P2,x,B1),RD[2]));
                                 else 
                                    a:= solve(i,x);
                                    P3:= factor(eval(P2,x= x + a));
                                    P3:= factor(eval(P3,x=1/x));
                                      if deg(P2,x) = 5 then
                                         Res:= Res union factor(eval(Find_F6a(numer(P3),x,B1),{x=1/(x-a),op(RD[2])}));
                                      else 
                                         Res:= Res union factor(eval(Find_F6a(numer(P3)*x,x,B1),{x=1/(x-a),op(RD[2])})); fi; fi; od;
    Res;
end:


Find_F6a := proc(P, x::name, B::set)
	local sx, f, a,b,c,d,p0,p1,p2, i, EQ12, res, EQ3, j, av, EQd, k,k1, dv, 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);
         # Following is the equation in terms of T = a, where a comes from FB. 
	EQ12 := 1048576*T^12+524288*p2*T^10+131072*p1*T^9+73728*T^8*p2^2-294912*T^8*p0+49152*p2*p1*T^7-21504*p1^2*T^6-18432*T^5*p0*p1+4608*T^5*p2^2*p1-6912*T^4*p0^2-432*T^4*p2^4-1920*T^4*p1^2*p2+3456*T^4*p0*p2^2-736*p1^3*T^3-288*T^2*p0*p1^2+72*T^2*p1^2*p2^2+16*p1^3*p2*T+p1^4;
	EQ12 := evala(Factors( evala(Primpart(EQ12,T)), B));
	res := NULL;
	     for i in EQ12[2] do if degree(i[1],T)=1 then
			av := evala( -coeff(i[1],T,0)/coeff(i[1],T,1) );
		        EQd := 48*a^2*d^2-48*a^2*(8*a^2+p2)*d+512*a^6+160*p2*a^4-40*p1*a^3+12*a^2*p2^2-4*p1*a*p2-p1^2;
			EQd := evala(Factors( evala(Primpart(eval(EQd,a=av),d)), B));
			   for k in EQd[2] do if degree(k[1],d)=1 then
				dv := evala( -coeff(k[1],d,0)/coeff(k[1],d,1) );
				eqns := {-3*p0*d+d^3+2*p0*b-c^2-3*p0*a^2, 6*a*d^2-3*p1*a^2+2*p1*b-3*p1*d-2*b*c, -3*p2*a^2+3*d^2+2*p2*b-3*p2*d+12*a^2*d-6*a*c-b^2, 12*a*d+8*a^3-6*a*b-2*c};
				# FB := (x^3+3*a*x^2+b*x+c)^2/(x^2+2*a*x+d)^3;
				eqns := factor(eval(subs(c = a*(6*d+4*a^2-3*b), eqns), {a=av,d=dv}));
				FB := (x^3+3*a*x^2+b*x+a*(6*d+4*a^2-3*b))^2/(x^2+2*a*x+d)^3;
				so := solve(eqns, {b}); # I should just gcd them.
				if so = NULL then
					next
				fi;
				F := eval(eval(FB, so[1]), {a = av, d=dv});  # either we do gcd or may have b's with multiplicity.
				F := traperror( evala(F) );
				 if F<>lasterror and degree(numer(F),x)=6 then
                                    # Now we need to adjust 0,1,infty.. Branching type: [3,3],[2,2,2],[1,1,1,1,4].
					res := res, sort(factor(1/(1-F)),x);
				 fi
		        fi od	
	    fi od;
    {res};
end: