# n = 4 in the following conditions:
# (a)  When f ramifies of order 2 above any one of {0,1,infty} with exponent difference odd integer/2.
# (b)  When f ramifies of order 2 above any two of {0,1,infty} with exponent difference <> odd integer/2. 


Deg2_4b:= proc(L::set,x,K)    

# L is the singularity structure of input differential operator, 
# i.e. a set of lists [pi,ei] where pi is an irreducible poly in x 
# and ei its corresponding exponent difference.
# K is the base field.

 local e0,e1,ei,i,j,k,k1,L1,Ps,A0,C1,C2,E0,Es,P,N,N1,f,f1,F,cand,c_val,Base_Field,Field_F;
 
if nargs = 2 then 
  Base_Field:= indets([args], {nonreal, RootOf, radical});
  return procname(args, Base_Field);
else Base_Field:= K;
fi;

# we are looking for the rational f which produces 4 non 
# removable singularities from 0,1,infinity.
# So, we must have at least 3 irreducible polys in this case and the polys may have at most degree 2.
  
L1:= L;
for i in L1 do if points(i[1],x) <> 1 and points(i[1],x) <> 2 then 
    return "Wrong input"; fi;
od;
if nops(L1) < 3 or add(points(i[1],x),i= L1) <> 4 then return "Not in Case Deg2_4b"; fi;

Ps:=[seq(i[1],i=L1)];  Es:= [seq(i[2],i=L1)];  N:= {seq(i,i=1..nops(L1))};

# Check exponent differences now: there should be a pair of points whose exponent differences are equal mod Z.

for i in N do 
E0:= Es[i];  A0:= points(Ps[i],x); N1:= N minus {i};
 for j in N1 while A0 <> 2 do
  if points(Ps[i],x) = 1 and (type(evala(Es[j]-E0),integer) or type(evala(Es[j]+E0),integer)) then 
    A0:= A0 + points(Ps[i],x);
  fi;
 od;
 if A0 = 2 then break; fi;
od;
if A0 <> 2 then return "Wrong input"; fi;

# Compute f now: let's put the ramified points at the root and pole of f:
# i.e, non ramified points lie above 1.
# f:= k1*(x-a0)^2/(x-b0)^2;  1-f:= k2*(x^2+c1*x+c0)/(x-b0)^2;

cand:= {};
for i in N do 
   if points(Ps[i],x) <> 1 then next; fi;
   if type(Es[i],integer) then e0:= normal(Es[i]/2); 
   else e0:= {normal(Es[i]/2), normal((Es[i]+1)/2)};
   fi;
   C1:= N minus {i}; 
   f:= c*(Ps[i])^2;
   for j in C1 do
     if points(Ps[j],x) <> 1 then next; fi;
     if type(Es[j],integer) then ei:= normal(Es[j]/2); 
     else ei:= {normal(Es[j]/2), normal((Es[j]+1)/2)};
     fi;         
     C2:= C1 minus {j};
     f1:= f/(Ps[j])^2; 
     if nops(C2) = 2 and (not type(evala(Es[C2[1]]- Es[C2[2]]),integer) and not type(evala(Es[C2[1]]+Es[C2[2]]),integer)) then next; fi;
     e1:= min(seq(abs(Es[k]), k= C2));
     P:= mul(Ps[k], k = C2);
     c_val:= {solve({coeffs(rem(numer(1-f1),P,x),x)})};
            if c_val <> {} and c_val <> {{}} then
              for k1 in c_val do F:= eval(f1,c=rhs(op(k1)));
                if max(degree(numer(F),x),degree(denom(F),x)) = 2 and type(rem(P, evala(Factor(numer(1-F))),x),numeric) then 
                Field_F:= indets(F, {RootOf, radical}) union `if`(has(F,I),{I},{});
                  if nops(Field_F minus Base_Field) <> 0 then next; fi;  
                    cand:= cand union {seq(seq([factor(F),[i1,e1,j1]],i1=e0),j1=ei)};
                fi;
              od;
           fi;     
       od;
   od; 
 cand;                   
end:



points:= proc(f,x)
if f = 1 then 1 else degree(f,x); fi;
end: