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:


# Input: two operators L1 and L2 and a basefiele
# Output: the parameter r in exp-product transformation from L1 to L2

findr := proc (L1, L2, basefield) 
         local S, singu, singu1, singu2, i, e1, e2, expo1, expo2, candi, j,  
               r, rpart, c, A, root1, root2, s1, s2, d, d1, d2, P1, P2; 
         S := {}; candi := {}; r := 0; 
         singu := `minus`(singularity(L1), removable(L1)); 
         singu2 := singularity(L2); 
         singu1 := `minus`(singu2, removable(L2));
if singu <> singu1 then
         return "they are not projective equivalent" 
end if; 
for i to nops(singu2) do 
          e1 := gen_exp(L1, T, x = singu2[i]);
          e2 := gen_exp(L2, T, x = singu2[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 := [singu2[i], findc(expo1, expo2)]; 
                             S := `union`(S, {c}) 
                 end if 
          elif nops(expo1) = nops(expo2) and nops(expo2) = 2 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 := [singu2[i], {d1}];
                               S := `union`(S, {c}) 
                       else "they are not projective equivalent" 
                end if; 
         elif nops(expo1) = nops(expo2) and nops(expo2) = 1 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 := [singu2[i], {d,d-1/3,d+1/3}]; S := `union`(S, 
                                                              {c}) 
                       else "they are not projective equivalent" 
                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:



# Input: two operators LB and Linp
# Output: the exp-product transformation and gauge transformation which 
          #send LB to Linp

projective := proc (LB, Linp,basefield) 
        local R, M, G, L, r, g, trans; 
        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:




# Input:  an operator Linp and a basefiled
# Output: the operator after applying a change of variables transformation 
          #on Linp

findM := proc (Linp, basefield) 
         local eq, sing, sing1, sing2, singu, i, e, e1, expo, candi, 
               candi1, candi0, c, m, R1, R2, R3, f2, f3, expo1, j, u1, k, 
               l, expo0, B, soln, b, infi, expoinfi, f, f1,g1, g2, g3, d, 
               a1, a2, a3, a4, a5, a6, candidate, singu10infi; 
         expo := {}; candi1 := {}; candidate := {}; u1 := 0; 
         singu10infi := {}; 
         sing := `minus`(singularity(Linp), removable(Linp));
if nops(sing) = 3 then
         singu := sing 
elif nops(sing) = 2 and has(sing, _Z^2) then
         for i in sing do
                if has(i, _Z^2) then
                           R1 := i 
                           else R3 := i 
                end if 
         end do; 
         R2 := evala(Trace(R1, `union`(basefield, {R1}), 
                      basefield)-R1); 
         singu := {R1, R2, R3}
elif nops(sing) = 1 and has(sing, _Z^3) then
          R1 := op(sing); 
          f1 := subs(_Z = x, op(R1)); 
          f2 := evala(Quo(f1, x-R1, x)); 
          f3 := evala(Factors(f2, `union`(basefield, {R1})))[2]; 
          if nops(f3) = 2 then
                    singu := {RootOf(f3[1][1], x), RootOf(f3[2][1], x), R1} 
                    else R2 := RootOf(f3[1][1], x); 
                         R3 := evala(Trace(R2, `union`(basefield, {R1, 
                               R2}), `union`(basefield, {R1}))-R2); 
                         singu := {R1, R2, R3} 
          end if 
else return "not solvable" 
end if; 
 
for i in singu do 
          e := gen_exp(Linp, T, x = i);
          if nops(e) = 3 or not has(e, _Z) then
                 e1 := [seq(op(j[1 .. -2]), j = e)] 
          elif nops(e) = 2 then
                 for j in e do
                        if has(j[1], _Z) then
                              a4 := j[1] 
                              else a5 := j[1] 
                        end if
                 end do; 
                 a6 := evala(Trace(a4, `union`(basefield, {a4}), 
                       basefield)-a4); 
                 e1 := [a4, a6, a5] 
          elif nops(e) = 1 then
                  a4 := e[1][1];
                  f1 := subs(_Z = x, op(a4)); 
                  f2 := evala(Quo(f1, x-a4, x)); 
                  f3 := evala(Factors(f2, `union`(basefield, {a4})))[2];
                  if nops(f3) = 2 then
                    singu := {RootOf(f3[1][1], x), RootOf(f3[2][1], x), a4}
                    else a5 := RootOf(f3[1][1], x); 
                         a6 := evala(Trace(a5, `union`(basefield, {a4, 
                               a5}), `union`(basefield, {a4}))-a5); 
                         e1 := [a4, a5, a6] 
                  end if 
         end if; expo := `union`(expo, {[i, e1]}) 
  end do; 
  for i in expo do
         if nops(modz(i[2])) < nops(i[2]) then
             candi1 := `union`(candi1, {i}) 
         end if 
  end do; 
  for i in candi1 do
         u1:=0;
         candi0 := `minus`(expo, {i});
         c := recur(i[2]);
         expo1 := i[2]+[-(`$`(c, 3))];
         for j in expo1 do
               if type(j, integer) then
                      u1 := u1 
                      else u1 := j+u1 
               end if 
         end do; 
         for j in candi0 do 
             for k in j[2] do 
                expo0 := j[2]+[-(`$`(k, 3))]; 
                B := deleteone0(expo0); 
                soln := solve({1-b1 = B[1], 1-b2 = B[2]}, {b1, b2});
                b := [op(soln[1])[2], op(soln[2])[2]]; 
                infi := `minus`(candi0, {j}); 
                expoinfi := op(infi)[2];
                f := MobiustoPQW(j[1], i[1], op(infi)[1]); 
                g1 := expoinfi[1]; 
                g2 := expoinfi[2];
                g3 := expoinfi[3]; 
                d := (1/3)*b[1]+(1/3)*b[2]-(1/3)*u1-(1/3)*g1-(1/3)*g2-(1/3)
                     *g3; 
                a1 := [g1+d, g2+d, g3+d]; 
                a2 := [g1+d+1/3, g2+d+1/3, g3+d+1/3]; 
                a3 := [g1+d+2/3, g2+d+2/3, g3+d+2/3]; 
                candidate := `union`(candidate, {[a1, b, f], [a2, b, f], 
                                      [a3, b, f]}) 
          end do 
    end do 
end do; 
candi := candidate; 
    for i in candidate do 
          if `or`(has(i[1], 0), has(i[2], 0)) then
                 candi := `intersect`(candi, `minus`(candidate, {i}))
          end if 
    end do; candi
end proc:



# Input: three different points
# Output: the mobius transforamtion which send  (P,Q,W) to (0,1,infinity)

MobiustoPQW := proc (P, Q, W)
  local a, b, c, d, eq1, eq2, eq3; 
  eq1 := `if`(P = infinity, a, a*P+b); 
  eq2 := `if`(Q = infinity, a-c, (a-c)*Q+b-d); 
  eq3 := `if`(W = infinity, c, c*W+d); 
  solve({eq1, eq2, eq3}, {a, b, c, d}); 
  subs(%, (a*x+b)/(c*x+d)); normal(%)
end proc:


# Input:  a list s
# Output: a list which deletes one 0 in s

deleteone0 := proc (s)
  if s[1] = 0 then
       return [s[2], s[3]] 
       elif s[2] = 0  then
            return [s[1], s[3]] 
            else return [s[1], s[2]]
  end if 
end proc:



# Input: an operator Linp
# Output: the parameters in 3F2 function and the parameters in 
          #transformations which send L to Linp, L is the minop of 3F2.

transformations := proc (Linp)
        local candidates, H, M, L, L1, i, j, rg, r, g, trans, basefield,a;
        basefield := indets([args], {RootOf, nonreal, radical}); 
        trans := {}; 
        candidates := findM(Linp, basefield);
for i in candidates do 
        if `or`(`and`(type(i[2][1], integer), i[2][1] < 0), `and`(type(i[2]
            [2], integer), i[2][2] < 0)) then
          if `and`(`and`(`not`(`and`(type(i[1][1], integer), i[1][1] < 0)), 
            `not`(`and`(type(i[1][2], integer), i[1][2] < 0))), `not`(`and`
                 (type(i[1][3], integer), i[1][3] < 0))) then
                   a := subs(i[2] = [i[2][1]-floor(i[2][1])+1, i[2][2]-
                        floor(i[2][2])+1], i); 
                   candidates := subs(i = a, candidates) 
          end if 
         end if 
end do;
for i in candidates do 
      H := hypergeom(i[1], i[2], i[3]); 
      M := de2diffop(hyperode(H, y(x)), y(x)); 
      rg := projective(M, Linp, basefield);
      for j in rg do 
          r := j[1];
          g := j[2]; 
          L := symmetric_product(M, Dx-r);
          L1 := rightdivision(LCLM(L, g), g)[1]; 
          if has(evala(L1/Linp), Dx) then
              trans := trans 
              else trans := `union`(trans, {[i[1], i[2], i[3], r, g]})
          end if
       end do
   end do; trans 
end proc:



## Example to test if two order 3 operators are projective equivalent

##  H := hypergeom([1, 2, 3], [1/2, 1/3], (x+1)/(x-5))
##  L1 := de2diffop(hyperode(H, y(x)), y(x))
##  M := symmetric_product(L1, Dx-1/2+x);
##  g := Dx
##  L2 := rightdivision(LCLM(M, g), g)[1]
##  projective(M, L2, {})
##  projective(L1, L2, {});

## And some other examples could be found in "eg.mw" file


## Example to find 3F2-type solutions of an irreducible order 3 operator with degree 1 ## pullback functions (rational function f in change of varaibles transformation).

###  L := Dx^3+(1/16*(85*x-61))*Dx^2/(x*(x-1))+(1/16*(78*x^2-131*x+29))*Dx/(x^2*(x-1)^2)+(7/64*(3*x-7))/(x^2*(x-1)^2);
###  Transformations(L);