## This deals with the general case: there may exist some generalized exponents with ramification > 1. (Then the operator is irregular singular)
## To test if two order 3 operators are projectively equivalent. 


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


# Input: a set or a list s
# Output: a list = s mod Z

modz := proc (s) 
        local i, s1, j; 
        s1 := s;
if nops(s) = 0 or nops(s) = 1 then
        return s
end if; 
for i to nops(s)-1 do 
     for j from i+1 to nops(s) do 
        if type(s[i]-s[j], integer) then
                if i = 1 then
                         s1 := [op(s[2 .. -1])] 
                         else s1 := [op(s[1 .. i-1]), op(s[i+1 .. -1])] 
                end if; s1 := procname(s1) 
        end if 
     end do 
end do; s1 
end proc:


# Input: a list or a set s
# Output: an element in s which is equal to another element in s mod Z

recur := proc (s) 
         local i, j; 
for i to nops(s)-1 do 
         for j from i to nops(s) do
              if type(s[i]-s[j], integer) then
                  return s[i]
              end if
         end do
end do 
end proc:



# Input: a list or a set s
# Output: the number of integers in s

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:


# Input: a differential operator L
# Output: the singularity set of L

singularity := proc (L) 
         local LL, v, i, alpha, sing; sing := {}; 
         LL := collect(L/lcoeff(L, Dx), Dx, factor); 
         v := factors(gcd(denom(LL), denom(LCLM(Dx-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:


# Input: a differential operator L
# Output: the removable singularity set of L

removable := proc (L) 
          local sing, i, rev, a, b, c; 
          rev := {}; 
          sing := singularity(L); 
for i in sing do 
          a := 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  # two exponents equal, then 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)}) then # non integer exponent difference,  
                                       # then nonremobable 
                rev := rev 
                elif formal_sol(L, `has logarithm?`, x = i) then # has log, 
                                                      # then nonremovable
                     rev := rev 
                     else rev := `union`(rev, {i}) 
   end if
end do; rev 
end proc:



# Input:  two lists or two sets
# Output: check if these two lists (sets) are equal mod Z.

equalmodz := proc (L1, L2) 
         local i, j, a, result, b, B, A; 
         A := []; 
         result := -1; 
if nops(L1) <> nops(L2) then
              result := 0; 
              return "not equal mod z" 
end if; 
for i to nops(L1) do 
         a[i] := []; B[i] := {} 
end do;
for i to nops(L1) do 
         for j to nops(L2) do 
               b[i][j] := L1[i]-L2[j]; 
               if type(b[i][j], integer) then
                     a[i] := [op(a[i]), L2[j]] 
               end if 
         end do 
end do; 
for i to nops(L1) do 
         if a[i] = [] then
               result := 0; 
               return "not equal mod Z" 
                     else A := [op(A), op(a[i])] 
         end if 
end do; 
if {op(A)} <> {op(L2)} then
          result := 0; 
          return "not equal mod Z" 
end if; 
for i to nops(L1) do 
          for j to nops(L1) do 
                if a[i] = a[j] then
                        B[i] := `union`(B[i], {j}) 
                end if 
          end do;
          if nops(a[i]) <> nops(B[i]) then
                      result := 0; 
                      return "not equal mod z" 
          end if 
end do; 
if result <> 0 then
           result := 1 
end if; result 
end proc:


# Input: two lists (sets) L1 and L2
# Output: the difference set {d} such that L1+d = L2 mod Z

findc := proc (L1, L2) 
          local i, j, c, T, result, L3, L4; 
          result := {}; 
          L3 := [seq(seq(i-j, i = L1), j = L1)]; 
          L4 := [seq(seq(i-j, i = L2), j = L2)];
if equalmodz(L3, L4) <> 1 then
          "error"
          else for i to nops(L2) do 
                   c[i] := L2[i]-L1[1]; 
                   T[i] := [L1[1]+c[i], L1[2]+c[i], L1[3]+c[i]]; 
                   if equalmodz(L2, T[i]) = 1 then
                           result := `union`(result, {c[i]}) 
                   end if 
          end do 
end if; 
          {op(modz([op(result)]))}
end proc:


# Input: the exponent of Dx-r at x=infinity
# Output: the polynomial part of r 

polypart := proc (e) 
          local r;
          r := collect(subs(T = 1/x, T*(coeff(e, T, 0)-e)), x) 
end proc:


# Input:  a set of {[singularity,[diffenrence of two operators at this singularity]]}
# Output: r in exp-product transformation which collect all differences at all singularities.

intsum := proc (S, basefield) 
          local s, i, j, u, p, candi; 
          candi := {}; 
if nops(S) = 0 then
          {}
          elif nops(S) = 1
               then for i in S[1][2] do 
                    candi := `union`(candi, {[[S[1][1], i], evala(Trace
                              (coeff(i, T, 0), `union`(basefield, {S[1]
                              [1]}), basefield))]}) 
               end do 
               else s := {S[1]}; u := `minus`(S, s); p := procname(u, 
                                                          basefield); 
                    for i in s[1][2] do 
                        for j in p do 
                           candi := `union`(candi, {[[s[1][1], i], seq(j
                                    [k], k = 1 .. nops(j)-1), evala(Trace
                                    (coeff(i, T, 0), `union`(basefield, {s
                                    [1][1]}), basefield))+j[-1]]}) 
                         end do 
                    end do 
end if; candi 
end proc:





findr := proc (L1, L2, basefield) 
 #options trace;
  local S, candi, r, singu1, singu2, singu, i, e1, e2, expo1, expo2, c, j, root1, root2, s1, s2, P1, P2, d1, d2, A, rpart, d;
  S := {}; candi := {}; r := 0; 
  singu1 := `minus`(singularity(L1), removable(L1)); 
  singu2 := `minus`(singularity(L2), removable(L2)); 
  singu := `union`(singularity(L1), singularity(L2)); 
  if singu1 <> singu2 then
      return "they are not projective equivalent" 
  end if; 
  for i to nops(singu) do 
      e1 := gen_exp(L1, T, x = singu[i]); 
      e2 := gen_exp(L2, T, x = singu[i]); 
      expo1 := [seq(op(j[1 .. -2]), j = e1)]; 
      expo2 := [seq(op(j[1 .. -2]), j = e2)];
      if nops(expo1) = nops(expo2) and nops(expo2) = 3 then
          if findc(expo1, expo2) = {} then
                  "they are not projective equivalent" 
                   else c := [singu[i], findc(expo1, expo2)]; 
                   S := `union`(S, {c}) 
          end if 
          elif nops(expo1) = nops(expo2) and nops(expo2) = 2 then
                 if has(expo1, _Z) and has(expo2, _Z) then
                       for j in e1 do 
                            if has(j[1], _Z) then
                                 root1 := j[1] 
                                      else s1 := j[1] 
                            end if
                       end do; 
                       for j in e2 do 
                             if has(j[1], _Z) then
                                  root2 := j[1] 
                                  else s2 := j[1] 
                             end if
                       end do; 
                       P1 := collect(evala(Norm(X-root1, {}, {})), X, evala); 
                       P2 := collect(evala(Norm(X-root2, {}, {})), X, evala); 
                       d1 := s2-s1; 
                       d2 := (1/2)*coeff(P1-P2, X, 1); 
                       if evala(subs(X = X+d2, P2)-P1) = 0 and type(evala(2*(d1-d2)), integer) then
                              c := [singu[i], {d1}];
                              S := `union`(S, {c}) 
                              else "they are not projective equivalent" 
                       end if; 
                                else for j in e1 do
                                         if degree(op(j[2])[1], T) = 2 then
                                             root1 := subs(T = RootOf(subs(T = _Z, j[2])), j[1]) 
                                             else s1 := subs(T = RootOf(subs(T = _Z, j[2])), j[1]) 
                                         end if 
                                      end do; 
                                      for j in e2 do
                                          if degree(op(j[2])[1], T) = 2 then 
                                              root2 := subs(T = RootOf(subs(T = _Z, j[2])), j[1]) 
                                              else s2 := subs(T = RootOf(subs(T = _Z, j[2])), j[1]) 
                                          end if 
                                      end do; 
                                      P1 := collect(evala(Norm(X-root1, {}, {})), X, evala); 
                                      P2 := collect(evala(Norm(X-root2, {}, {})), X, evala); 
                                      d1 := s2-s1; 
                                      d2 := (1/2)*coeff(P1-P2, X, 1); 
                                      if evala(subs(X = X+d2, P2)-P1) = 0 and type(evala(2*(d1-d2)), integer) then
                                           if singu[i] = infinity then
                                                  d1 := subs(x = 1/T, d1) 
                                                  else d1 := subs(x = T+singu[i], d1) 
                                           end if;
                                           c := [singu[i], {d1}]; 
                                           S := `union`(S, {c}) 
                                           else "they are not projective equivalent" 
                                      end if
            end if 
            elif nops(expo1) = nops(expo2) and nops(expo2) = 1 then
                    if has(expo1, _Z) and has(expo2, _Z) then
                          root1 := op(e1)[1]; 
                          root2 := op(e2)[1]; 
                          P1 := collect(evala(Norm(X-root1, {}, {})), X, evala); 
                          P2 := collect(evala(Norm(X-root2, {}, {})), X, evala); 
                          d := (1/3)*coeff(P1-P2, X, 2); 
                          if evala(subs(X = X+d, P2)-P1) = 0 then
                                 c := [singu[i], {d,d+1/3,d-1/3}]; 
                                 S := `union`(S, {c}) 
                                 else "they are not projective equivalent" 
                          end if 
                                    else root1 := subs(T = RootOf(subs(T = _Z, op(e1[1])[2])), op(e1[1])[1]);
                                         root2 := subs(T = RootOf(subs(T = _Z, op(e2[1])[2])), op(e2[1])[1]); 
                                         P1 := collect(evala(Norm(X-root1, {}, {})), X, evala);
                                         P2 := collect(evala(Norm(X-root2, {}, {})), X, evala); d := (1/3)*coeff(P1-P2, X, 2); 
                                         if evala(subs(X = X+d, P2)-P1) = 0 then
                                               if singu[i] = infinity then
                                                    d := subs(x = 1/T, d) 
                                                    else d := subs(x = T+singu[i], d) 
                                               end if; 
                                               c := [singu[i], {d,d+1/3,d-1/3}];
                                               S := `union`(S, {c}) 
                                                    else "they are not projective equivalent" 
                                         end if 
            end if
      end if
  end do;
  A := intsum(S, basefield); 
  for i in A do 
      if not type(i[-1], integer) then
          A := `minus`(A, {i}) 
      end if 
  end do; 
  for i in A do
       for j to nops(i)-1 do 
            if i[j][1] = infinity then
                 rpart := polypart(i[j][2]) 
                 else rpart := evala(Trace(subs(T = x-i[j][1], i[j][2]/T), `union`(basefield, {i[j][1]}), basefield)) 
            end if; r := r+rpart
       end do; candi := `union`({r}, candi) 
  end do; candi 
end proc:





projective := proc (LB, Linp) 
   local R, M, G, L, r, g, trans, basefield; 
   basefield := indets([args], {RootOf, nonreal, radical}); 
   trans := {}; 
   R := findr(LB, Linp, basefield); 
   for r in R do 
       M := symmetric_product(LB, Dx-r); 
       G := Homomorphisms(M, Linp); 
       for g in G do 
          L := rightdivision(LCLM(M, g), g)[1]; 
          if has(evala(L/Linp), Dx) then
               trans := trans 
               else trans := `union`(trans, {[r, g]})
          end if 
       end do 
   end do; trans 
end proc:



## Example, more examples in "eg2.mw" file
## L2 := (x-1)^2*Dx^3+x*Dx
## L := symmetric_product(L2, Dx+x^5)
## LL:=rightdivision(LCLM(L, Dx^2-Dx), Dx^2-Dx)[1]:
## projective(L2,LL);
## projective(L2,L);