# This program finds a set of mobius transformations on the sets of 5 
# points:
# To find the mobius transformation in degree 5 polynomials (5 point
# case),we have the following cases:
# n:=degree of irreducible factors.

# Case 1: 5 <----> 5.
# Case 2: 4,1 <-----> 4,1.
# Case 3: 3,2 <-----> 3,2.
# Case 4: 3,1,1 <-----> 3,1,1.
# Case 5: 2,2,1 <-----> 2,2,1.
# Case 6: 2,1,1,1 <-----> 2,1,1,1.
# Case 7: 1,1,1,1,1 <-----> 1,1,1,1,1.


MobiusTR5:= proc(f,g)   # f,g are the polynomials of degree 5 (degree 4 when
                        # they have infinity).
local  i,f1,g1,S1,S2,RD,deg1,deg2,deg3,Fctf,Fctg,Base_Field;

Base_Field:= indets([f,g],{radical,nonreal,RootOf});
RD:= radfield(Base_Field);

#Now rewrite everything in terms of RootOf (if any):
Base_Field:= eval(Base_Field,RD[1]);
f1:= eval(f,RD[1]);  g1:= eval(g,RD[1]);

Fctf:=  factors(f1,Base_Field);    Fctg:=  factors(g1,Base_Field);
   if not {degree(f1,x),degree(g1,x)} subset {4,5} then
     error "Wrong input.";
   fi;
S1:= {seq(i[1],i=Fctg[2])};   S2:= {seq(i[1],i=Fctf[2])};
deg1:={};  deg2:={};  deg3:= {};
  for i in S1 do
     if degree(i,x)=1 then deg1 := deg1 union {i}; 
     elif  degree(i,x)=2 then deg2 := deg2 union {i};
     else deg3 := deg3 union {i}; 
     fi;
  od;
           if  nops(deg3) <> 0 then  eval(Case_a(f1,g1),RD[2]);
           elif  nops(deg2) = 2  then  eval(Case_b(f1,g1),RD[2]);
           elif  nops(deg2) = 1  then  eval(Case_c(f1,g1),RD[2]);
           else  eval(Case_d(f1,g1),RD[2]);
           fi;
end:


# Case_a covers the following  Cases 1-4:

# Case 1:   5  <------> 5.
# Case 2:   4,1  <------>  4,1.
# Case 3:   3,2  <------>  3,2.
# Case 4:   3,1,1  <------>  3,1,1.

Case_a:=proc(f,g)    # f,g polys. of degree 5.
local m,ms,p,q,RTs,eqns,Sol,res,M,G,S,S1,i,j,g1,Q;
# options trace;
G:= factors(g,indets(g,RootOf)); S:= {seq(i[1],i=G[2])}; g1:={};
  for i in S do if degree(i,x) > 2 then  g1:= g1 union {i}; else next; fi;
  od;
  if nops(g1) <> 1 then  return "Not in Case_a"; fi;
 
 G:=op(g1);  m:= (a*x+b)/(c*x+d);   p:= RootOf(G);
 RTs:= evala(Roots(f,p)); Q:= {seq(i[1],i=RTs)}; q:={};
   for i in Q do if has(i,p) then  q:= q union {i};  else next;  fi;  od;

# Now we need (a*p+b)/(c*p+d)=q.
 Sol:={};
     for i in q do eqns:= {coeffs(IsZero(m,p,i),p)};
                   Sol:= Sol union {solve(eqns)};
     od;
 S1:={};
    for i in Sol do
      ms:= traperror(evala(eval((a*x+b)/(c*x+d),i)));
        if ms <> lasterror and has(ms,x)  and  IsMobiusTr(f,g,ms) then  
           S1:= S1  union  {ms}; 
        fi;
           
    od; 
 S1;
end:


# Case 5: 2,2,1 <-----> 2,2,1.

Case_b:=proc(f,g)    # f,g polys. of degree 5.
local m,p,p1,q,q1,eqn1,eqn2,eqns,RTs,F,G,S,S1,S2,i,j,k,f1,f2,g1,g2,Q,Q1,Sol,M,ms;
# options trace;
G:= factors(g,indets(g,RootOf));  F:= factors(f,indets(f,RootOf));
S1:= {seq(i[1],i=G[2])}; S2:= {seq(i[1],i=F[2])};
g1:={}; g2:={};   f1:= {};  f2:= {};
  
      for i in S1 do
         if degree(i,x) = 2 then  g1:= g1 union {i};
         elif degree(i,x) = 1 then g2:= g2 union {i};
         else  return "Not in case_b.";
         fi; 
      od;
        
      for i in S2 do
          if degree(i,x)=2 then f1:= f1 union {i};
          elif degree(i,x)=1 then  f2:= f2 union {i};
          else return "Not in case_b.";
          fi;
      od;

      if nops(g1) <> 2  or nops(f1) <> 2 then  return "Not in Case_b"; fi;

 S:={}; p:= RootOf(g1[1]);  RTs:= evala(Roots(f,p));
 Q:= {seq(j[1],j=RTs)};  Q1:={};
  for j in Q do if has(j,p) then  Q1:= Q1 union {j} ; else next; fi; od;        
  if  nops(g2)=1 then  
     p1:= -evala(coeff(g2[1],x,0)/coeff(g2[1],x,1));  #solve(g2[1]);
  else p1:= infinity ;
  fi;

  if  nops(f2)=1 then  q1:= -evala(coeff(f2[1],x,0)/coeff(f2[1],x,1));
  else q1:= infinity ;
  fi;
 m:= (a*x+b)/(c*x+d);
  Sol := {};
  eqn1:= {IsZero(m,p1,q1)};
    for k in Q1 do  eqn2:={coeffs(IsZero(m,p,k),p)}; 
                    eqns:= eqn1  union  eqn2 ; 
                    Sol:= Sol  union {solve(eqns)}; 
    od;

 M:={};
   for j in Sol do
     ms := traperror(evala(Normal( eval(m, j),expanded)));
       if ms<>lasterror and has(ms,x)  and  IsMobiusTr(f,g,ms) then  
           M := M union {ms}; 
       fi;
   od; 
M;
end:

#Case 6: 2,1,1,1 <-----> 2,1,1,1.

Case_c:= proc(f,g)    # f,g polys. of degree 5.
local f1,f2,g1,g2,g3,i,i1,j,k,m,p,p1,p2,q,q1,F,G,M,Q,Q1,S1,S2,RTs,
eqn1,eqn2,Eqns,Sol,ms;
# options trace;
G:= factors(g,indets(g,RootOf));  F:= factors(f,indets(f,RootOf));
S1:= {seq(i[1],i=G[2])}; S2:= {seq(i[1],i=F[2])};
g1:={};  g2:={};  f1:={};  f2:= {}; 
       for i in S1 do
         if degree(i,x) = 2 then  g1:= g1 union {i};
         elif degree(i,x) = 1 then g2:= g2 union {i};
         else return "Not in Case_c";
         fi; 
       od;
 
       for i in S2 do 
           if degree(i,x)=2 then  f1:= f1 union {i};
           elif degree(i,x)=1 then f2:= f2 union {i};
           else return "Not in Case_c";
           fi;
       od;
       if nops(g1) <> 1  or nops(f1) <> 1 then return "Not in Case_c"; fi;

   p1:= {seq(-evala(coeff(i,x,0)/coeff(i,x,1)),i=g2)};
   q1:= {seq(-evala(coeff(i,x,0)/coeff(i,x,1)),i=f2)}; 
        if nops(p1) <> 3  then  p1:= p1 union {infinity} ;  fi;
        if nops(q1) <> 3  then  q1:= q1 union {infinity} ;  fi;
  p:= RootOf(op(g1));  RTs:= evala(Roots(op(f1),p));
  Q:= {seq(j[1],j=RTs)};   Q1:= select(has,Q,p);  
  if Q1={} then return {}; fi;
  m:= (a*x+b)/(c*x+d);  M:={};   Sol:= {};

for i in Q1 do
   eqn1:= {coeffs(IsZero(m,p,i),p)};   p2:= p1[1];   
      for j in q1 do eqn2:= {IsZero(m,p2,j)}; 
                     Eqns:= eqn1  union  eqn2 ; 
                     Sol:= Sol  union  {solve(Eqns)}; 
      od;
od;
 
      for i1 in Sol do
            ms := traperror(evala(Normal(eval(m, i1), expanded)));
            if ms<>lasterror and has(ms,x) and  IsMobiusTr(f,g,ms) then
               M := M union {ms};
            fi;  
      od; 
M;
end:


#Case 7: 1,1,1,1,1 <-----> 1,1,1,1,1.

Case_d:= proc(f,g)    # f,g polys. of degree 5.
local  F,G,S1,S2,f1,g1,i,i1,j,j1,j2,j3,p1,p2,p3,m,M,P,Q,eqn1,eqn2,eqn3,eqns,sol,ms;
# options trace;
  G:= factors(g,indets(g,RootOf));  F:= factors(f,indets(f,RootOf));
  S1:= {seq(i[1],i=G[2])}; S2:= {seq(i[1],i=F[2])};
  g1:={};  f1:={}; 
    for i in S1 do 
       if degree(i,x) = 1 then g1:= g1 union {i}; 
       else return "Not in Case_d";
       fi; 
    od; 
    for i in S2 do 
       if degree(i,x)=1 then  f1:= f1 union {i}; 
       else return "Not in Case_d";  
       fi; 
    od;

   P:= {seq(-evala(coeff(i,x,0)/coeff(i,x,1)),i=g1)};
   Q:= {seq(-evala(coeff(i,x,0)/coeff(i,x,1)),i=f1)};
     if nops(P) <> 5 then P:= P union {infinity} ;  fi;
     if nops(Q) <> 5 then Q:= Q union {infinity} ;  fi;

m:= (a*x+b)/(c*x+d);

# The Mobius transformation carries elements of P to elements of Q. 
# In this case we can find all such transformations using 3 pts. in P and 
# evaluating their all possible images via m. 

  sol:={};
  p1:=P[1];  p2:=P[2];  p3:=P[3];
    for j1 in Q do  eqn1 := {IsZero(m,p1,j1)};
      for j2 in Q minus {j1} do  eqn2 := {IsZero(m,p2,j2)}; 
          for j3 in Q minus {j1,j2} do
            eqn3 := {IsZero(m,p3,j3)};   
            eqns:=eqn1 union eqn2 union eqn3;  
            sol:=  sol  union {solve(eqns)}; od; od; od;        
        M:={};    
          for i1 in sol do
             ms := traperror(evala(Normal(eval(m, i1), expanded)));
             if ms<>lasterror and has(ms,x) and  IsMobiusTr(f,g,ms)  then 
                M := M union {ms}; 
             fi;  
          od;  
  M;
end:


# This program checks,for degree 5 polynomials,if they are Mobius 
# equivalent with the given Mobius Transformation or not.

IsMobiusTr:=proc(f,g,m)  # f,g are degree 5 polys and m is mobius 
                         #  transformation.
local f1,g1,IsZero;
# options trace;
    # if has([f,g],infinity) then  return 
    # procname(op(eval([f,g], infinity = x-1)),m); fi;
    g1:=g;
    f1:= numer(evala(eval(f/(x-dummy)^5,x=m)));
    if degree(g1,x) < degree(f1,x) then g1:= g1*denom(m); fi;
    IsZero:= evala(Expand(f1- coeff(f1,x,degree(g))/lcoeff(g)*g));
    evalb(IsZero = 0);
end:


# This program evaluates the mobius transformation at a given point and 
# produces a linear equation in terms of a,b,c,d:
# m:= (a*x+b)/(c*x+d)
# Let eval(m,x=p) = q .
# Then  evala(eval(m,x=p)-q)=0.

IsZero := proc(m,p,q)    # eval(m,x=p)=q.
    if p = infinity then       # when m evaluated at infinity.
      if q = infinity then  c ;    #  c = 0.
      else  evala(a-q*c) ;   #  a/c = q. 
      fi;
    else 
      if q = infinity  then evala(c*p+d); 
      else  evala((a*p+b)-(c*p+d)*q);
      fi;
    fi;
end: