# This program finds 2F1 type solutions: 
# y(x) = exp(int(r,x))*(r_0*S(f) + r_1*diff(S(f),x)), of a second order 
# linear differential operator L with 5 non removable singularities which 
# lies in Class(H^(1/12,5/12)_(1,x)) where r,r_0,r_1 are rationalfunctions 
# and S(f)=2F1(a,b;c|f). 
# The tables for Belyi and Belyi_one are also included as separate files in # this folder. So they can be uploaded from the folder or can be read from 
# the path given:


read"Belyi_320":
read"Belyi_420":
read"Belyi_620":
read"Belyi_one_320":
read"Belyi_one_420":
read"Belyi_one_620":
read"MobiusTR5":


_Envdiffopdomain:= [Dx,x]:
with(DEtools):


Solver5:= proc(L,x,B::set)   # B:= Base Field.
local  L1,sln1,sln2,sln3,sings,RD,n,Base_Field;

if nargs = 2 then 
  Base_Field:= indets([args], {RootOf, radical, nonreal}); 
  return procname(args,Base_Field);
else Base_Field:= B;
fi;

  RD:= radfield(indets(Base_Field));
 # Rewrite L and Base_Field in terms of RootOf (if any) only:
  L1:= eval(L,RD[1]);
  Base_Field:= eval(Base_Field,RD[1]);
  sings:= eval(printsing(L1),infinity = 1);
  n:= add(nPoints(i[1],x),i = sings);
    if n <> 5 then return "We only deal with 5 singularities"; fi;

# Call Solver for [e0,e1,ei] = [0,1/2,1/6] first:
 sln1:= Solver5_02k(L1,x,6,Base_Field);
   if sln1 <> NULL then return {eval(sln1,RD[2])}; fi;

# Call Solver for [e0,e1,ei] = [0,1/2,1/4] second:
 sln2:= Solver5_02k(L1,x,4, Base_Field);
   if sln2 <> NULL then return {eval(sln2,RD[2])}; fi;

# Call Solver for [e0,e1,ei] = [0,1/2,1/3] finally:
 sln3:= Solver5_02k(L1,x,3, Base_Field);
   if sln3 <> NULL then return {eval(sln3,RD[2])};
   else return {}; 
   fi;

end:


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


# This program computes a solution (if any) of a differential operator with # 5 regular singularities in terms of 2F1 with exp diff 0,1/2,1/k, k=3,4,6 
# (or its decomposition) at 0,1,infty resp.

Solver5_02k:= proc(L,x,n::integer,C::set)
#n:= 3,4,6, C:: Base field. 

local a,b,c,i,j,k,k1,A,B,G,LH,LH1,f_dec,f_dec1,g_dec,g_dec1,f_final,
soln,Candidates,ExpDiff;
#option trace;

if n=3 then A:=Belyi_320: B:=Belyi_one_320:
            ExpDiff:= [1/3,1/2,0];
elif n=4 then A:=Belyi_420: B:=Belyi_one_420:
            ExpDiff:= [1/4,1/2,0];
elif n=6 then A:=Belyi_620: B:=Belyi_one_620: 
            ExpDiff:= [1/6,1/2,0];
else return "Wrong input";
fi;

f_final:= {}; soln:= NULL;
LH:= Dx^2-(-2*x+x*e0+x*e1+1-e0)*Dx/(x*(x-1))+(1/4)*(e0+e1-ei-1)*(e0+e1+ei-
     1)/(x*(x-1));
# Compute candidate maps first
Candidates:= FindCandidates_k20(L,x,A,B,ExpDiff,C);
# Now check if they have decompositions:

if A = Belyi_320 then 
   g_dec:= [[-4*x*(x-1),[1/3,1/3,0]]];
elif A = Belyi_420 then 
   g_dec:= [[-4*x*(x-1),[1/4,1/4,0]],[(1/4)*x^2/(x-1),[1/2,0,0]]];
   g_dec1:= [[(2*x-1)^2,[0,0,0]],[(1/4)*(2*x-1)^2/((x-1)*x),[0,0,0]]];
elif A = Belyi_620 then 
   g_dec:= [[-4*x*(x-1),[1/6,1/6,0]],[(1/4)*x^2/(x-1),[1/3,0,0]]];
fi;

 for i in Candidates do
   for j in g_dec do 
     f_dec:= DecompDeg2(op(i),op(j),C);
     if A <> Belyi_420 then f_final:= f_final union f_dec;
     else  
#We have to check second round of decomposition:
      for k in f_dec do 
       if k[2][1] <> 1/2 then f_final:= f_final union {k};
       else 
         for k1 in g_dec1 do 
            f_dec1:= DecompDeg2(op(k),op(k1),C);
            f_final:= f_final union f_dec1;
         od;
       fi;
      od;
     fi;
   od;
 od;


for i in f_final do
  a := -i[2][3]/2-i[2][2]/2-i[2][1]/2+1/2; 
  b := -i[2][2]/2-i[2][3]/2+i[2][1]/2+1/2; c := 1-i[2][3];
  LH1:= collect(transfo(eval(LH,{e0=i[2][3],e1=i[2][2],ei=i[2][1]}),
                        1/i[1]),Dx,factor);
  G := equiv(LH1, L);
    if G = 0 then next; fi;
    G:=eval(diffop2de(G,y(x)),y(x)=hypergeom([a,b],[c],factor(1/i[1])));
      if not(has(indets(G, function),hypergeom)) then next; fi;
	 G := collect(G, select(has,indets(G, function),hypergeom));
      soln := soln,G;
od;
if soln <> NULL then return sort([op({soln})], length)[1]; fi;
  
end:


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


# This program computes the candidate Belyi and near Belyi maps f such that
# H^(a,b)_(c,x) and input differential operator have the same singularity
# structure.

FindCandidates_k20:= proc(L,x,A::list,B::list,C::list,E::set)  
#k:= 3,4,6, A:=Belyi_k20,B:=Belyi_one_k20,C:=[1/k,1/2,0],E:: base field 

local i,i1,j,k,f,Fs,fs,L1,B_2,Sing_L,Sing_L2,Sing_L3,Inv1_L,Inv2_L,
Inv1_L2,Inv2_L2,t_val,t_fact,Mobs,expdiff_L,MinPoly_Inv1_L;

L1:=L;
Sing_L:= eval(printsing(L1),infinity=1);
Inv1_L:= I5(mul(i[1],i=Sing_L),x);
MinPoly_Inv1_L:= sqrfree(evala(Norm(x-Inv1_L)))[2][1][1];
Inv2_L:= I5tilde(mul(i[1],i=Sing_L),x);
expdiff_L:= sort([seq(`if`(i[1]=1,i[2],i[2]$degree(i[1],x)),i=Sing_L)]);

Fs:= {};
# compute candidate Belyi maps:
for i in A do
  if MinPoly_Inv1_L <> i[3] then next fi;
  if not has(Inv1_L, RootOf) then
     Fs:= Fs union {[i[1], C]};
  elif has(Inv1_L, RootOf) then
     i1:= [op(i), op(indets(i[2],RootOf)),
          evala(Norm(x-op(indets(i[2],RootOf))))];
     # i1:= [Belyi map, I5, minpoly of I5, RootOf(...), minpoly of RootOf]
        for j in Roots(i1[-1], indets(L1,RootOf)) do
          if evala(eval(i1[2], i1[4]= j[1])-Inv1_L) = 0 then
             #lprint(i[1]);
             Fs:= Fs union {[eval(i1[1], i1[4] = j[1]),C]};
          fi;
        od;
  fi;
od;

# compute candidate Belyi-1 maps:
for i in B do 
  if EqualModZ(expdiff_L, i[2]) = 1 then next; fi; 
    Sing_L2:= eval(DetermineSings(i[1],x,C,E),infinity=1);
    Inv1_L2:= I5(mul(i[1],i=Sing_L2),x);
    Inv2_L2:= I5tilde(mul(i[1],i=Sing_L2),x);
    t_fact:= factors(evala(Gcd(numer(Inv1_L2 - Inv1_L),
                     numer(Inv2_L2 - Inv2_L))),E)[2];

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


# compute Mobius transformation to match singularities from these maps to 
# the singularities of L:

fs:= {};
for i in Fs do 
  Sing_L3:= eval(DetermineSings(i[1],x,i[2],E),infinity=1);
   if add(nPoints(j[1],x),j = Sing_L3) <> 5 then next; fi;
     Mobs:= MobiusTR5(mul(j[1],j=Sing_L3), mul(j[1],j=Sing_L),E);
   if Mobs = {} or type(Mobs,string) then next; fi;
   for k in Mobs do fs:= fs union {[factor(eval(i[1],x=k)),i[2]]}; od;
od;

# Compute Belyi-2 maps for the case [E0,E1,Ei]=[1/3,1/2,0]:

if A = Belyi_320 then 
   B_2:= Find_F4(Sing_L,x,E) union Find_F6(Sing_L,x,E);
   if B_2 <> {} then
     for i in B_2 do fs:= fs union {[i,C]}; od;
   fi;
fi;

fs;

end:

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


# The following program computes the number of points (singularities)
# in a polynomial in x:

nPoints:= proc(P,x)
  if P = 1 then 1 else degree(P,x); fi;
end:


#############################################################
 
      
# This program finds the decomposition of degree two(if any)
# f = f_dec(gs) where f_dec with degree two transforms the 2F1 
# with exponent differences E1 to another 2F1 with exponent differences
# E2. Find more details in the paper (Check Fig 1): 

DecompDeg2:= proc(f,E1::list,g,E2::list,A::set) # f:=g(h), where g=f_dec.
local h,t,i,gs;

 if irem(max(degree(numer(f),x),degree(denom(f),x)),2) <> 0
   # non-decomposable 
   then return {[f,E1]};
 fi;
  gs:= NULL;
  h:= evala(Factors(numer(eval(g,x=t)-f),A))[2];
 
 for i in h do
  if degree(i[1],t) = 1 then  
    gs:= gs, [-evala(coeff(i[1],t,0)/coeff(i[1],t,1)),E2];
  fi;
 od;
 if gs = NULL then return {[f,E1]}; else return {gs}; fi;
end:


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


#This program computes the singularity structure of a differential #operator:

printsing:=proc(L)
      option remember;
	local S, sing, i,j,true_sing,v;
	S := lcoeff(primpart(L,Dx),Dx);
	S := gcd(S, lcoeff(primpart(LCLM(L, Dx) ,Dx),Dx) );
	sing := factors(S)[2];
        sing:=[seq(i[1], i=sing), infinity];
        true_sing:=NULL;
        for i in sing do
		j := `if`(i=infinity, i, RootOf(i,x));
           v := gen_exp(L,T,x = j );
           if nops(v)=2 or v[1][1]=v[1][2] or 
              (not type(v[1][1]-v[1][2], integer)) or 
              formal_sol(L,`has logarithm?`,x=j) then
                 true_sing := true_sing, [i, `if`(nops(v)=1, 
                              v[1][2]-v[1][1], v[2][1]-v[1][1]) ]
           fi;
        od;
	{true_sing}
end:


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


#This program applies change of variable (x --> f) to a differential #operator:

transfo:=proc(L,a)
#option trace;
local i,f;         
f:= add(mult(subs(x=a,coeff(L,Dx,i)),(1/diff(a,x)*Dx)$i),i=0..degree(L,Dx));        
sort(collect(f/lcoeff(f,Dx),Dx,factor),Dx);
end:


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


# This program finds the singularities after applying the change of variable # x-->F. The program assumes that the base GHE has exponent differences 
# 1/k,1/2,0 (k <> -1,0,1) at 0,1,infinity respectively.

DetermineSings := proc(F,x,S::list,B:: set) # S: exponent diff of GHE
                                            # B: Base field.
        local V,V1,F1,d,d1,n,n1,W,S1,i;
        S1:= S;
        d1:= denom(S1[1]);
        if nops(S1)<>3 or S1[3]<>0 or denom(S1[2]) <> 2 or 
           not member(d1,{3,4,6}) then 
           return "Wrong list of exponent differences";
        fi;
        V := factors(F,B)[2];  #  [ [poly, n], ...
                             #  n>0 means root,  n<0 means pole.
        n := degree(denom(F),x) - degree(numer(F),x);
        if n<>0 then
           V := [ [x-infinity, n], op(V)]    # infinity --> x-infinity.
        fi:
        W := factors(numer(evala(1-F)),B)[2];
        F1:= factor(evala(1-F));
        n1:= degree(denom(F1),x) - degree(numer(F1),x);
          if n = 0 and n1 <> 0 then  W:= [[x-infinity,n1], op(W)]; fi;
        V1 := {};
          for i in V do 
            if has(i[1],x) and i[2]>0 and irem(i[2],d1)<>0 then
               V1:= V1 union {[i[1],1/d1*i[2]]};
            elif has(i[1],x) and i[2]<0 then 
               V1:= V1 union {[i[1],0]};
            fi;
          od;
             for i in W do
               if has(i[1],x) and irem(i[2],2)<>0 then
                  V1:= V1 union {[i[1],1/2]}; 
               fi;
             od;
     eval(V1,infinity=x-infinity);
end:


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

# The following program checks if any two lists of 5 
# entries (exponent differences of five singularities) are equal mod Z. 
# Returns 0 if they are equal, returns 1 otherwise.

EqualModZ:= proc(l1,l2)
# l1,l2 are lists of 5 numbers.
local i,L1,L2;
  L1:= sort([seq(min(frac(abs(i)),1-frac(abs(i))),i=l1)]); 
  L2:= sort([seq(min(frac(abs(i)),1-frac(abs(i))),i=l2)]);
  if nops(L1) = 5 and L1 = L2 then return 0; else return 1; fi;
 end:


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


# The following program computes a five point invariant (taking the sum
# of j-invariants of any combination of four points).
# Input: a square-free polynomial P of degree 4 or 5, representing a 5 point # subset S
# of C union {infinity}.   Here infinity is in S iff P has degree 4.
#
# Output: An element of C that is invariant under Mobius transformations 
# applied to S.
#
# Note: The input must be square-free, otherwise the program causes a 
# division by 0.

I5 := proc(P,x)
	local c0,c1,c2,c3,i;
	if lcoeff(P,x)<>1 then
		return procname(collect(P/lcoeff(P,x),x,evala),x)
		# From now on P will be monic
	elif degree(P,x)<>5 then
		# Handle the case that the degree is not 5:
		return I4_and_infinity(P,x)
		# From now on P is monic of degree 5.
	elif coeff(P,x,4) <> 0 then
		return procname(collect(eval(P,x=x-coeff(P,x,4)/5),x,evala),x);
		# After this, P = x^5 + c3*x^3 + c2*x^2 + c1*x + c0 for some c.i
	fi;
	c0,c1,c2,c3 := seq(coeff(P,x,i), i=0..3);

	# The denominator of the following expression is discrim(P,x)
	# which is non-zero iff P is square-free.
	#
	# See the end of the file for an explanation how this formula was 
     # found.

	evala(
256*(-180*c1^2*c2^4+c3^4*c2^4+15000*c2^2*c1*c0^2-14000*c0*c1^3*c2+5400*c3^2*c1^2*c2*
c0-555*c3^4*c2*c0*c1-5700*c3*c2^3*c1*c0+945*c0*c2^5+260*c3^4*c1^3-720*c3^2*c1^4+6750
*c3^2*c2^2*c0^2-120*c3^3*c2^2*c1^2+115*c3^3*c2^3*c0-6000*c3^3*c1*c0^2-25000*c3*c2*c0
^3+960*c3*c2^2*c1^3+30000*c3*c1^2*c0^2+125000*c0^4-6*c3^5*c1*c2^2+945*c3^5*c0^2+2240
*c1^5+9*c3^6*c1^2+15*c3^2*c2^4*c1)/(108*c3^5*c0^2-72*c3^4*c2*c0*c1+16*c3^4*c1^3-900*
c3^3*c1*c0^2+16*c3^3*c2^3*c0-4*c3^3*c2^2*c1^2+825*c3^2*c2^2*c0^2+560*c3^2*c1^2*c2*c0
-128*c3^2*c1^4-3750*c3*c2*c0^3+2000*c3*c1^2*c0^2-630*c3*c2^3*c1*c0+144*c3*c2^2*c1^3+
3125*c0^4+2250*c2^2*c1*c0^2+108*c0*c2^5-1600*c0*c1^3*c2-27*c1^2*c2^4+256*c1^5)
	)
end:

# The 5 point invariant, with 4 points given by a monic square-free
# polynomial P, and the 5'th point is infinity.
#

I4_and_infinity := proc(P,x)
	local c0,c1,c2,c3,i;
	if degree(P,x)<>4 or lcoeff(P,x)<>1 then
		error "Input must be monic of degree 4"
	fi;
	c0,c1,c2,c3 := seq(coeff(P,x,i), i=0..3);

	# Remark: The formula used here was found by the following 
     # computation:
	#
	# unassign('c0','c1','c2','c3');
	# P4 := x^4 + c3*x^3 + c2*x^2 + c1*x + c0;
	# R := RootOf(P4,x);
	# P3 := evala(Normal( P4/(x-R) ));
	# F := algcurves[j_invariant](y^2 - P3, x,y);
	# F := evala(Trace(F));  # this is the sum of the j's of the 4 subsets 
     #of S that contain infinity
	# F := factor( F + algcurves[j_invariant](y^2 - P4, x,y) ); 
     # this added j( S minus {infinity} )

	evala(
-256*(260*c0*c2^4-720*c0^2*c2^2+37*c3^3*c1*c2^3-8*c1^2*c2^3-6*c2^5*c3^2+c3^4*c2^4-6*c3^5*c1*c2^2-45*c3^5*c1*c0-
57*c2^4*c1*c3+72*c1^3*c2*c3-120*c3^2*c2^3*c0-120*c3^2*c1^2*c0+960*c3^2*c2*c0^2-1680*c1*c3*c0^2+420*c3^3*c2*c1*c0
-1080*c3*c2^2*c1*c0+9*c2^6+2240*c0^3-216*c1^4+1440*c1^2*c2*c0+9*c1^2*c3^6-8*c1^3*c3^3-180*c3^4*c0^2+90*c1^2*c3^2
*c2^2-57*c1^2*c3^4*c2+15*c3^4*c0*c2^2)/(-256*c0^3+4*c3^2*c2^3*c0+6*c3^2*c1^2*c0-144*c1^2*c2*c0+4*c1^2*c2^3-144*
c3^2*c2*c0^2-18*c1^3*c2*c3+27*c3^4*c0^2+192*c1*c3*c0^2+128*c0^2*c2^2-18*c3^3*c2*c1*c0+80*c3*c2^2*c1*c0-16*c0*c2^
4+27*c1^4+4*c1^3*c3^3-c1^2*c3^2*c2^2)
	)
end:

# Remark: the formula in I5 was computed as follows:
#
# P5 := x^5 + c3*x^3 + c2*x^2 + c1*x + c0;
# c0 := solve( subs(x=p, P5), c0);
# P5 := collect(subs(x = x+p, P5),x,normal);
# P5 := factor(subs(x=1/x, P5));
# P4 := numer(P5);
# P4 := collect(P4/lcoeff(P4,x), x, factor);
# INV := I4_and_infinity(P4,x);
# INV := factor(evala(subs(p = RootOf(c0-c00,p), INV)));
# unassign('c0');
# subs(c00 = c0, INV);


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


# This program computes another five point invariant (taking the product
# of j-invariants of any combination of four points)
# Input: a square-free polynomial P of degree 4 or 5, representing a 5 point # subset S
# of C union {infinity}.   Here infinity is in S iff P has degree 4.
#
# Output: An element of C that is invariant under Mobius transformations 
# applied to S.
#
# Note: The input must be square-free, otherwise the program causes a 
# division by 0.

I5tilde := proc(P,x)
	local c0,c1,c2,c3,i;
	if lcoeff(P,x)<>1 then
		return procname(collect(P/lcoeff(P,x),x,evala),x)
		# From now on P will be monic
	elif degree(P,x)<>5 then
		# Handle the case that the degree is not 5:
		return I4_and_infinity_tilde(P,x)
		# From now on P is monic of degree 5.
	elif coeff(P,x,4) <> 0 then
		return procname(collect(eval(P,x=x-coeff(P,x,4)/5),x,evala),x);
		# After this, P = x^5 + c3*x^3 + c2*x^2 + c1*x + c0 for some c.i
	fi;
	c0,c1,c2,c3 := seq(coeff(P,x,i), i=0..3);

	# The denominator of the following expression is discrim(P,x)
	# which is non-zero iff P is square-free.
	#
	# See the end of the file for an explanation how this formula was 
      # found.

	evala(
(-6*c3^5*c2^2*c1+1200*c0^2*c3^3*c1-660*c2^3*c3*c1*c0+920*c3^2*c1^2*c2*c0+21*c3^4*c1*c2*c0-1200*c1^3*c2*c0-13*c2^3*c3^3*c0-88*c2^2*c3^3*c1^2-192*c2^2*c3*c1^3-3000*c0^2*c2^2*c1+150*c0^2*c2^2*c3^2+5000*c3*c0^
3*c2+14000*c3*c0^2*c1^2+81*c2^5*c0+36*c2^4*c1^2+304*c3^2*c1^4+132*c3^4*c1^3+81*c0^2*c3^5+c3^4*c2^4+9*c3^6*c1^2+192*c1^5+100000*c0^4+15*c3^2*c2^4*c1)/(-900*c0^2*c3^3*c1-630*c2^3*c3*c1*c0+560*c3^2*c1^2*c2*c0
-72*c3^4*c1*c2*c0-1600*c1^3*c2*c0+16*c2^3*c3^3*c0-4*c2^2*c3^3*c1^2+144*c2^2*c3*c1^3+2250*c0^2*c2^2*c1+825*c0^2*c2^2*c3^2-3750*c3*c0^3*c2+2000*c3*c0^2*c1^2+108*c2^5*c0-27*c2^4*c1^2-128*c3^2*c1^4+16*c3^4*c1^
3+108*c0^2*c3^5+256*c1^5+3125*c0^4)
	)
end:
# This program is not exactly the product of the j invariants, but rather,
# the cube root of this product, scaled
# by a suitable constant.


# The 5 point invariant, with 4 points given by a monic square-free 
# polynomial P, and the 5'th point is infinity.
#

I4_and_infinity_tilde := proc(P,x)
	local c0,c1,c2,c3,i;
	if degree(P,x)<>4 or lcoeff(P,x)<>1 then
		error "Input must be monic of degree 4"
	fi;
	c0,c1,c2,c3 := seq(coeff(P,x,i), i=0..3);

	# Remark: The formula used here was found by the following 
     # computation:
	#
	# unassign('c0','c1','c2','c3');
	# P4 := x^4 + c3*x^3 + c2*x^2 + c1*x + c0;
	# R := RootOf(P4,x);
	# P3 := evala(Normal( P4/(x-R) ));
	# F := algcurves[j_invariant](y^2 - P3, x,y);
	# F := evala(Norm(F));  # this is the product of the j's of the 4 
     # subsets of S that contain infinity
	# F := factor( F * algcurves[j_invariant](y^2 - P4, x,y) ); 
     # multiplied by j( S minus {infinity} )
	#
	# We then removed a constant factor from F and took the cube-root, 
     # which lead to:

	evala(

(-3*c3^5*c1+c3^4*c2^2+3*c3^4*c0+19*c3^3*c2*c1-6*c3^2*c2^3-16*c3^2*c2*c0-8*c3^2*c1^2-30*c3*c2^2*c1+9*c2^4-8*c1*c3*c0+24*c0*c2^2+24*c1^2*c2+16*c0^2)*(c2^2+12*c0-3*c1*c3)/(256*c0^3-4*c3^2*c2
^3*c0+144*c3^2*c2*c0^2+144*c1^2*c2*c0-6*c3^2*c1^2*c0-192*c1*c3*c0^2+18*c3^3*c2*c0*c1+18*c1^3*c2*c3+c1^2*c2^2*c3^2-27*c3^4*c0^2-128*c0^2*c2^2+16*c0*c2^4-27*c1^4-80*c0*c2^2*c3*c1-4*c3^3*c1^3-4*c1^2*c2^3)

	)
end:

# Remark: the formula in I5 was computed as follows:
#
# P5 := x^5 + c3*x^3 + c2*x^2 + c1*x + c0;
# c0 := solve( subs(x=p, P5), c0);
# P5 := collect(subs(x = x+p, P5),x,normal);
# P5 := factor(subs(x=1/x, P5));
# P4 := numer(P5);
# P4 := collect(P4/lcoeff(P4,x), x, factor);
# INV := I4_and_infinity_tilde(P4,x);
# INV := factor(evala(subs(p = RootOf(c0-c00,p), INV)));
# unassign('c0');
# subs(c00 = c0, INV);
# lprint(%);

# Test:
# f := randpoly(x,degree=5, dense);
# g := numer(normal(subs( x = (3*x+5)/(7*x-11), f)));
# if I5tilde(f,x) <> I5tilde(g,x) then
#	error "program does not work!"
# fi;


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


# This program computes degree 4 Belyi-2 map (which produces 5 
# singularities from 0,1,infty) from given "Singularity structure".
# The program sends the singularity with exp. diff. 1/3 to infinity and 
# then 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.]
                                    # B:= Base field. 
local A1,A2,P1,P2,P3,a,i;        
  # option trace;          
 
  P1:= P;  A1:= {};  A2:= {};
  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 nPoints(A2[1],x) <> 1 or add(nPoints(i,x),i = A1) <> 4 
   then return {} 
fi;
P2:= 1;
for i in P1 do if i[1] <> 1 then P2:= sort(expand(P2*i[1]),x) fi; od;                 
if A2[1] = 1 then  return Find_F4a(P2,x,B);
else a:= -evala(coeff(A2[1],x,0)/coeff(A2[1],x,1)); #solve(A2[1],x);
  P3:=factor(eval(eval(P2,x= x+a),x=1/x)); #P3:=factor(eval(P3,x=1/x));
  if degree(P2,x) = 4 then
     return  factor(eval(Find_F4a(numer(P3)*x,x,B),x=1/(x-a)));
  else  return  factor(eval(Find_F4a(numer(P3),x,B),x=1/(x-a)));
  fi; 
fi;
                        
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 
# singularities from 0,1,infty.
# This program takes "singularity structure" of a Belyi-2 map of 
# degree 6 and returns the map. All singularities come from exp diff 0.
# This program sends a 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,a,i,Res;
  P1:= P;  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(nPoints(i,x),i=A1) <> 5 then return {}; fi;
           
  P2:= 1;
  for i in P1 do if i[1] <> 1 then P2:= sort(expand(P2*i[1]),x) fi; od;                 
  Res:= {};
  for i in A1 do 
    if nPoints(i,x) <> 1 then next;
    elif i = 1 then  Res:= Res union Find_F6a(P2,x,B);
    else 
        a:= -evala(coeff(i,x,0)/coeff(i,x,1));  #solve(i,x);
        P3:= factor(eval(eval(P2,x= x + a),x=1/x));
        # P3:= factor(eval(P3,x=1/x));
         if nPoints(P2,x) = 5 then
           Res:= Res union factor(eval(Find_F6a(numer(P3),x,B),
                                       x=1/(x- a)));
         else 
           Res:= Res union factor(eval(Find_F6a(numer(P3)*x,x,B),x=1/(x-
                 a)));
         fi;
     fi;
   od;

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:


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

# This program computes the projective equivalence:  y ----> exp(int(r,x))# (r_0 y + r_1 diff(y,x)) between two operators:
# Example to show the syntax:
#
# _Envdiffopdomain := [Dx,x]:  # (see help page of DEtools[mult] for this)
#
# L1 := x^2*Dx^2-4*x^6+4*x^4-3/4;
# L2 := Dx^2 + 2-10*x+4*x^2-4*x^4;
# equiv(L1,L2);
# --> Gives the map from the solution space of L1 to the solution space
#  of L2, showing effectively that solving L1 is equivalent to solving L2.
#     To get the inverse map, call equiv(L2,L1).

# Checks if M2 can be solved in terms of M1, and if so, finds the map.

equiv:=proc(M1,M2,d)
	local DF,x,T1,L1,T2,L2,C0,D0,S,C1,D1,dd,v,v1;
  if nargs>2 then
	if type(d,list(name)) and nops(d)=2 then
	  _Envdiffopdomain:=d
	elif type(d,function) then
	  _Envdiffopdomain:=[DF,x];
	  return `DEtools/diffop2de`(procname(seq(
		  `DEtools/de2diffop`(args[i],d),i=1..2)),d)
	else
		error "wrong number or type of arguments"
	fi
  elif not assigned(_Envdiffopdomain) then
	error "domain is not assigned"
  fi;
	DF,x := op(_Envdiffopdomain);
	if not type([M1,M2],list(polynom(ratpoly(anything,x),DF))) then
		error "wrong arguments or dependent variable not specified"
	elif degree(M1,DF)<>2 or degree(M2,DF)<>2 then
		error "case not handled yet"
	fi;

        if indets([args],{RootOf,radical,nonreal}) = {} then
                Normalizer := normal
        else
                Normalizer := evala
        fi;

	T1:=Normalizer(coeff(M1,DF,1)/coeff(M1,DF,2)/2);
	L1:=`DEtools/symmetric_product`(M1,DF-T1);
	T2:=Normalizer(coeff(M2,DF,1)/coeff(M2,DF,2)/2);
	L2:=`DEtools/symmetric_product`(M2,DF-T2);
	C0:=Normalizer(coeff(L1,DF,0)/coeff(L1,DF,2));
	D0:=Normalizer(coeff(L2,DF,0)/coeff(L2,DF,2)-C0);
  if D0=0 then
	S:=1
  else
	C1:=Normalizer(diff(C0,x));
	D1:=Normalizer(diff(D0,x));
	dd:=Normalizer(D1/D0);
	v:=map(Normalizer,[2*diff(C1,x)+diff(D1,x)-2*dd*C1-D1*dd+D0^2,
	 6*C1+D1-4*dd*C0, 2*(D0+2*C0), -dd, 1]);
	v:=`DEtools/Expsols`(v,0,x,`use Int`,`no algext`,radical,denom,2);
	if nops(v)=0 then
		return 0
	elif nops(v)>1 then
		userinfo(1,'dsolve',"input is reducible");
		return 0
	else
		v:=v[1];
		v1:=Normalizer(diff(v,x));
		S:=collect(v*DF +
		  (diff(v1,x,x)+(4*C0+5*D0)*v1+v*(2*C1+D1))/2/D0-2*v1
		,[exp,DF],Normalizer)
	fi
  fi;
	S:=combine(collect(exp(Int(-T2,x))*`DEtools/mult`(
	 S,exp(Int(T1,x))),[exp,DF],Normalizer),exp);
	v:=seq(`if`(op(0,i)=exp,i,NULL),i=indets(S,function));
	if nops([v])<>1 then
		error "exp not collected"
	fi;
	S := subs(v=`DEtools/kovacicsols/e_int`(Normalizer(diff(op(v),x)),x),S);
    if has(S,RootOf) and not has(S,DF) then add(`DEtools/kovacicsols/combi`(normal(coeff(S,DF,i)),x)*DF^i, i=0..degree(S,DF));
    else  collect(subs(v=`DEtools/kovacicsols/e_int`(Normalizer(diff(op(v),x)),x),S),DF,Normalizer);
    fi;	
end: