with(DEtools):

# Given an operator L,then candiuv(SetofMob(singularity(L))) will give all candidates of [u,v].


# first step is to find singularities of the given operator L



singularity := proc (L) 
     local LL, v, i, alpha, sing, DD,X; sing := {}; DD:=_Envdiffopdomain[1]; X:=_Envdiffopdomain[2];
     LL := collect(L/lcoeff(L, DD), DD, factor); 
     v := factors(gcd(denom(LL), denom(DEtools[LCLM](DD-17, LL))))[2]; 
     for i to nops(v) do 
           alpha := RootOf(v[i][1], X); 
           sing := `union`(sing, {alpha}) 
     end do; 
     sing := `union`(sing, {infinity}) 
end proc:

# A:=singularity(L);

numberofint := proc (s) 
   local i, a; a := 0; 
   for i to nops(s) do 
       if type(s[i], integer) = false 
            then a := a 
                 else a := a+1 
       end if 
   end do; a 
end proc:

removable := proc (L) 
     local sing, i, rev, a, b, c,DD,X; 
     DD:=_Envdiffopdomain[1]; X:=_Envdiffopdomain[2]; 
     rev := {}; 
     sing := singularity(L); 
     for i in sing do 
          a := DEtools[gen_exp](L, T, X = i); 
          b := {seq(op(j[1 .. -2]), j = a)}; 
          c := [seq(op(j[1 .. -2]), j = a)]; 
          if nops(b) <> nops(c) then  # there are two exponents equal, so the solution has log, so nonremovable
                 rev := rev 
                 elif numberofint({seq(seq(m-n, m = b), n = b)}) <> nops({seq(seq(m-n, m = b), n = b)}) # there exists nonzero exponent diff, nonremovable
                         then rev := rev 
                                 elif DEtools[formal_sol](L, `has logarithm?`, X = i) # has log, then nonremovable
                                         then rev := rev 
                                                else rev := `union`(rev, {i}) 
          end if 
     end do; rev 
end proc:

#A:=A minus removable(L);
# from these singularities, find possible Mobius transformations u and v.

# pick 3 different element from given set A.
# A is the set of all singularites of the given operator.
# output is like {[a,b,c],[d,e,f],...}
threesing:=proc(A)    
   local i, j, k, B; B:={};
   for i in A do
     for j in A minus {i} do
        for k in A minus {i,j} do
          B:= B union {[i,j,k]}
        od;
     od;
   od; B;
end:



# Mob gives Mobius transformation which send {p,q,r} to {0,1,infinity}.
# input are 3 points

Mob:=proc(p,q,r)    
   #options trace;
   local s;
   if {p,q,r} intersect {infinity} = {} 
      then solve({a*p+b, c*r+d,a*q+b-c*q-d},{a,b,c,d});
           s:=subs(%,(a*x+b)/(c*x+d));           
           elif p=infinity 
                then solve({b-c*q-d, c*r+d}, {b,c,d});
                     s:=subs(a=0,%,(a*x+b)/(c*x+d));
                     elif q=infinity 
                          then solve({a*p+b, a*r+d}, {a,b,d});
                               s:=subs(c=a,%,(a*x+b)/(c*x+d));
                               else solve({a*p+b, a*q+b-d}, {a,b,d});
                                    s:=subs(c=0,%,(a*x+b)/(c*x+d));
   fi; normal(s);
end:



# SetofMob gives the set of all possible Mobius transformations whose p,q,r are in singularities.
# input is the set of singularities of the given operator L
# output is like {f1,f2,f3,...}
SetofMob:=proc(A)   
   local B,C,i;C:={};
   B:=threesing(A);
   for i in B do
      C:=C union {Mob(op(i))};
   od;
end:



# find roots of f,1-f and 1/f.
# input is a function of degree 1.
# output is like {a,b,c}
sing:=proc(f)   
   if degree(numer(f))>degree(denom(f))
      then {solve(numer(f),x),solve(f-1,x),infinity}
      elif degree(numer(f))<degree(denom(f))
           then {solve(denom(f),x),solve(f-1,x),infinity}
           elif lcoeff(denom(f))=lcoeff(numer(f))
                then {solve(numer(f),x),infinity,solve(denom(f),x)}
                else {solve(numer(f),x),solve(f-1,x),solve(denom(f),x)}
   fi;
end:


# find roots of u-v (u not equal to v since in F1 is reducible if u equals to v)
# input are 2 functions of degree 1.
# output is like {a}
rootof:=proc(u,v)
   local n, d; if simplify(u-v)=0 then return "error u v equal"; fi;
   d:=u-v; d:=normal(d); n:=degree(numer(d))-degree(denom(d));
   if solve(d,x)=NULL then {infinity}
        elif n<0 then {solve(d,x)} union {infinity} else {solve(d,x)}
   fi;
end:


rootof:=proc(u,v)
   #options trace;
   local n, d,s, S;  
   if indets(simplify(u-v))={} then S:={};  # u-v is a non-zero constant. then only x=infinity may be root of u-v=0. 
                                                    # x=infinity being root of u-v only happens if u=infinity when x=infinity
                                                    # so x=infinity is already in the singularity of L since we collect the root of u=infinity    
   else 
   d:=u-v; d:=normal(d); n:=degree(numer(d))-degree(denom(d));
   if solve(d,x)=NULL then S:={infinity}  # d = const/(x+b) or const/(x+b)/(x+c)
        elif n<0 then S:={solve(d,x)} union {infinity} # d= (ax+b)/(x^2+...)
        else S:={solve(d,x)}; 
             
    fi;  fi;
   S;
end:


# find all candidates for u and v s.t. roots of u, 1-u, 1/u, v,1-v,1/v and u-v are in the singularities and union of them and infinity equals singularity set.
# input is a set of functions from "SetofMob"
# output is like {[u,v],[u1,v1],...}
CandiUV:=proc(M,A)   
   #options trace;
   local C,B,c, SS, i, j, n; 
   C:={}; 
   for i in M do
       for j in M minus {i} do
           C:=C union {[i,j]};
       od;
   od;
   for c in C do
       if sing(c[1]) union sing(c[2]) union rootof(c[1],c[2])  <> A then C:=C minus {c}; fi;
   od; C;
 end: