# This program finds hypergeometric solution of degree 2 of a second order linear differential operator;
# i.e. solution of the form: y(x) = exp(int(r,x))*(r_0*S(f) + r_1*diff(S(f),x)),
# where r,r_0,r_1 are rational functions and S(f)=2F1(a,b;c|f) with degree(f,x) = 2:



_Envdiffopdomain := [Dx,x]:
with(DEtools):


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

# Here, L is a second order differential operator and K is the specified (base) field.

local a,b,c,d,i,G,L1,LH,LH1,Cand,Base_Field,Sing_L1,soln;

 if nargs = 2 then 
     Base_Field:= indets([args], {RootOf, nonreal, radical}) union `if`(has([args],I),{I},{});
     return procname(args, Base_Field);
 else Base_Field:= K;
 fi;

LH:= Dx^2-(-2*x+x*e0+x*e1+1-e0)*Dx/(x*(x-1))+(1/4)*(e0+e1-ei-1)*(e0+e1+ei-1)/(x*(x-1));
L1:= L;
Sing_L1:= eval(printsing(L1,Base_Field), infinity = 1);  
d:= add(points(i[1],x), i=Sing_L1);
 
 if d < 3 or d > 6 then return "Does not have a degree 2 solver"; fi;

if d = 3 then Cand:= Deg2_3(Sing_L1,x,Base_Field);
elif d = 4 then Cand:= Deg2_4(Sing_L1,x,Base_Field);
elif d = 5 then Cand:= Deg2_5(Sing_L1,x,Base_Field);
else Cand:= Deg2_6(Sing_L1,x,Base_Field);
fi;

soln:= NULL;
 if Cand <> {} and not(type(Cand,string)) then 
  for i in Cand do
      if normal(add(makemin(i[2][j]),j=1..3)) >=1 then  lprint("This selection leads to Liouvillian solution"); next; fi;
      # We don't want c in 2F1(a,b;c|f) to be 0,-1,-2,... We can convert such c to 1 upto equivalence; 
       if denom(i[2][1]) = 1 and i[2][1] >= 1 then 
          if irem(i[2][1],2) = 0 or has({denom(i[2][2]),denom(i[2][3])},2) then i[2][1] := 0;
          else if i[2][2] < 0 then  i[2][1]:= 0; i[2][2]:= i[2][2]+1; else i[2][1]:= 0; i[2][2]:= i[2][2]-1; fi;
          fi;
       fi;          
    LH1:= transfo(eval(LH,{e0=i[2][1],e1=i[2][2],ei=i[2][3]}),i[1]);
     a:=normal((1-i[2][1]-i[2][2]-i[2][3])/2); b:= normal((1-i[2][1]-i[2][2]+i[2][3])/2); c:=normal(1-i[2][1]); 
       G := equiv(LH1, L);
        if G = 0 then next; fi;
	  G := eval(diffop2de(G, y(x)),  y(x) = hypergeom([a,b],[c],i[1]));
          if not(has(indets(G, function),hypergeom)) then next; fi;
	  G := collect(G, select(has,indets(G, function),hypergeom));
          soln := soln,G; od; fi;
	 if soln <> NULL then {sort([op({soln})], length)[1]}; else {}; fi;
end:

###################################################################################################################

# The main task is to find f. We have the following choices(cases):
# Let d be the number of singularities after applying the change of variables  x--> f. Then 2 <= d <= 6.
# By Riemann-Hurwitz's formula: total amount of ramification = 2n-2 = 2. (n = 2).

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

# d = 3 
# Deg2_3: When f ramifies of order 2 above any two of {0,1,infty} with exponent difference odd integer/2, and <> odd integer/2.

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

# d = 5 :
# Deg2_5: When f ramifies of order 2 above any one of {0,1,infty} with exponent difference <> odd integer/2.

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

# Each case (Deg2_3..Deg2_6) below take the singularity structure L of the input differential operator
# and return the candidate f and [e0,e1,ei].

###########################################################################################################################

# d = 3 when f ramifies of order 2 above any two of {0,1,infty} with exponent difference odd integer/2, and <> odd integer/2.

Deg2_3:= proc(L::set,x,Base_Field)    

local L1,a2,i,i1,j,j1,k,k1,P,Ps,Es,N,C1,C2,e0,e1,ei,f,f1,F,Field_F,cand,c,c_val;

# we are looking for the rational f which produces 3 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.
# at least one of them must have degree 1:
  
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) <> 3 then return "Not in Case Deg2_3"; fi;

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

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

cand:= {};  e1:= 1/2;
  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;
    f:= c*(Ps[i])^2;
    C1:= N minus {i};
    if nops(C1) = 1 then 
      P:= Ps[C1[1]]; ei:= Es[C1[1]];
    else 
      if type(evala(Es[C1[1]]- Es[C1[2]]),integer) or type(evala(Es[C1[1]]+Es[C1[2]]),integer) then
          P:= mul(Ps[j], j= C1); ei:= min(seq(abs(Es[j]), j= C1));
      else next;
      fi;
    fi;
     f1:= f/P;
     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 {seq([factor(F),[k2,e1,ei]], k2 = e0)};
                fi;
              od;
           fi;       
       od; 
 cand;                   
end:


##########################################################################################################

# d = 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 two of {0,1,infty} with exponent difference <> odd integer/2.

Deg2_4:= proc(L::set,x,Base_Field)

local sol1,sol2;

sol1:= Deg2_4a(L,x,Base_Field);   sol2:= Deg2_4b(L,x,Base_Field);

return `if`(type(sol1, string), {}, sol1) union `if`(type(sol2, string), {}, sol2);

end:


##########################################################################################################

# d = 4:
# (a)  When f ramifies of order 2 above any one of {0,1,infty} with exponent difference odd integer/2.

Deg2_4a:= proc(L::set,x,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,Field_F,c_val,cand;

# 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:


#################################################################################################################################


# d = 4 in the following conditions:
# (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,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,Field_F;

# 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:


####################################################################################################################################


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

Deg2_5:= proc(L::set,x,Base_Field)    

local L1,a2,i,i1,j,j1,k,k1,P,Ps,Es,N,C1,C2,e0,e1,ei,f,f1,F,Field_F,cand,c,c_val;
 
# we are looking for the rational f which produces 5 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.
# at least one of them must have degree 1:
  
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) <> 5 then return "Not in Case Deg2_5"; fi;

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

# Compute f now: let's put the ramified point at the root of f:
# f:= k1*(x-a0)^2/(x^2+b1*x+b0);  1-f:= k2*(x^2+c1*x+c0)/(x^2+b1*x+b0);
# where other polynomials have non repeated roots.

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         
         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 next; 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 {seq([factor(F),[k2,e1,ei]], k2 = e0)};
                fi;
              od;
           fi;       
       od; 
  od;
 cand;                   
end:


######################################################################################################################################


# d = 6:
# When f does not ramify above 0,1,infinity.

Deg2_6:= proc(L::set,x,Base_Field)    

local L1,a1,a2,i,i1,j,j1,k,k1,b,P,Ps,Es,N,C1,C2,e0,e1,ei,f,f1,F,Field_F,cand,c,c_val;

# 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:

#######################################################################################################################

# Compute places of f:

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

########################################################################################################################

# Compute the singularity structure of a differential operator:

printsing:=proc(L, B::set)

# Here, L is a second order differential operator and B is the base field or some specified field.

	#global Dx, x;
	local S, sing, i,j,Field_L, B1, true_sing,v;

	S := lcoeff(primpart(L,Dx),Dx);
	S := gcd(S, lcoeff(primpart(LCLM(L, Dx) ,Dx),Dx) );
	sing := factors(S, B)[2];
        sing:=[seq(i[1], i=sing), infinity];
        true_sing:=NULL;
        for i in sing do
		j := `if`(i=infinity, i, RootOf(i,x));
                v := gen_exp(L,T,x = j );
                if nops(v)=2 or v[1][1]=v[1][2] or 
                   (not type(v[1][1]-v[1][2], integer)) or formal_sol(L,`has logarithm?`,x=j) then
                        true_sing := true_sing, [i, `if`(nops(v)=1, v[1][2]-v[1][1], v[2][1]-v[1][1]) ]
                fi;
        od;
	{true_sing}
end:


#############################################################################################################################

# Apply the change of variables x ---> f to a differential operator:

transfo:=proc(L,a)
#option trace;
local i,f;        
 
f:=add(mult( subs(x=a,coeff(L,Dx,i)),(1/diff(a,x) * Dx)$i),i=0..degree(L,Dx));        
sort(collect(f/lcoeff(f,Dx),Dx,factor),Dx)
end:

##########################################################################################################################

# Compute the projective equivalence between two differential operators:

# This program computes the projective equivalence:  y ----> exp(int(r,x))# (r_0 y + r_1 diff(y,x)) between two operators:
# Example to show the syntax:
#
# _Envdiffopdomain := [Dx,x]:  # (see help page of DEtools[mult] for this)
#
# L1 := x^2*Dx^2-4*x^6+4*x^4-3/4;
# L2 := Dx^2 + 2-10*x+4*x^2-4*x^4;
# equiv(L1,L2);
# --> Gives the map from the solution space of L1 to the solution space
#  of L2, showing effectively that solving L1 is equivalent to solving L2.
#     To get the inverse map, call equiv(L2,L1).

# Checks if M2 can be solved in terms of M1, and if so, finds the map.

equiv:=proc(M1,M2,d)
	local DF,x,T1,L1,T2,L2,C0,D0,S,C1,D1,dd,v,v1;
  if nargs>2 then
	if type(d,list(name)) and nops(d)=2 then
	  _Envdiffopdomain:=d
	elif type(d,function) then
	  _Envdiffopdomain:=[DF,x];
	  return `DEtools/diffop2de`(procname(seq(
		  `DEtools/de2diffop`(args[i],d),i=1..2)),d)
	else
		error "wrong number or type of arguments"
	fi
  elif not assigned(_Envdiffopdomain) then
	error "domain is not assigned"
  fi;
	DF,x := op(_Envdiffopdomain);
	if not type([M1,M2],list(polynom(ratpoly(anything,x),DF))) then
		error "wrong arguments or dependent variable not specified"
	elif degree(M1,DF)<>2 or degree(M2,DF)<>2 then
		error "case not handled yet"
	fi;

        if indets([args],{RootOf,radical,nonreal}) = {} then
                Normalizer := normal
        else
                Normalizer := evala
        fi;

	T1:=Normalizer(coeff(M1,DF,1)/coeff(M1,DF,2)/2);
	L1:=`DEtools/symmetric_product`(M1,DF-T1);
	T2:=Normalizer(coeff(M2,DF,1)/coeff(M2,DF,2)/2);
	L2:=`DEtools/symmetric_product`(M2,DF-T2);
	C0:=Normalizer(coeff(L1,DF,0)/coeff(L1,DF,2));
	D0:=Normalizer(coeff(L2,DF,0)/coeff(L2,DF,2)-C0);
  if D0=0 then
	S:=1
  else
	C1:=Normalizer(diff(C0,x));
	D1:=Normalizer(diff(D0,x));
	dd:=Normalizer(D1/D0);
	v:=map(Normalizer,[2*diff(C1,x)+diff(D1,x)-2*dd*C1-D1*dd+D0^2,
	 6*C1+D1-4*dd*C0, 2*(D0+2*C0), -dd, 1]);
	v:=`DEtools/Expsols`(v,0,x,`use Int`,`no algext`,radical,denom,2);
	if nops(v)=0 then
		return 0
	elif nops(v)>1 then
		userinfo(1,'dsolve',"input is reducible");
		return 0
	else
		v:=v[1];
		v1:=Normalizer(diff(v,x));
		S:=collect(v*DF +
		  (diff(v1,x,x)+(4*C0+5*D0)*v1+v*(2*C1+D1))/2/D0-2*v1
		,[exp,DF],Normalizer)
	fi
  fi;
	S:=combine(collect(exp(Int(-T2,x))*`DEtools/mult`(
	 S,exp(Int(T1,x))),[exp,DF],Normalizer),exp);
	v:=seq(`if`(op(0,i)=exp,i,NULL),i=indets(S,function));
	if nops([v])<>1 then
		error "exp not collected"
	fi;
	S := subs(v=`DEtools/kovacicsols/e_int`(Normalizer(diff(op(v),x)),x),S);
    if has(S,RootOf) and not has(S,DF) then add(`DEtools/kovacicsols/combi`(normal(coeff(S,DF,i)),x)*DF^i, i=0..degree(S,DF));
    else  collect(subs(v=`DEtools/kovacicsols/e_int`(Normalizer(diff(op(v),x)),x),S),DF,Normalizer);
    fi;	
end:


########################################################################################################################

# The following program returns the normal form (1/denom(ei) or 0) of exponent difference
# which helps to check whether the GHE has Liouvillian solution or not:

makemin:= proc(a)
 if denom(a) = 1 then return 0;  #logarithmic case
 else return 1/denom(a); fi;
end: