# 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_4a:= 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 L1,a1,a2,c,f,f1,F,i,i1,j,j1,k,k1,e0,e1,ei,Ps,Es,A0,Ai,E0,Ei,N,N1,C1,C2,Base_Field,Field_F,c_val,cand;
 
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 2 irreducible polys 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 add(points(i[1],x),i= L1) <> 4 then return "Not in Case Deg2_4a"; 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 two pair of points whose exponent differences are equal mod Z.

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

Ei:= Es[N1[1]];  Ai:= points(Ps[N1[1]],x);
if Ai <> 2 then
  if not type(evala(Es[N1[2]]-Ei),integer) and not type(evala(Es[N1[2]]+Ei),integer) then 
    return "Wrong input"; 
  fi;
fi;

# Compute f now: let's put the non ramified points at the root and pole of f:
# i.e, disappearing ramified point lies above 1.
# f:= k1*(x^2+a1*x+a0)/(x^2+b1*x+b0);  1-f:= k2*(x-c0)^2/(x^2+b1*x+b0); we do not know c0:

cand:= {};  e1:= 1/2;
for i in N do 
    C1:= N minus {i}; a1:= points(Ps[i],x);
    f:= c*Ps[i];
    if a1 = 2 then e0:= Es[i];
    else 
      for i1 in C1 while a1 <> 2 do if points(Ps[i1],x) = 1 and (type(evala(Es[i1]- Es[i]),integer) or type(evala(Es[i1]+Es[i]),integer)) then
           a1:= a1+points(Ps[i1],x);
           C1:= C1 minus {i1};
           f:= f*Ps[i1];
           e0:= min(seq(abs(Es[k]), k= {i,i1}));  fi; od; 
     fi;

     for j in C1 do         
       C2:= C1 minus {j}; a2:= points(Ps[j],x);
       f1:= f/Ps[j];
       if a2 = 2 then ei:= Es[j];
       else
         j1:= C2[1]; 
         if points(Ps[j1],x) <> 1 and not type(evala(Es[j1]- Es[j]),integer) and not type(evala(Es[j1]+Es[j]),integer) then
           return "Wrong input"; 
         fi;
         f1:= f1/Ps[j1];
         ei:= min(seq(abs(Es[k]), k= {j,j1}));   
       fi;
           c_val:= {solve({discrim(numer(1-f1),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 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 {[factor(F),[e0,e1,ei]]};
                fi;
              od;
           fi;       
       od;
   od; 
 cand;                   
end:



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