# To use this program in Maple, the file "Conic" must also be read
# into Maple. If your Maple version is older than Maple 10, then the
# file "pyth_LLL" needs to be read into Maple as well.

# Find a transformation that sends a third order differential operator
# to an operator of the form symmetric_power(L2, 2) for some second
# order operator L2, if such transformation exist over the rationals.
#
# Author: Mark van Hoeij.

macro(
	symmetric_power = DEtools['symmetric_power'],
	diffop2de = DEtools['diffop2de'],
	symmetric_product = DEtools['symmetric_product'],
	ExpSolsOverQ = `DEtools/Expsols`,
	LCLM = DEtools['LCLM'],
	mult = DEtools['mult'],
	rightdivision = DEtools['rightdivision'],
	SolveLinear = SolveTools:-Linear
):
if not assigned(infolevel[dsolve]) then
	infolevel[dsolve] := 1
fi:

# Input:
#   *)   L: a third order irreducible differential operator
#   *)   (optional) domain: a list containing two names.
#
# If domain = [Dx, x] then that means that the independent variable
# is x, and that Dx is d/dx. For more details on this syntax see
# the help page of DEtools[mult]. If the variable _Envdiffopdomain 
# has been assigned then the argument domain may be omitted.
#
# L must have rational coefficients (no algebraic numbers).
#
# The output will be, if it exists, a list [L2, R] where L2, R have
# rational coefficients, L2 has order 2, and R sends the solutions of
# symmetric_power(L2, 2) to the solutions of L.  If such L2, R do not
# exist then the output will be FAIL.
#
ReduceOrder := proc(L3, domain)
	local Dx,x,L6,a,b,r,L,a0,a1,a2,c,cn,cd,i,j,C,s0,s1,pt,UsedPt,
	 b0,b1,b2,R,LR,S,eq,L2;

	# Type checking and preprocessing.
	#
	if nargs>1 and type(domain,list(name)) and nops(domain)=2 then
		_Envdiffopdomain := domain
	elif not type(_Envdiffopdomain,list(name)) then
		error "differential algebra not specified"
	fi;
	Dx, x := op(_Envdiffopdomain);
	if not type(L3, polynom(ratpoly(rational),Dx)) then
		error "Only the rational coefficients case is implemented"
	elif degree(L3,Dx)<>3 then
		error "Only third order is implemented"
	elif lcoeff(L3,Dx)<>1 then
		# Make input monic
		return procname(collect(L3/lcoeff(L3,Dx),Dx,normal))
	fi;
	if not assigned(_Env_DF_sq) then _Env_DF_sq := table([x, denom(L3)]) fi;

	# Compute the symmetric square of L3, which is of order 6
	# or of order 5 if the problem is trivial.
	#
	L6 := symmetric_power(L3, 2);
	if degree(L6, Dx) = 5 then
		userinfo(1,'dsolve',"input is already a symmetric square");
		a := coeff(L3,Dx,2)/3;
		b := normal((coeff(L3,Dx,1)-2*a^2-diff(a,x))/4);
		return [Dx^2 + a*Dx + b, 1]
	fi;

	if coeff(L6,Dx,0)=0 then
		r := 0; L := L3
	else
		# Compute a 1st order factor Dx-r of L6:
		#
		# Note: this step can be speeded up by using the
		# generalized exponents data that would have already
		# been computed anyway when trying to solve L3 with
		# DFactor or `DEtools/LiouvSols3`.
		#
		# Also note that if no ReduceOrder is possible, then
		# in most examples this would be already visible in the
		# generalized exponents.  So to make the most efficient use
		# of this code, one should first look at the generalized
		# exponents to see if it makes sense to call this code.
		#
		r := ExpSolsOverQ( diffop2de(L6,y(x)),
			y(x), `use Int`, `no algext`);

		if r=[] then
			return FAIL
		fi;
		r := normal(diff(r[1],x)/r[1]);

		# Scale L3 so that it's symmetric square has 1 as a solution
		L := symmetric_product(L3, Dx + r/2) # Sends Dx to Dx + r/2
	fi;
	L := collect(L, Dx, normal):
	if coeff(L,Dx,0)=0 then
		userinfo(1,'dsolve',"input is reducible");
		return FAIL # the reducible case is not treated.
	fi;

	# Notation: V(L) = solution space of L.
	# We'll search an operator R that acts on V(L) in such a way that
	# the induced action on V(L6) (= the symmetric square of V(L) )
	# will send 1 to 0. That means that this induced action sends
	# V(L6) to something of dimension < 6 which implies that R itself
	# sends V(L) to V(symmetric_power(L2,2)) for some L2.
	#
	# If we write R = b0 + b1*Dx + b2*Dx^2 then the induced action
	# on V(L6) depends quadratically on b0,b1,b2. So the condition
	# that 1 is sent to 0 translates into a quadratic equation (a conic)
	# C(b0,b1,b2)=0.  We only have code to solve conics over
	# fields Q(t1..tk), so algebraic coefficients can not be handled. A
	# computation lead to the following formula for conic C
	#
	a0,a1,a2 := seq(normal(coeff(L,Dx,i)),i=0..2);
	cn := diff(a0,x,x)+7/3*a2*diff(a0,x)+(4/9*a2^2+4*a1-diff(a2,x)/3)*a0;
	cd := diff(a2,x,x)+6*a0-3*diff(a1,x)+2*(diff(a2,x)-a1+2/9*a2^2)*a2;
	c := factor(cn/cd);
	userinfo(2,'ReduceOrder',"c = ",c);
	# Note: c' = 2/3*(a0-a2*c).
	C := b0^2 + c*b1^2 - 2*c*b0*b2 + diff(c,x)*b1*b2
	 + (diff(c,x,x)/2 + a2*diff(c,x)/2 + a1*c) * b2^2;
	C := collect(C,{b0,b1,b2},distributed,normal);
	userinfo(2,'ReduceOrder',"C =", C);

	# The coefficient of b0*b1 in C is zero. Now clear the coefficients
	# of b0*b2 and b1*b2.
	s0 := -c * b2;
	s1 := normal((a0/c-a2)/3) * b2;
	C := subs({b0 = b0-s0, b1 = b1-s1},C):
	C := collect(C,{b0,b1,b2},distributed,normal);
	userinfo(2,'ReduceOrder',"C after completing square =", C);
	if {coeff(C,b0),coeff(C,b1),coeff(C,b2)}<>{0} then
		error "formula incorrect"
	fi;

	# Find a rational nonzero point on C, if it exists, and read off R
	pt := Conic( coeff(C,b0,2),coeff(C,b1,2),coeff(C,b2,2) );
	userinfo(2,'ReduceOrder',"Point on C after completing square  =", pt);
	if pt = FAIL then
		return pt
	fi;

	# Multiplying by +/- 1 maps a point to another point on the conic.
	# Try 4 such points and select the one with smallest b0+b1*Dx+b2*Dx^2:
	for i in {1,-1} do for j in {1,-1} do
		b0,b1,b2 := i*pt[1], j*pt[2], pt[3];
		b0, b1 := b0-s0, b1-s1;
		S := eval(b0+b1*Dx+b2*Dx^2);
		S := collect(1/lcoeff(S,Dx) * S, Dx, normal);
		# Reverse the earlier step Dx -> Dx + r/2:
		S := symmetric_product(S, Dx - r/2);
		# Since Dx -> Dx + r/2 is an automorphism of C(x)[Dx] it
		# follows that the new S now sends L3 to a symmetric square
		if not assigned(R) or length(R) > length(S) then
			# Let R be the smallest S found:
			R := S;
			UsedPt := eval([b0, b1, b2])
		fi
	od od;
	userinfo(2,'ReduceOrder',"Used this point on C", UsedPt);
	userinfo(2,'ReduceOrder',"Gauge transformation R = ", R);

	# R sends the solutions of V(L3) to V(LR) where LR is
	LR := Dx^3 + S2*Dx^2 + S1*Dx + S0;
	LR := eval(LR, SolveLinear(
	 {coeffs(rightdivision(mult(LR, R), L3)[2],Dx)},  {S2,S1,S0}));
	userinfo(2,'ReduceOrder',"LR = ", LR);

	# Usually LR gets smaller if we multiply R by some ratpoly S
	# determined this way:
	j := convert(coeff(LR,Dx,2),parfrac,x);
	j := [`if`(type(j,`+`),op(j),j)];
	S := 1;
	for i in j do if has(denom(i),x) then
		r := normal( numer(i)/diff(denom(i),x) );
		if type(r,rational) then
			S := S * denom(i)^round(r/3)
		fi
	fi od;
	R := collect(S*R, Dx, normal);
	LR := symmetric_product(LR, Dx-normal(diff(S,x)/S));
	userinfo(2,'ReduceOrder',"Multiplied R by:", S);
	userinfo(2,'ReduceOrder',"New gauge transformation R =", R);
	userinfo(2,'ReduceOrder',"New LR =", LR);

	# L3 was assumed irreducible, so GCRD(R,L3)=1, so S*R + T*L3 = 1
	# for some S,T.  We can find S by the extended Euclidean algorithm,
	# or by equating remainder(S*R-1, L3) to 0.
	S := S0 + S1*Dx + S2*Dx^2;
	eq := rightdivision(mult(S,R)-1, L3)[2];
	eq := SolveLinear({coeffs(eq,Dx)}, {S0,S1,S2}); # equate remainder to 0
	if eq=NULL then
		userinfo(1,'dsolve',"input is reducible");
		return FAIL
	fi;
	S := eval(S, eq);
	# Now S: V(LR) --> V(L3) is the inverse of R: V(L3) --> V(LR).
	userinfo(2,'ReduceOrder',"Inverse gauge transformation S =", S);

	# The induced action of R on V(L6) had a kernel (the constants)
	# Therefore dim( V(symmetric_power(LR, 2)) ) < 6, which implies
	# that LR is of the form symmetric_power(L2, 2) for some L2
	a := coeff(LR,Dx,2)/3;
	b := normal((coeff(LR,Dx,1)-2*a^2-diff(a,x))/4);
	L2 := Dx^2 + a*Dx + b;

	# return L2 and S which maps V( symmetric_power(L2,2) ) to V(L3)
	[L2, S]
end: