# The following algorithm takes a polynomial P in k[x] of degree 5 (where k is a specified field)
# whose roots are Qs= {q1,q2,q3,q4,q5}, qi = infinity if P has degree 4, and finds a list of all 
# (near) Belyi maps f in k(x) from our table whose (0,3,i)-exceptional points 
# (up to Mobius transformation)are Qs, i=3,4,6.


read"C:/Users/NiruViju/Desktop/DessinDatabase/Belyi_320":
read"C:/Users/NiruViju/Desktop/DessinDatabase/Belyi_420":
read"C:/Users/NiruViju/Desktop/DessinDatabase/Belyi_620":
read"C:/Users/NiruViju/Desktop/DessinDatabase/Belyi_one_320":
read"C:/Users/NiruViju/Desktop/DessinDatabase/Belyi_one_420":
read"C:/Users/NiruViju/Desktop/DessinDatabase/Belyi_one_620":
read"C:/Users/NiruViju/Desktop/FiveSings/MobiusTR5":
read"C:/Users/NiruViju/Desktop/I5":
read"C:/Users/NiruViju/Desktop/I5tilde":
read"C:/Users/NiruViju/Desktop/FindExceptionalPoints":


FindAllF:= proc(k::set,P)
local A, P1, P2, Inv1P, Inv2P, Inv1, Inv2, f,F1,i,i1,j,j1,k1,k2,Cand,Cand1,FinalCand, Mobs, TB, TB1, t_fact, t_val;   
 P1:= P;  Inv1P:=  I5(P1,x);   Inv2P:= I5tilde(P1,x);
 TB:= [Belyi_320, Belyi_420, Belyi_620]:
 TB1:= [Belyi_one_320, Belyi_one_420, Belyi_one_620]:

Cand:= {};

# Compute Candidate Belyi maps first:

for i to 3 do
  if i = 1 then A:= [3, 2, 0];
  elif i = 2 then A:= [4, 2, 0];
  else A:= [6, 2, 0];
  fi;

    for j in TB[i] do  
      if normal(j[2] - Inv1P) <> 0 then next; fi;
      P2:= FindExceptionalPoints(j[1],A);
      Cand:= Cand union {[j[1],P2,A]};
    od;

# Compute Candidate Belyi-1 maps now:

for j in TB1[i] do 
  P2:= FindExceptionalPoints(j[1],A);
    Inv1:= I5(P2,x);
    Inv2:= I5tilde(P2,x);
    t_fact:= factors(gcd(numer(Inv1 - Inv1P),
                     numer(Inv2 - Inv2P)),k)[2];

      for k1 in t_fact do 
        if degree(k1[1],t) > 1 then next; fi;
         t_val:= -evala(coeff(k1[1],t,0)/coeff(k1[1],t,1)); 
         #lprint(t_val);
        f:= factor(eval(j[1], t=t_val)); 
        if max(degree(numer(f),x),degree(denom(f),x)) <> 
           max(degree(numer(j[1]),x),degree(denom(j[1]),x))
             then next;
        fi;
        P2:= eval(P2, t=t_val);
        Cand:= Cand union {[f,P2,A]}; 
      od; 
 od;

# Now find Mobius transformation (if any) that carries the 
# roots of input polynomial to the exceptional points of 
# Candidates:

FinalCand:= {};
  for j in Cand do 
     Mobs:= MobiusTR5(j[2], P1);
     if Mobs = {} or type(Mobs,string) then next; fi;
       for j1 in Mobs do 
        F1:= factor(eval(j[1],x=j1));
         FinalCand:= FinalCand
                     union{[F1,j[3]]};
          
       od;
   od;

# Now compute Belyi-2 maps:

 
   Cand1:= Find_F4(P1,x,k)union Find_F6(P1,x,k);
   if Cand1 <> {} then
     for i1 in Cand1 do FinalCand:= FinalCand 
                       union {[i1,[3,2,0]]}; od;
   fi;
od;

{seq(`if`(P1 = FindExceptionalPoints(op(k2)),k2,NULL),k2=FinalCand)};


end;


#####################################################################


# This program computes degree 4 Belyi-2 map with five
# Exceptional points given in the input.
# The program sends exceptional point which comes from the
# linear factor in the input to infinity and 
# then calls FindF4.
# Branching pattern:  (1,3),(2,2),(1,1,1,1) above 0,1,infty 
# In this program, P is a polynomial of degree 5 
# (4 if it contains infinity) and B:= Base field.

Find_F4:= proc(P,x::name,B::set)     
local A1,A2,P1,P2,a,i,Cand;        
  # option trace;          
  Cand:= {};
  P1:= factors(P,B);
  A2:= {seq(`if`(degree(i[1],x) = 1, i[1], NULL),i=P1[2])};

if nops(A2) = 0 and degree(P,x) = 5 then return {} fi;

if degree(P,x) = 4 then Cand:= Cand union Find_F4a(P,x,B);
else 
 for i in A2 do 
      a:= -evala(coeff(i,x,0)/coeff(i,x,1));
      #solve(A2[1],x);
      P2:=factor(eval(eval(P,x= x+a),x=1/x)); 
   #P3:=factor(eval(P3,x=1/x));
  if degree(P2,x) = 4 then
     Cand:= Cand union factor(eval(Find_F4a(numer(P2)*x,x,B),x=1/(x-a)));
  else  Cand:= Cand union
    factor(eval(Find_F4a(numer(P2),x,B),x=1/(x-a)));
  fi;
 od;
fi; 
 Cand;
                        
end:


#####################################################################


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

Find_F4a := proc(P, x::name,B::set)  # The case when infty is a zero of f.
# option trace;
local sx, f,F,a,b,c,d,p0,p1,p2,i,k,sol,av,FB,EQ9,EQa,EQb1,eqns,res,b1v;

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 := evala(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;
	   sol := solve(eqns, {b0}); # we should just gcd them.
	   if sol = NULL then  next  fi;
	   F := eval(eval(FB, sol[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:


##############################################################

# This program computes a degree 6 Belyi minus 2 map which produces 5 
# exceptional points as in input.
# This program sends a linear factor to infinity and calls Find_F6a.
# Branching pattern:  (3,3),(2,2,2),(1,1,1,1,4).
# Here, P is a polynomial of degree 5 (4 if it contains infty).


Find_F6:= proc(P,x::name,B::set)   
local A1,P1,P2,P3,a,i,Res;
  
  P1:= factors(P,B);
  A1:= {seq(`if`(degree(i[1],x) = 1, i[1], NULL),i=P1[2])};

if nops(A1) = 0 and degree(P,x) = 5 then return {} fi;
           
  Res:= {};
  if degree(P,x) = 4 then  Res:= Res union Find_F6a(P,x,B);
  else 
   for i in A1 do  
        a:= -evala(coeff(i,x,0)/coeff(i,x,1));  #solve(i,x);
        P2:= factor(eval(eval(P,x= x + a),x=1/x));
        # P3:= factor(eval(P3,x=1/x));
         if degree(P,x) = 5 then
           Res:= Res union factor(eval(Find_F6a(numer(P2),x,B),
                                       x=1/(x- a)));
         else 
           Res:= Res union factor(eval(Find_F6a(numer(P2)*x,x,B),x=1/(x-
                 a)));
         fi;
     od;
     fi;

Res;
end:


#####################################################################

#This program computes Belyi-2 map from the given singularity
#structure of input differential operator assuming that the singularity
#with exponent difference int/3 is at infinity.


Find_F6a := proc(P, x::name, B::set)
local sx,f,F,a,b,c,d,p0,p1,p2,i,k,EQ12,res,av,EQd,dv,eqns,FB,sol;	
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;
	    sol := solve(eqns, {b}); # I should just gcd them.
	    if sol = NULL then next fi;
	    F := eval(eval(FB, sol[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: