with(DEtools):
read "D:/Homomorphisms.txt";

# Input: two C(x,y)[Dx,Dy]-modules M1 and M2
# Output: one dimensional module OM s.t. tensorProduct(M1,OM) homomorphic to M2

# Step 1. Find cyclic vector of M1 with respect to x, then find its minimal operator w.r.t x Say L1x. Same for M2: L2x. Likewise for y: L1y and L2y.
# Step 2. (Do not use command in Maple) For L1x, multiple lcm(denoms of L1x) to make it in C[x,y][Dx], take its coeff of highest order Dx^n as anx.
#          Do same to L2x--bnx; L1y--any;  L2y--bny.  Now factor anx*any*bnx*bny. 
# Step 3. For factors including x, find exponent difference between L1x and L2x at x=RootOf(factor,x). (also include infinity)
#         For factors including y, find exponent difference between L1y and L2y at y=RootOf(factor,x). (also include infinity)
# Step 4. For factors including both x and y, intersect exp-diff set from x and that from y.
# Step 5. Candidate for functions f should be a multiple of any factor with power in diff-exp set.
# Step 6. For any such candidate f, obtain one dimension module of f, say OM. If tensor product of M1 and OM is homomorphic to M2, then f is good.

projHom:=proc(M1,M2,opt:=0)
        local k,vars,h, j, L1, L2, i, d, an, bn, product, fac,f,p,gen1,gen2,s1,s2,dset,newS, s, checkS, A, Differ, newDiff, fSet, C, SS, ss, Candi, candi,F,M,den,newL1,newL2,S,Differ1,Differ2,OM;        
        vars:=M1[1];      # variables
        for j to nops(vars) do   # for each variable, do
                L1[j]:=CycVec(M1,vars[j])[2];      # minimal operator of the cyclic vector of M1
                for i from 0 to degree(L1[j],M1[2][j]) do
                    d[i]:=denom(coeff(L1[j],M1[2][j],i))   # denominators in the coefficient of the operator 
                od;
                den:=lcm(seq(d[i],i=0..degree(L1[j],M1[2][j])));  
                newL1[j]:=simplify(den*L1[j]);             # all coefficients of the new operator are rational functions 
                an[j]:=coeff(newL1[j],M1[2][j],degree(L1[j],M1[2][j]));  # leading coefficient of the new operator

                L2[j]:=CycVec(M2,vars[j])[2];              # same to the module M2
                for i from 0 to degree(L2[j],M1[2][j]) do
                           d[i]:=denom(coeff(L2[j],M1[2][j],i))
                od;
                den:=lcm(seq(d[i],i=0..degree(L2[j],M1[2][j])));
                newL2[j]:=simplify(den*L2[j]);
                bn[j]:=coeff(newL2[j],M1[2][j],degree(L2[j],M1[2][j]));
        od;
        p:=mul(an[j],j=1..nops(vars))*mul(bn[j],j=1..nops(vars));  # multiply these leading coefficients together
        fac:=factors(p)[2]; # factors of these leading coeffs
        for j to nops(vars) do
            S[j]:={};  
            k:={1,2} minus {j}; # here we assume only two variables in D
            _Envdiffopdomain := [M1[2][j],M1[1][j]];
               for f in fac do
                     if has(f[1],vars[j]) then  # get the exponents of the operator at these singularities (roots of leading coefficients)
                           gen1:=gen_exp(L1[j],T,vars[j]=RootOf(f[1],vars[j]));
                           gen2:=gen_exp(L2[j],T,vars[j]=RootOf(f[1],vars[j]));
                  if opt=0 then   # option = 0, general case, not adding any options 
                           s1:=ram(gen1); # exponents with ramifications
                           s2:=ram(gen2);
                           dset:=[f[1],findExp(s1,s2)]; # factor in f and its possible exponents
                  else # adding some options, now only consider all ramifications = 1 (most cases)  
                           s1:={seq(seq(m[i],i=1..nops(m)-1),m=gen1)};
                           s2:={seq(seq(m[i],i=1..nops(m)-1),m=gen2)};
                     if opt = 1 then  # restrict to surjective
                           dset:=[f[1],diffSurj(s1,s2)];
                     elif opt = 2 then # restrict to injective 
                           dset:=[f[1],diffInj(s1,s2)];
                     else "error: type a right option value from 0 to 2";
                     fi;
                  fi; 
                           S[j]:=S[j] union {dset}; fi;
                     
               od;     
          #     gen1:=gen_exp(L1[j],T,vars[j]=infinity);
          #     gen2:=gen_exp(L2[j],T,vars[j]=infinity);
          #     s1:=ram(gen1);
          #     s2:=ram(gen2);
          #     dset:=[1/vars[j],findExp(s1,s2)]; S[j]:=S[j] union {dset};
               newS[j]:=S[j];
               for s in S[j] do
                     if has(s[1],vars[k[1]])   # factor has both variables, then subtract that element from S[j]
                           then newS[j]:=newS[j] minus {s} fi;
               od;
               checkS[j]:=S[j] minus newS[j];  # checkS[j] set consists of all factors which have both variables, they need to be checked furthur
        od;
        A:=newS[1] union newS[2];     # factors containing only one variable
        for s in checkS[1] do         # check if the candidate exponents difference obtained from x part match difference obtained from y part
               Differ1:=s[2];         # candidates exponents difference obtained from x part    
               for ss in checkS[2] do   
                      if s[1]=ss[1] 
                           then Differ2:=ss[2]; # obtain candidate exponents difference from y part
                           break; fi;
               od;
               newDiff:={};      
               for i in Differ1 do
                     for j in Differ2 do
                           if type(i-j,integer)  # exponents difference match if they differ by integers
                               then newDiff:=newDiff union {i}; fi;
                     od;
               od;
               newDiff:=modz(newDiff);
               if newDiff<>{} 
                    then A:=A union {[s[1],newDiff]}; # now A contains all factors and their candidate exponents 
               fi;
        od;
        
        for i to nops(A) do
               fSet[i]:={};
               for j to nops(A[i][2]) do
                     fSet[i]:=fSet[i] union {(A[i][1])^(A[i][2][j])};
               od;
        od;
        SS:={seq(fSet[i],i=1..nops(A))};  # sets of all candidate factors f^e for r  
        C:=combi(SS);  # a set of sets, each set contains each factor^e once
        Candi:=simplify({seq(combii(mul(c[i],i=1..nops(c)),{op(vars)}),c=C)}); # all candidate r 
        F:={}; lprint(nops(Candi));
        for candi in Candi do
               OM:=OneDiModule(candi,[x,y]);  # 1-dimensional module for r
               M:=TProModule(M1,OM);          # tensor product 
               h:=hom(M,M2);
               if h<>{}               
                     then return [candi,h]; 
                     break; 
               fi;
         od;

end:


# Input: a set of sets say A={S1,S2,...,Sn}, Si is a set.
# Output: all possible combinations of {{s1,s2,...,sn},...} si in Si.
combi:=proc(A)
  local B;
     if nops(A)=0 then return {{}}
        elif nops(A)=1 then return {seq({s},s=A[1])};
        else B:=A minus {A[1]};
             return {seq(seq(i union j, i=procname(B)), j=procname({A[1]}))};
     fi;
end:

     
               
       
combii:=proc(a,vars)
        local b,v,i,j,B,ra;
if type(a,`*`) then
        ra:=proc(j) `if`(has(j,{'RootOf','Int','int'}),j,radnormal(j)) end:
        B:=evalb(nargs=2 and nops(v)=1);
        b:=[seq(`if`(type(i,`^`),[op(i)],[i,1]),i=select(has,[op(a)],vars))];
        b:=[seq(`if`(type(i[1],`^`),  [op(1,i[1]), op(2,i[1])*i[2]], i), i=b)];
        for i from 1 to nops(b)-1 do for j from i+1 to nops(b) do
                if not has(normal(b[i][1]/b[j][1]),vars) then
                        b := subsop(i = [1,1], j = [b[j][1], b[i][2] + b[j][2]], b)
                fi
        od od;
        v:={seq([i[2],-i[2]],i=b)};
        mul(`DEtools/kovacicsols/simpl_prod`(ra(mul(`if`(j[2]=i[1],j[1],`if`(j[2]=i[2]
         ,1/j[1],1)),j=b)),vars,ra(i[1]),B),i=v)
else
        a
fi
end:

      

# Input:   a list or set of numbers
# Output:  the set which is equal to the input mod Z
# modz([1,2,1/2,3/2])={1,1/2}

modz := proc (s)
 local i, s1, j, news; 
 news:=[op(s)]; # list of all elements in s
 s1:=s;
 if nops(s) = 0 or nops(s) = 1 
    then return {op(s)}; 
    else for i to nops(s)-1 do  # compare all pairs 
              for j from i+1 to nops(s) do
                    if type(news[i]-news[j], integer) 
                       then if i=1 then
                         s1 := {op(news[2..-1])}; 
                         else s1 := {op(news[1..i-1]),op(news[i+1..-1])}
                       fi; s1:=procname(s1);
                    fi; 
               od;
    od; s1; 
fi;
end proc:



# Input:  two sets of exponents lists(output of "ram"-- with ramifications). 
# Output: the exponent difference between input 
# findExp({[0, 1], [1-b1, 1], [1-b2, 1]},{[0, 1], [1, 1], [2, 1]})={0, b1-1, b2-1}
findExp:=proc(s1,s2)
    local l1,l2,i,j,e,f,d,P1,P2,candiExp; 
    candiExp:={};
    l1:=[seq(i,i=s1)]; l2:=[seq(i,i=s2)]; # convert sets into lists
    for i to nops(l1) do
        e[i]:=l1[i];
    od;
    for i to nops(l2) do
        f[i]:=l2[i];
    od;
    for i to nops(l1) do
        for j to nops(l2) do
            if f[j][2]=e[i][2] and f[j][2]=1 then  # there exists an exponents in l2 which has same ramification as in l1 and equal to 1
                  d:=f[j][1]-e[i][1];  candiExp:=candiExp union {d}; 
                  elif f[j][2]=e[i][2] then        # there exists an exponents in l2 which has same ramification as in l1 but greater than 1
                          P1 := collect(evala(Norm(X-e[i][1], {}, {})), X, evala);  # minpoly of the exponent in the field extension from C(tp) to C(T) 
                          P2 := collect(evala(Norm(X-f[i][1], {}, {})), X, evala);  # minpoly
                          if f[j][2]=2 then d := (1/2)*coeff(P1-P2, X, 1);   candiExp:=candiExp union {d};   # ramification = 2
                             elif f[j][2]=3 then d:=(1/3)*coeff(P1-P2, X, 2);  candiExp:=candiExp union {d}; # ramification = 3
                          fi;
            fi;
        od;
    od;
    modz(candiExp);
end:





# Input: generalized exponent list at some singularity; (output of gen_exp(L,T,x=xi)).
# Output: a set of lists, each list consists of two entries, first is the generalized exponent(in terms of T), second is its ramification.
# ram([[0, T = x], [1-b1, T = x], [1-b2, T = x]])={[0, 1], [1-b1, 1], [1-b2, 1]}

ram:=proc(l)
   #options trace;
   local i, e, d, r, ll, Tvalue,s; 
   ll:={};
   for i to nops(l) do
       e[i]:=l[i];
       d[i]:=degree(op(e[i][-1])[1],T);  # degree of T^i part
       r[i]:=1/d[i];                     # ramification
       if d[i]=1 then                    # ramification is 1
              ll:=ll union {seq([e[i][j],1],j=1..nops(e[i])-1)};
              else  Tvalue:=RootOf(subs(T=_Z,e[i][2]));  # ramification not 1
                    ll:=ll union {[simplify(subs(T=Tvalue,e[i][1])),r[i]]};
       fi;
   od;
   return ll;
end:



# Input:  two gen-exp sets A and B.e.g {0,1,c-b1} and {0,b1-c}. Now we only consider regular singular case, i.e., ramification=1
# Output:  candidate exponents differences s.t. there exists a surj from A to B.
# diffSurj({0,1,c-b1},{0,b1-c})={-1+b1-c}

with(DEtools):
with(combinat):

diffSurj:=proc(A,B)
      local candi, newA, newB,n,m,lset,l,ll,diff;
      candi:={};
      if nops(B)>nops(A) # there is no surjective
         then return {} 
      fi;
      newA:=modz(A); newB:=modz(B);   # sets mod Z
      n:=nops(newA); m:=nops(newB);
      if m>n then return {}; fi;
      lset:=permute(A,nops(B));     # in combinat package, gives all combinations of nops(B) elements in A
      l:=[seq(i,i=B)];       
      for ll in lset do
                diff:=modz(l-ll);   
                if nops(diff)=1 then candi:=candi union {op(diff)}; fi;
      od;
      modz(candi);
end:


# Input: same as input in diffSurj
# Output: candidates of gen_exp diff s.t. there exists a injective from A to B.
# This is same as diffSurj(B,A) with elements inversed.

diffInj:=proc(A,B)
      if diffSurj(B,A)={} then return {}; else return {seq(-i,i=diffSurj(B,A))}; fi;
end: