# This program finds 2F1 type solutions of degree 1, i.e; solutions of the form: 
# y(x) = exp(int(r,x))*(r_0*S(f) + r_1*diff(S(f),x)), of a second order 
# linear differential operator L where f has degree 1.
# In other words L has 3 non removable singularities.
# r,r_0,r_1 are rationalfunctions and S(f)=2F1(a,b;c|f). 


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


HypergeomDeg1:= proc(L,x,B::set)   # B:= Base Field.
local a,b,c,i,j,k,G,L1,LH,LH1,sln1,sln2,sln3,sings,RD,m,n,Base_Field,MTs,E,p,Candidates;

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

  RD:= radfield(indets(Base_Field));
 # Rewrite L and Base_Field in terms of RootOf (if any) only:
  L1:= eval(L,RD[1]);
  Base_Field:= eval(Base_Field,RD[1]);
  sings:= printsing(L1);
  # sings1:= [[seq(`if`(i[1] = infinity, i[1], solve(i[1])),i=sing)],[seq(i[2],i=sing)]];
  n:= add(nPoints(i[1],x),i = sings);
    if n <> 3 then return "We only deal with 3 singularities"; fi;
    if nops(sings) <> 3 then return "We do not handle this case here."; 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));

# Now find Mobius transformations from singularities of L to singularities of LH (0,1,infinity);
MTs:= MobiusTR(x*(x-1), mul(`if`(i[1] = infinity, 1, i[1]), i = sings), Base_Field);
Candidates:= {};
for i in MTs do 
  for j from 1 to 3 do
    if i[3][j] = 0 then p[0] := i[2][j];
    elif i[3][j] = 1 then p[1] := i[2][j];
    elif i[3][j] = infinity then p[inf] := i[2][j];
    else "wrong input"; 
    fi;
  od;
    for k in sings do 
      for m in {0,1,inf} do       
       if k[1] = p[m] or evala(eval(k[1], x=p[m])) = 0 
         then E[m]:= k[2]; break; 
         # now we get [e0, e1, ei].
       fi;
      od;
   od;
#E0:= ExpDiff[1]; E1:= ExpDiff[2]; Einf:= ExpDiff[3];
# The parameter c in 2F1(a,b;c|f) which is determined by E0 can not be 0,-1,-2,...
# We can convert such c to 1 upto equivalence; 
if denom(E[0]) = 1 and E[0] >= 1 then 
   if irem(E[0],2) = 0 or has({denom(E[1]),denom(E[inf])},2) then  E[0]:= 0;
   else 
     if E[1] < 0 then  E[0]:= 0; E[1]:= E[1]+1;
     else E[0]:= 0; E[1]:= E[1]-1;
     fi;
   fi;
fi;
a:= normal((1-E[0]-E[1]-E[inf])/2); b:= normal((1-E[0]-E[1]+E[inf])/2); c:= 1-E[0];
LH1 := transfo(eval(LH,{e0=E[0],e1=E[1],ei=E[inf]}),i[1]);
G := equiv(LH1, L1);
    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));
         Candidates := {op(Candidates),G};
od;

return sort([op(Candidates)], length)[1];

end:



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


# The following program computes the number of points (singularities)
# in a polynomial in x:

nPoints:= proc(P,x)
  if P = infinity then 1 else degree(P,x); fi;
end:


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


#This program computes the singularity structure of a differential #operator:

printsing:=proc(L)
      option remember;
	local S, sing, i,j,true_sing,v;
	S := lcoeff(primpart(L,Dx),Dx);
	S := gcd(S, lcoeff(primpart(LCLM(L, Dx) ,Dx),Dx) );
	sing := factors(S)[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:


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


#This program applies change of variable (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:

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


# This program finds a set of mobius transformations 
# (a*x+b)/(c*x+d) which carry the points {p1,p2,p3} to {q1,q2,q3}
# The program is designed to solve 2F1-solvable differential 
# operator with three singularities; which means these
# singularities can be moved to 0,1,infinity using Mobius 
# transformation. So this algorithm assumes that the points
# are not conjugates; i.e. the min poly has three (2 if one 
# point is infinity) linear factors over the base field K.

# We have the following case:

# Case 1: 1,1,1 <-----> 1,1,1.


MobiusTR:= proc(f,g,K::set)   
# f,g are the polynomials of degree 3 in x 
# (degree 2 when one root is infinity).
# K is the base field.

local  i,i1,i2,i3,f1,g1,m,p1,p2,p3,F,G,M,S1,S2,P,Q,RD,ms,eqn1,eqn2,eqn3,eqns,sol,Fctf,Fctg,Base_Field;

if nargs = 2 then
Base_Field:= indets([f,g],{radical,nonreal,RootOf});
return procname(f,g,Base_Field);
else Base_Field:= K;
fi;
RD:= radfield(Base_Field);

#Now rewrite everything in terms of RootOf (if any):
Base_Field:= eval(Base_Field,RD[1]);
f1:= eval(f,RD[1]);  g1:= eval(g,RD[1]);

if not {degree(f1,x),degree(g1,x)} subset {2,3} then
     error "Wrong input.";
fi;

  G:= factors(g,Base_Field);  F:= factors(f,Base_Field);
  S1:= {seq(`if`(has(i[1],x),i[1], NULL),i=G[2])}; 
  S2:= {seq(`if`(has(i[1],x),i[1], NULL),i=F[2])};
for i in S1 union S2 do if degree(i,x) <> 1 then return "We do not handle this case now."; fi; od;

   P:= {seq(-evala(coeff(i,x,0)/coeff(i,x,1)),i=S1)};
   Q:= {seq(-evala(coeff(i,x,0)/coeff(i,x,1)),i=S2)};
     if nops(P) = 2 then P:= P union {infinity};  fi;
     if nops(Q) = 2 then Q:= Q union {infinity};  fi;

m:= (a*x+b)/(c*x+d);

# The Mobius transformation carries elements of P to elements of Q. 
# In this case we can find all such transformations using 3 pts. in P and 
# evaluating their all possible images via m. 

  sol:={};
  p1:=P[1];  p2:=P[2];  p3:=P[3];
    for i1 in Q do  eqn1 := {IsZero(m,p1,i1)};
      for i2 in Q minus {i1} do  eqn2 := {IsZero(m,p2,i2)}; 
        i3:= op(Q minus {i1,i2});
        eqn3 := {IsZero(m,p3,i3)};   
        eqns:=eqn1 union eqn2 union eqn3;  
        sol:=  sol union{[solve(eqns),[p1,p2,p3],[i1,i2,i3]]}; 
      od;
   od;        
  M:={};    
  for i in sol do
    ms := traperror(evala(Normal(eval(m, i[1]), expanded)));
    if ms<>lasterror and has(ms,x) and IsMobiusTr(f,g,ms)  then 
       M := M union {[ms,i[2],i[3]]}; 
    fi;  
  od;  
  M;
end:


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


# This program checks,for degree 3 polynomials,if they are Mobius 
# equivalent with the given Mobius Transformation or not.

IsMobiusTr:= proc(f,g,m)  
# f,g are degree 3 polys and m is mobius transformation.
local f1,g1,IsZero;
# options trace;
g1:=g;
f1:= numer(evala(eval(f/(x-dummy)^3,x = m)));
if degree(g1,x) < degree(f1,x) then g1:= g1*denom(m); fi;
IsZero:= evala(Expand(f1- coeff(f1,x,degree(g1,x))/lcoeff(g1)*g1));
    evalb(IsZero = 0);
end:


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


# This program evaluates the mobius transformation at a given point and 
# produces a linear equation in terms of a,b,c,d:
# m:= (a*x+b)/(c*x+d)
# Let eval(m,x=p) = q .
# Then  evala(eval(m,x=p)-q)=0.

IsZero := proc(m,p,q)    # eval(m,x = p) = q.
    if p = infinity then       # when m evaluated at infinity.
      if q = infinity then  c ;    #  c = 0.
      else  evala(a-q*c) ;   #  a/c = q. 
      fi;
    else 
      if q = infinity  then evala(c*p+d); 
      else  evala((a*p+b)-(c*p+d)*q);
      fi;
    fi;
end:


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


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