# Liouvillian solutions for order 3, imprimitive case.
#
# Main features that are implemented:
#
#  1).  The theorem in [HRUW] that gives the Riccati polynomial.
#
#  2).  The recurrence relations mentioned in Section 4.4 in [HW].
#
#  3).  These recurrence relations are combined with denominator bounds in order
#       to get some preliminary linear equations. This way we get good efficiency
#       without using the formal solutions method given in [HW], which is harder
#       to implement.
#
#  4).  With tricks specific to the order 3 imprimitive case we can give sharper
#       denominator bounds than the more general bounds given in [HW].
#
#  5).  We're searching for a semi-invariant, i.e. an exponential function
#       times something rational. This exponential function is:
#           sqrt(discriminant of Riccati poly).
#       With local analysis we can generally determine this discriminant
#       up to a square (in rare occasions more than 1 case needs to be
#       checked) by looking at where the branch points of index 2 could be,
#       which means that in the order 3 imprimitive situation, computing
#       semi-invariants is not more expensive than computing invariants (the
#       latter is otherwise considered an easier problem) because there's
#       rarely more than one type to consider.  The local analysis is similar
#       to that in DEtools[kovacicsols] and uses some unpublished tricks.
#
#  6).  To determine a semi-invariant that factors, see [SU], we use just 1 Brill
#       equation, but a very special one, namely one that must either give the
#       zero equation or the correct equation (no spurious solutions) for the parameter.
#       The zero equation can be avoided by a linear transformation of the coordinates,
#       and thus this 1 Brill equation is the only one that we need.
#
#
# Mark van Hoeij, June 2003.
#
#
#
# Citations:
#
#  [HW] M. van Hoeij, J-A. Weil, "An algorithm for computing invariants
#  of differential Galois groups" J. Pure Appl. Algebra, 117&118, 353-379 (1997)
#
#  [HRUW] M. van Hoeij, J-F. Ragot, F. Ulmer, J-A. Weil "Liouvillian solutions of
#  linear differential equations of order three and higher" J. Symb. Comput., 28,
#  589-609 (1999)
#
#  [SU] M. Singer, F. Ulmer "Linear Differential Equations and Products of Linear
#  Forms. J. of Pure and Applied Algebra, 117-118,  549-563  (1997).


macro(
	solve_linear = SolveTools:-Linear, # Use in Maple >7.
	sort_poly	= PolynomialTools:-Sort,
	g_factors       =`DEtools/kovacicsols/factors`,
	eqns = `order3/eqns`,
	find_var = `order3/find_var`
):

eqns := proc(m30, m12, vars, a0, a1, a2, X1, X2, X3, x, mden, ii, jj, nth)
	local m03,m21,m21d,m20,m12d,m11,m02,m01,m10,m00,vergl1,vergl2,i,sol,d;

	# The variable name "mij" refers to:  y^i * (y')^j * (y'')^(3-i-j)
	# The formulas we use to determine the "mij" in this algorithm from m30 and m12
	# are the recurrence relations mentioned in Section 4.4 in [HW].
	# In essence, a "recurrence relation" between the "mij" is simply an
	# equation in the symmetric power system (viewed as 10 equations in 10 unknowns).
	# Out of the 10 equations, we use eight as "recurrence relations", i.e. we use
	# those eight to express eight "mij" in terms of m30, m12.
	# The two remaining equations in the symmetric power system are then necessary
	# and sufficient conditions for m30, m12.
	#
	# As input we require the monomials m30 = y^3 and m12 = y * (y')^2
	# These monomials are rational functions of the form:
	#   "a known rational function" times
	#   "an unknown polynomial for which we have a degree-bound"
	#
	# This unknown polynomial contains a number of unknown constants
	# (if we use the formal solutions method presented in [HW]
	#  then the number of unknown constants will be very small, and if we simply use
	#  an unknown for each coefficient then the number of unknown constants appearing
	#  in m30 is degree-bound + 1,  and likewise for m12).  These unknowns,
	# however many there may be, must be given in the input "vars".
	#
	m21 := Normalizer( 1/3*diff(m30,x) );
	m21d := Normalizer(diff(m21,x));
	m20 := Normalizer( m21d-2*m12 );
	m03 := Normalizer( -a0*m30 - a1*m21 - a2*m20 + 3*diff(m12,x) - diff(m21d,x) );
      if nth < 1 then
	# get some preliminary linear equations using a denominator bound.
	# This helps reduce the number of unknowns, which reduces expression size,
	# which in turn improves the running time.
	#
	# Note that the code between "if nth < .. then" and  "fi;" is only a
	# speedup, i.e., this procedure still works if you remove these if..fi's.
	#
	d:=numer(Normalizer(denom(m03)/eval(subs(ii=0,jj=3,mden))));
	if has(d,x) then
		sol := solve_linear({coeffs(collect(evala(Rem(numer(m03),d,x)),x),x)},vars);
		d := map(Normalizer,subs(sol, [m30,m12]));
		return procname(op(d), indets(d) intersect vars, args[4..-2],1)
	fi
      fi;
	m12d := Normalizer( diff(m12,x) );
	m11 := Normalizer( 1/2*m12d - 1/2*m03 );
	m02 := Normalizer( 1/3*diff(m03,x) );
	m01 := Normalizer( 1/2 * (a0*m12 + a1*m03 + a2*m02 + diff(m02,x) ));
	m10 := Normalizer( 1/2*diff(m12d,x) - 5/2*m02 + a0*m21 + a1*m12 + a2*m11);
      if nth < 2 then
	d:=numer(Normalizer(denom(m01)/eval(subs(ii=0,jj=1,mden))));
	if has(d,x) then
		sol := {coeffs(collect(evala(Rem(numer(m01),d,x)),x),x)}
	else
		sol := {}
	fi;
	d:=numer(Normalizer(denom(m10)/eval(subs(ii=1,jj=0,mden))));
	if has(d,x) then
		sol := sol union {coeffs(collect(evala(Rem(numer(m10),d,x)),x),x)}
	fi;
	if sol<>{} then
		sol := solve_linear(sol,vars);
		d := map(Normalizer,subs(sol, [m30,m12]));
		return procname(op(d), indets(d) intersect vars, args[4..-2],2)
	fi
      fi;
	m00 := Normalizer( 2*a0*m11 + 2*a1*m02 + 2*a2*m01 + diff(m01,x) );
      if nth < 3 then
	d:=numer(Normalizer(denom(m00)/eval(subs(ii=0,jj=0,mden))));
	if has(d,x) then
		sol := solve_linear({coeffs(collect(evala(Rem(numer(m00),d,x)),x),x)},vars);
		d := map(Normalizer,subs(sol, [m30,m12]));
		return procname(op(d), indets(d) intersect vars, args[4..-2],3)
	fi
      fi;
	vergl1 := Normalizer( -2*a0*m20 -2*a1*m11 -2*a2*m10 + m01 - diff(m10,x) );
	vergl2 := Normalizer( -3*a0*m10 -3*a1*m01 -3*a2*m00 - diff(m00,x) );

	sol := solve_linear({seq(coeffs(collect(numer(i),x),x),i=[vergl1,vergl2])}, vars);
	#
	# This solved vergl1 = 0, vergl2 = 0. Then the following vector will
	# be a solution of the symmetric power system:
	#   [m30,m21,m12,m03,m20,m11,m02,m10,m01,m00]
	#
	# Translate this solution to a polynomial using: 
	#    collect( (X1*y + X2*yp + X3*ypp)^3, {y,yp,ypp},distributed);
	# and replace each monomial y^i*yp^j*ypp^(3-i-j) by the monomial mij,
	# which leads to the polynomial:

	subs(sol, m30*X1^3 + 3*m21*X1^2*X2 + 3*m12*X1*X2^2 + m03*X2^3 + 3*m20*X1^2*X3  
	 + 6*m11*X1*X2*X3 + 3*m02*X2^2*X3 + 3*m10*X1*X3^2 + 3*m01*X2*X3^2 + m00*X3^3 )

	# We'll then determine the Riccati poly with the theorem in [HRUW]
	# This theorem simply says:   subs(X2 = -1, X3 = 0, polynomial).
end:


# Input: L must be an irreducible third order operator.
# Output: A Liouvillian solution iff L has imprimitive Galois group.
#
order3 := proc(L, ext::set)
	local T,Dx,x,pos_set,sing,cm,dens,n3,n2,q,p,dp,qq,g,i,j,S,ld,gg,k
	 ,ram,Nbr,U,N,B0,TR,LL,V,SV,sv,nb,b0,L2,m30,m12,m_ijk,b12,var,Pol,
	 Pol_over_C,X1,X2,X3,ii,jj,kk;
	if nargs = 1 then
		return procname(L, indets(A,{'nonreal', 'radical', 'algext'}))
	fi;
	Dx,x := op(_Envdiffopdomain);
	Normalizer:=`if`(ext={},normal,evala);
	if not type(L,polynom(ratpoly(anything,x),Dx)) or degree(L,Dx)<>3 then
		error "wrong input"
	elif lcoeff(L,Dx)<>1 then
		return procname(collect(L/lcoeff(L,Dx),Dx,Normalizer), ext)
	fi;
	pos_set:={_IMPR_}; # So far only the imprimitive case implemented.
	sing := sort_poly([infinity,seq(i[1],i=g_factors(
	 denom(Normalizer(L)),x,[op(ext)]))],x);
	cm:=NULL;
	dens:={}: # set of denominators. Include x when not rational
	n3 := 0; # number of singularities with branch index 3.
	n2 := 0; # number of singularities with branch index 2.
for q in sing while pos_set<>{} do
	if q=infinity then
		p:=q;
		dp[p]:=1
	else
		p:=RootOf(q,x);
		dp[p]:=degree(q,x);
		qq[p]:=q
	fi;
	cm:=cm,p;
	g[p]:=DEtools['gen_exp'](L,T,x=p);

	# If gen. exp. has multiplicity > 1 then stop.
	for i in g[p] do j:=i[1..-2]; if nops(j)>nops({op(j)}) then
		return NULL
	fi od;

	if nops(g[p])=2 and nops(g[p][1])+nops(g[p][2])=5 then
		# Means that two (not three) gen.exp differ by an integer.
		# Replace the "bigger" one by the smaller one:
		g[p] := [seq(seq([i[1]], j=1..nops(i)-1),i=g[p])]
	fi;

	if nops(g[p])=1 then
		if nops(g[p][1])=2 then
			# ramification or algebraic extension.
			if lhs(g[p][1][2])=T then
				# algebraic extension
				S[p] := evala(Trace(g[p][1][1],
				 ext union indets(g[p][1][1],RootOf),
				 ext))/3;
				ld[p] := ldegree(collect(g[p][1][1]-S[p],T,evala),T)
			else
				# ramification
				n3 := n3 + dp[p];
				ld[p] := ldegree(g[p][1][1],T);
				j := T^3/lhs(g[p][1][2]);
				S[p] := add(Normalizer(coeff(g[p][1][1],T,3*i)*j^i)*T^i,
				 i=ld[p]..0);
				while ld[p] mod 3 = 0 or coeff(g[p][1][1],T,ld[p]) = 0 do
					ld[p] := ld[p] + 1
				od;
				ld[p] := ld[p]/3
			fi
		elif nops(g[p][1])=4 then
			# if p<>infinity and Normalizer(g[p][1][3]-g[p][1][1])=2 then
			# 	return NULL
			# fi;
			# (semi)-apparant singularity
			S[p] := g[p][1][1];
			if p=infinity then
				ld[p] := 0
			else
				ld[p] := 4 # Indicates apparant singularity.
			fi
		else
			error "should not happen"
		fi
	elif nops(g[p])=3 then
		S[p] := collect(add(i[1],i=g[p])/3,T,Normalizer);
		gg[p] := [seq(collect(i[1]-S[p],T,Normalizer),i=g[p])];
		ld[p] := min(op(map(ldegree,gg[p],T)));
		k := 0;
		for i from 1 to 2 do for j from i+1 to 3 do if k=0 then
			k := 2*Normalizer(gg[p][i]-gg[p][j]);
			if type(k,integer) and type(k,odd) then
				# possible ramification 2
				n2 := n2 + dp[p];
				gg[p] := abs(k/2);
				S[p] := S[p] - gg[p]/3;
				ram[p] := 2; # Indicates a possible, not certain, ramification 2.
				if nops({op(g[p])})=2 and p<>infinity
				and min(seq(Normalizer(i[1]-S[p]),i=g[p]))=0 then
					ld[p] := 2
				fi
			else
				k:=0
			fi
		fi od od;
		k := {seq(3*i,i=gg[p])};
		if type(k,set(integer)) and nops(k mod 3)=3 then
			# Possible ramification 3
			n3 := n3 + dp[p];
			j := min(op(gg[p]));
			S[p] := S[p] + j;
			gg[p] := sort([seq(i-j,i=gg[p])])[2];
			if p<>infinity then
				ld[p] := 3 # Indicates a possible, not certain, ramification 3.
			fi
		fi
	elif nops(g[p])=2 then
		if indets(g[p][2],RootOf) minus ext <> {}
		or lhs(g[p][2][2])<>T then
			g[p] := [ g[p][2], g[p][1] ]
		fi;
		if indets(g[p][1],RootOf) minus ext <> {} then
			# algebraic extension
			S[p] := 1/3 * (g[p][2][1] + evala(Trace(g[p][1][1],
			 indets(g[p][1],RootOf) union ext, ext)) );
			ld[p] := min(seq(ldegree(collect(i[1]-S[p],T,evala),T),i=g[p]))
		elif lhs(g[p][1][2])<>T then
			# ramification
			n2 := n2 + dp[p];
			# Not necessarily optimal:
			ld[p] := min(ldegree(g[p][1][1],T)/2, ldegree(g[p][2][1],T));
			j := T^2/lhs(g[p][1][2]);
			S[p] := (add(Normalizer(coeff(g[p][1][1],T,2*i)*j^i)*2*T^i,
			 i=ceil(ld[p])..0) + g[p][2][1])/3
		else
			error "should not happen"
		fi
	fi
od;
	Nbr := 2*n3 + n2;
	if Nbr < 4 then
		# Not enough branch points.
		return NULL
	fi;
	cm := {cm};
	p := infinity;
	# constant_term minus gen_exp at infinity:
	U[p] := collect(subs(T = 1/x, coeff(S[p],T,0)-S[p] )/x ,x);
	N[p] := coeff(S[p],T,0);
	for p in cm do if p<>infinity then
		U[p] := subs(T = x-p, S[p]/T);
		N[p] := coeff(S[p],T,0);
		if dp[p]>1 then
			U[p], N[p] := seq(evala(Trace(i, ext union {p}, ext)), i = [U[p], N[p]])
		fi
	fi od;
	B0 := -3*Normalizer(add(N[p],p=cm));  # 3* because 3'rd symmetric power
	if not type(2*B0,integer) then
		return NULL
	fi;
	TR := Normalizer(add(U[p],p=cm));
	LL := DEtools['symmetric_product'](Dx+TR, L);
	V := {seq(`if`(ram[p] = 2, p, NULL),p=cm)};
	for p in V do
		if p=infinity then U[p] := 0
		else
			U[p] := gg[p]/(x-p);
			if dp[p]>1 then
				U[p] := evala(Trace(U[p], ext union {p}, ext))
			fi
		fi
	od;
	SV := combinat['subsets'](V);
	while not SV[finished] do
		sv := SV[nextvalue]();
		nb := Nbr - add(dp[p],p=sv);
		b0 := B0 - add(dp[p]*gg[p], p = sv);
		if nb<4 or type(nb,odd) or not type(b0,integer) then
			next
		fi;
		L2 := DEtools['symmetric_product'](Dx + add(U[p],p=sv)/3, LL);
		# lprint(L2);
		m30 := add(c30[i]*x^i,i=0..b0);
		m12 := mul(qq[p]^`if`(ld[p] <= 0, 2*(1-ld[p]),
				 `if`(ld[p] = 2, `if`(member(p,sv),2,`if`(gg[p]=1/2, 1,0)),
				 `if`(ld[p] = 3, max(0,floor(2*(1-gg[p])))
					, 0))), p=cm minus {infinity});
		m_ijk := mul(qq[p]^`if`(ld[p] <= 0, (jj+2*kk)*(1-ld[p]),
			`if`(ld[p] = 2 or ld[p]=3, `if`(member(p,sv),jj+2*kk,
			 floor( jj*max(0, 1-gg[p]) + kk*max(0,2-gg[p]) ))
			, 0)), p=cm minus {infinity});
		if Normalizer(eval(subs(jj=2,kk=0,m_ijk))/m12) <> 1 then
			error "should not happen"
		fi;
		b12 := b0 + degree(m12,x) - ceil(2*(1+ld[infinity]));
		m12 := add(c12[i]*x^i,i=0..b12)/m12;
		if lcoeff(L2,Dx)<>1 then L2 := collect(L2/lcoeff(L,Dx),Dx) fi;
		var := indets([m30,m12]) minus indets([L,x]);
		Pol := eqns(m30,m12,var, seq(coeff(L2,Dx,i),i=0..2), X1,X2,X3, x,
		 subs(kk=3-ii-jj, m_ijk),ii,jj,0);
		if Pol<>0 then
			var := indets(Pol) intersect var;
			if nops(var)=1 then
				Pol := subs(var[1]=1,Pol)
			elif nops(var)=2 then
				i:=2;
				# Using an idea by Bronstein to replace x by a non-singular point
				# to get the polynomial form (over the constants) of the
				# invariant:
				while Normalizer(subs(x=i,denom(L2)))=0 do i:=i+1 od;
				Pol_over_C := subs(x=i,Pol);

				# Use a Brill equation:
				Pol := subs(find_var(Pol_over_C,X1,X2,X3,var,ext), Pol)
			else
				return NULL # Should not happen if L is irreducible
			fi;

			# Now use the main theorem in [HRUW] to get the Riccati Polynomial:
			Pol := subs(X2=-1,X3=0,Pol);

			TR := Normalizer(TR+add(U[p],p=sv)/3);
			# Make it monic, shift by TR, return the answer.

			Pol := collect(subs(X1=X1-TR,Pol/lcoeff(Pol,X1)), X1, evala);
			return exp(Int( RootOf(Pol,X1) ,x))
		fi;
	od;
	NULL
end:


find_var := proc(Pol,X1,X2,X3,var,ext)
	local res,P,r,i,j,b;
	res := 0;
	P := Pol;
	r := rand(-10^3..10^3);
	while res=0 do
		P := collect(P, {X1,X2,X3}, distributed);
		for i from 0 to 3 do for j from 0 to 3-i do
			b[i,j] := coeff(coeff(coeff(P,X1,i),X2,j),X3,3-i-j)
		od od;
		res := collect(
		# Brill equation.
		#
		# Note: this equation is of degree 4, which is the correct degree because there
		# are exactly 4 correct parameter values, see also the paper by
		# Hendriks and van der Put, JSC, 19, 1995.  By looking at the factorable semi-invariants
		# of the possible groups here, one can calculate that the following equation
		# must be a Mobius transformation of var[1] * (var[2]^3 - var[1]), which implies
		# that the Galois group over the constants of this equation must be a subgroup
		# of D_4 in general,  and must be a subgroup of C_4 if sqrt(-3) in constants.
		#
		 b[3,0]*b[2,1]*b[0,2]^2-b[3,0]*b[1,1]*b[0,2]*b[1,2]+b[3,0]*b[0,1]*b[1,2]^2
		 -3*b[3,0]*b[0,1]*b[2,1]*b[0,3]+b[0,3]*b[2,0]*b[2,1]*b[1,1]-b[2,1]^2*b[0,3]*b[1,0]
		 -b[0,3]*b[2,0]^2*b[1,2]+3*b[3,0]*b[1,0]*b[0,3]*b[1,2], var, Normalizer);

		# If res=0 (of course this never happens, but we include this case because we want
		# the code to be complete) then we'll use this change of coordinates:
		P := subs(X1 = X1+r()*X2+r()*X3, X2 = X2+r()*X3+r()*X1, X3 = X3+r()*X1+r()*X2, Pol)
	od;
	if coeff(res,var[1],0)=0 then
		var[1]=0, var[2]=1
	else
		res := g_factors( subs(var[1]=1, res), var[2], [op(ext)] );
		# Use sort_poly to take a factor of minimal degree.
		var[1]=1, var[2]=RootOf(sort_poly([seq(i[1],i=res)],var[2])[1],var[2])
	fi
end: