# Hypergeomdeg2:
# This program finds the solutions of second order linear differential equation in terms of 
# 2F1(a,b;c|f), if any, where f is a rational function of degree 2, i.e. the solution of the form y = exp(int(r,x))*(r0*S(x) + r1*diff(S(x),x)).
# where S(x) = 2F1(a,b;c|f), and r,r0,r1 are rational functions.
 
# The main task is to find f. We have the following choices(cases):
# Let n be the number of singularities after applying the change of variables  x--> f. Then 2 <= n <= 6.
# By Riemann-Hurwitz's formula: total amount of ramification = 2d-2 = 2.

# n = 2 when f ramifies of order 2 above any two of {0,1,infty} with exponent difference odd integer/2. That is Liouvillian Case!

# n = 3 when f ramifies of order 2 above any two of {0,1,infty} with exponent difference odd integer/2, and <> odd integer/2.
# This case will be done in another algorithm!

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

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

# n = 6 When f does not ramify above {0,1,infty}.

# This program (Deg2_6) takes the singularity structure (set of 
# lists of irreducible polys and their exp. differences) of input 
# and returns the set of f's of degree 2 which produce those 
# 6 points (the roots of those polys) from (0,1,infty).

Deg2_6:= 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,i,i1,j,j1,k,k1,b,P,Ps,Es,N,C1,C2,e0,e1,ei,Base_Field,f,f1,F,Field_F,cand,c,c_val;
 
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 6 non 
# removable singularities from 0,1,infinity.
# So, we must have at least 3 irreducible polys and the polys may have at most degree 2.
  
L1:= L; 
if nops(L1) < 3 then return "Not in Case Deg2_6"; fi;
for i in L1 do if degree(i[1],x) > 2 then return "Not in Case Deg2_6"; fi; od;
b:= add(points(i[1],x),i= L1);
  if b <> 6 then return "Not in Case Deg2_6"; fi;

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

# Compute f now:
# f:= k1*(x^2+a1*x+a0)/(x^2+b1*x+b0);  1-f:= k2*(x^2+c1*x+c0)/(x^2+b1*x+b0);
# where non of the polynomials have repeated root.

cand:= {};
  for i in N do 
        e0:= Es[i];
        C1:= N minus {i}; a1:= points(Ps[i],x);
    f:= c*Ps[i];
    if a1 = 1 then  
      for i1 in C1 while a1 <> 2 do if points(Ps[i1],x) = 1 and (type(evala(Es[i1]-e0),integer) or type(evala(Es[i1]+e0),integer)) then
           a1:= a1+points(Ps[i1],x);
           C1:= C1 minus {i1};
           f:= f*Ps[i1]; fi; od;
        if a1 <> 2 then next; fi; 
        fi;
      for j in C1 do         
         ei:= Es[j]; C2:= C1 minus {j}; a2:= points(Ps[j],x);
         f1:= f/Ps[j];
         if a2 = 1 then 
           for j1 in C2 while a2 <> 2 do if points(Ps[j1],x) = 1 and (type(evala(Es[j1]-ei),integer) or type(evala(Es[j1]+ei),integer)) then
           a2:= a2+ points(Ps[j1],x);
           C2:= C2 minus {j1};
           f1:= f1/Ps[j1]; fi; od;
           if a2 <> 2 then next; fi;
           fi;
           e1:= Es[C2[1]];
           if nops(C2) = 2 then 
             if not type(evala(Es[C2[2]]-e1),integer) and not type(evala(Es[C2[2]]+e1),integer) then
                return "Wrong input"; 
             fi;
           fi;
           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 {[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: