# case2 deals with the case of A={rootof(u,1-u,1/u)}<>{rootof(v,1-v,1/v)}=B

case2:=proc(LL,u,v)
       local GG,l1,l2,g,e,C,B1,B2,A,L1,Lp,a1,a2,r,i,l,e1,L2,s,e2,S,SS,L,sb1,sb2,SB1,SB2;
       S:={}; SS:={};  
       L:=collect(LL/coeff(LL,Dx^(degree(LL,Dx))),Dx,factor); 
       GG:=Gen(L); l1:=sing2(u); l2:=sing2(v); lprint(l1); lprint(l2); L1:=EQN(F1(a,b1,b2,c,u,v),x);
 # case1: rootof(1/u) not in B and rootof(1/v) not in A
                   
       if not has({seq([i],i=l1)},{[l2[3]]}) and not has({seq([i],i=l2)},{[l1[3]]})  then lprint("case1");    
            # has([infinity, 1, RootOf(_Z^2+_Z-1)],-1) returns true, which should be avoided.
             for g in GG do
                 if g[-1]=l1[3] then e1:=MatchExp1(g[1],g[2],g[3],a,b1,1+b1); GG:=GG minus {g}; break;fi; 
             od; if e1={} then  return "no solution" fi; 
             for g in GG do
                 if g[-1]=l2[3] then e2:=MatchExp1(g[1],g[2],g[3],a,b2,1+b2); break;fi; 
             od; if e2={} then return "no solution" fi; 
             e1:={seq(`if`(indets(i)={},NULL,i),i=op(op(e1)))};
             e2:={seq(`if`(indets(i)={},NULL,i),i=op(op(e2)))};
             B1:=solve(op(e1),b1);
             B2:=solve(op(e2),b2);
             L2:=subs(b2=B2,b1=B1,L1);
             L2:=collect(normal(L2),{Dx,Dx^2,Dx^3},factor);
             a1:=coeff(L,Dx^2); a2:=coeff(L2,Dx^2);
             r:=1/3*(a1-a2);
             Lp:=L2-DEtools[symmetric_product](L,Dx-r);
             Lp:=collect(Lp,Dx,normal):
             s:=[solve({coeffs(numer(Lp),{Dx,x})},{a,c})]; 
             if s<>[] then for i in s do 
                 A:=op(op(i)[1])[2]; C:=op(op(i)[2])[2]; B1:=subs(a=A,c=C,B1); B2:=subs(c=C,a=A,B2); S:=S union {[A,B1,B2,C,u,v]};
             od;fi;  
  # case2: rootof(u) not in B and rootof(v) not in A 
             
       elif not has({seq([i],i=l1)},{[l2[1]]}) and not has({seq([i],i=l2)},{[l1[1]]}) then lprint("case2");
            for g in GG do
                 if g[-1]=l1[1] then e1:=MatchExp1(g[1],g[2],g[3],0,1,1-c+b2); GG:=GG minus {g}; break;fi; 
            od; if e1={} then  return "no solution" fi; 
            for g in GG do
                 if g[-1]=l2[1] then e2:=MatchExp1(g[1],g[2],g[3],0,1,b1-c+1); break;fi;            
            od; if e2={} then return "no solution" fi; 
            e1:={seq(`if`(indets(i)={},NULL,i),i=op(op(e1)))};   
            e2:={seq(`if`(indets(i)={},NULL,i),i=op(op(e2)))};    
            B2:=solve(op(e1),b2);
            B1:=solve(op(e2),b1);
            L2:=subs(b2=B2,b1=B1,L1);
            L2:=collect(normal(L2),{Dx,Dx^2,Dx^3},factor);
            a1:=coeff(L,Dx^2); a2:=coeff(L2,Dx^2);
            r:=1/3*(a1-a2);
            Lp:=L2-DEtools[symmetric_product](L,Dx-r);
            Lp:=collect(Lp,Dx,normal):
            s:=[solve({coeffs(numer(Lp),{Dx,x})},{a,c})]; 
            if s<>[] then for i in s do 
                 A:=op(op(i)[1])[2]; C:=op(op(i)[2])[2]; B1:=subs(a=A,c=C,B1); B2:=subs(c=C,a=A,B2); S:=S union {[A,B1,B2,C,u,v]};
            od; fi;  
  # case3: rootof(u) not in B and rootof(1-v) not in A
         
       elif not has({seq([i],i=l1)},{[l2[2]]}) and not has({seq([i],i=l2)},{[l1[1]]}) then lprint("case3");
                 for g in GG do
                     if g[-1]=l1[1] then e1:=MatchExp1(g[1],g[2],g[3],0,1,1-c+b2); GG:=GG minus {g}; break;fi; 
                 od; if e1={} then  return "no solution" fi; 
                 for g in GG do
                     if g[-1]=l2[2] then e2:=MatchExp1(g[1],g[2],g[3],0,1,c-a-b2); break;fi;            
                 od;   if e2={} then return "no solution" fi; 
                 e1:={seq(`if`(indets(i)={},NULL,i),i=op(op(e1)))}; 
                 e2:={seq(`if`(indets(i)={},NULL,i),i=op(op(e2)))};  
                 B2:=solve(op(e1),b2);
                 A:=solve(subs(b2=B2,op(e2)),a);
                 L2:=subs(b2=B2,a=A,L1);
                 L2:=collect(normal(L2),{Dx,Dx^2,Dx^3},factor);
                 a1:=coeff(L,Dx^2); a2:=coeff(L2,Dx^2);
                 r:=1/3*(a1-a2);
                 Lp:=L2-DEtools[symmetric_product](L,Dx-r);
                 Lp:=collect(Lp,Dx,normal):
                 s:=[solve({coeffs(numer(Lp),{Dx,x})},{b1,c})]; 
                 if s<>[] then for i in s do 
                    B1:=op(op(i)[1])[2]; C:=op(op(i)[2])[2]; B2:=subs(c=C,B2); S:=S union {[A,B1,B2,C,u,v]};
                 od; fi;            
    # case4: rootof(u) not in B and rootof(1/v) not in A
      
         elif not has({seq([i],i=l1)},{[l2[3]]}) and not has({seq([i],i=l2)},{[l1[1]]})  then lprint("case4");
                      for g in GG do
                          if g[-1]=l1[1] then e1:=MatchExp1(g[1],g[2],g[3],0,1,1-c+b2); GG:=GG minus {g}; break;fi; 
                      od; if e1={} then  return "no solution" fi; 
                      for g in GG do
                          if g[-1]=l2[3] then e2:=MatchExp1(g[1],g[2],g[3],a,b2,1+b2); break;fi;            
                      od;   if e2={} then return "no solution" fi; 
                      e1:={seq(`if`(indets(i)={},NULL,i),i=op(op(e1)))}; 
                      e2:={seq(`if`(indets(i)={},NULL,i),i=op(op(e2)))}; 
                      C:=solve(op(e1),c);
                      A:=solve(op(e2),a);
                      L2:=subs(a=A,c=C,L1);
                      L2:=collect(normal(L2),{Dx,Dx^2,Dx^3},factor);
                      a1:=coeff(L,Dx^2); a2:=coeff(L2,Dx^2);
                      r:=1/3*(a1-a2);
                      Lp:=L2-DEtools[symmetric_product](L,Dx-r);
                      Lp:=collect(Lp,Dx,normal):
                      s:=[solve({coeffs(numer(Lp),{Dx,x})},{b1,b2})]; 
                      if s<>[] then for i in s do 
                         B1:=op(op(i)[1])[2]; B2:=op(op(i)[2])[2]; A:=subs(b2=B2,A); C:=subs(b2=B2,C); S:=S union {[A,B1,B2,C,u,v]};
                      od; fi;  
    # case5: rootof(1-u) not in B and rootof(v) not in A    
              
         elif not has({seq([i],i=l1)},{[l2[1]]}) and not has({seq([i],i=l2)},{[l1[2]]}) then lprint("case5");
                          for g in GG do
                             if g[-1]=l2[1] then e1:=MatchExp1(g[1],g[2],g[3],0,1,b1-c+1); GG:=GG minus {g}; break;fi; 
                          od; if e1={} then return "no solution" fi; 
                          for g in GG do
                             if g[-1]=l1[2] then e2:=MatchExp1(g[1],g[2],g[3],0,1,c-a-b1); break;fi;            
                          od;   if e2={} then  return "no solution" fi; 
                          e1:={seq(`if`(indets(i)={},NULL,i),i=op(op(e1)))};  
                          e2:={seq(`if`(indets(i)={},NULL,i),i=op(op(e2)))};   
                          B1:=solve(op(e1),b1);
                          A:=solve(subs(b1=B1,op(e2)),a);
                          L2:=subs(a=A,b1=B1,L1);
                          L2:=collect(normal(L2),{Dx,Dx^2,Dx^3},factor);
                          a1:=coeff(L,Dx^2); a2:=coeff(L2,Dx^2);
                          r:=1/3*(a1-a2);
                          Lp:=L2-DEtools[symmetric_product](L,Dx-r);
                          Lp:=collect(Lp,Dx,normal):
                          s:=[solve({coeffs(numer(Lp),{Dx,x})},{b2,c})]; 
                          if s<>[] then for i in s do 
                             B2:=op(op(i)[1])[2]; C:=op(op(i)[2])[2]; B1:=subs(c=C,B1);  S:=S union {[A,B1,B2,C,u,v]};
                          od; fi;  

     # case6: rootof(1/u) not in B and rootof(v) not in A
             elif not has({seq([i],i=l1)},{[l2[1]]}) and not has({seq([i],i=l2)},{[l1[3]]})  then lprint("case6");
                          for g in GG do
                             if g[-1]=l2[1] then e1:=MatchExp1(g[1],g[2],g[3],0,1,b1-c+1); GG:=GG minus {g}; break;fi; 
                          od; if e1={} then  return "no solution" fi; 
                          for g in GG do
                             if g[-1]=l1[3] then e2:=MatchExp1(g[1],g[2],g[3],a,b1,b1+1); break;fi;            
                          od;   if e2={} then  return "no solution" fi; 
                          e1:={seq(`if`(indets(i)={},NULL,i),i=op(op(e1)))};  
                          e2:={seq(`if`(indets(i)={},NULL,i),i=op(op(e2)))};  
                          C:=solve(op(e1),c);
                          A:=solve(op(e2),a);
                          L2:=subs(a=A,c=C,L1);
                          L2:=collect(normal(L2),{Dx,Dx^2,Dx^3},factor);
                          a1:=coeff(L,Dx^2); a2:=coeff(L2,Dx^2);
                          r:=1/3*(a1-a2);
                          Lp:=L2-DEtools[symmetric_product](L,Dx-r);
                          Lp:=collect(Lp,Dx,normal):
                          s:=[solve({coeffs(numer(Lp),{Dx,x})},{b1,b2})]; 
                          if s<>[] then for i in s do 
                             B1:=op(op(i)[1])[2]; B2:=op(op(i)[2])[2]; A:=subs(b1=B1,A); C:=subs(b1=B1,C); S:=S union {[A,B1,B2,C,u,v]};
                          od; fi; 
      # other cases
                          else     
                              if not has({seq([i],i=l2)},{[l1[2]]})  then lprint("case7");
                                 for g in GG do
                                     if g[-1]=l1[2] then e1:=MatchExp1(g[1],g[2],g[3],0,1,c-a-b1); GG:=GG minus {g}; break;fi; 
                                 od; if e1={} then  return "no solution" fi; 
                                 else 
                                 for g in GG do
                                     if g[-1]=l1[3] then e1:=MatchExp1(g[1],g[2],g[3],a,b1,1+b1); GG:=GG minus {g}; break;fi; 
                                 od; if e1={} then  return "no solution" fi; 
                              fi;    
                              if not has({seq([i],i=l1)},{[l2[2]]})  then lprint("case8");
                                 for g in GG do
                                     if g[-1]=l2[2] then e2:=MatchExp1(g[1],g[2],g[3],0,1,c-a-b2); break;fi; 
                                 od;  if e2={} then  return "no solution" fi; 
                                 else 
                                 for g in GG do
                                     if g[-1]=l2[3] then e2:=MatchExp1(g[1],g[2],g[3],a,b2,1+b2); break;fi; 
                                 od; if e2={} then  return "no solution" fi; 
                              fi;
                              e1:=simp(e1);   if e1={} then lprint("a"); return "no solution" fi; 
                              e2:=simp(e2);   if e2={} then lprint("a"); return "no solution" fi; 
                              SB1:=sol(e1,b1);
                              SB2:=sol(e2,b2);
                     for sb1 in SB1 do  for sb2 in SB2 do 
                              L2:=subs(b2=sb2,b1=sb1,L1);
                              L2:=collect(normal(L2),{Dx,Dx^2,Dx^3},factor);
                              a1:=coeff(L,Dx^2); a2:=coeff(L2,Dx^2);
                              r:=1/3*(a1-a2);
                              Lp:=L2-DEtools[symmetric_product](L,Dx-r);
                              Lp:=collect(Lp,Dx,normal):
                              s:=[solve({coeffs(numer(Lp),{Dx,x})},{a,c})]; 
                              if s<>[] then for i in s do 
                                 A:=op(op(i)[1])[2]; C:=op(op(i)[2])[2]; B1:=subs(a=A,c=C,sb1); B2:=subs(c=C,a=A,sb2); S:=S union {[A,B1,B2,C,u,v]};
                              od;fi;     od; od;                        
 fi; 

  for s in S do 
        L1:=EQN(F1(a,b1,b2,c,u,v),x); L2:=subs(a=s[1],b1=s[2],b2=s[3],c=s[4],L1);
        a1:=coeff(L,Dx^2); a2:=coeff(L2,Dx^2); r:=simplify(1/3*(a1-a2));
        SS := SS union {[op(s),r]};
  od; SS:
end:


simp:=proc(e1)
  local A, i, j, B; B:={}; if nops(e1)>1 then
for i to nops(e1) do A:={};
     for j in op(op(e1)[i]) do
           if indets(j)<>{} then A:=A union {j}; fi;
     od; B:=B union {A};
od; B;
else {seq(`if`(indets(i)={},NULL,i),i=op(op(e1)))}; fi; 
end:

sol:=proc(e1,var)
  local i,B, A,C; A:={}; B:={};
   if nops(e1)>1 then
       for i to nops(e1) do A:=solve(op(e2)[i],var); B:=B union A; od; C:={seq(op(i)[2],i=B)};
       else C:={solve(op(e1),var)};
   fi; C;     
   
end: