# This file computes for a number of examples the conic equation one
# would get in the Singer approach and in the vHvdP approach. In most
# examples Singer's conic equation has smaller coefficients.

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

N := 7;  # to select an example, replace this by a number 1 - 7


unassign('den1');
unassign('den2');

if N=1 then
	# Example 1
	L := Dx^3-4*x*Dx-18;

elif N=2 then
	# Example 2
	L := Dx^3-4*x^3*Dx-26*x^2;

elif N=3 then
	# Example 3:
	# Only example I found so far where the second approach gave a
	# slightly smaller conic:
	L := Dx^3+1/x*Dx^2-1/9*(1+36*x^2)/x^2*Dx-4/x;
	den1 := x^2:

elif N=4 then
	# Example 4:
	L := Dx^3-(3*x^2+2)/x/(x^2+2)*Dx^2-(4*x^4+9*x^2+6)/x^2/(x^2+2)*Dx+4*x/(x^2+2);
	den2 := (x^2+2)^2 * x^2:

elif N=5 then
	# Example 5 (took the adjoint of example 4, then rescaled it)
	L := Dx^3-(10+9*x^2)/x/(x^2+2)*Dx^2-(-13*x^4-36+4*x^6-44*x^2)/(x^2+2)^2/x^2*Dx
		+4*(3*x^4+4*x^2+4)/x/(x^2+2)^2;
	den2 := (x^2+2)^4 * x^6:

elif N=6 then
	# Example 6
	L := Dx^3-9*x*Dx-27/2;

else
	# Example 7
	L := Dx^3-3/x*Dx^2-9*x*Dx-9/2;
	den2 := x^3:
fi;

# This is the 6x6 system in section 2.2 of the paper:
sys := {
diff(b00(x),x) = 2*b01(x),
diff(b01(x),x) = b11(x) + b02(x),
diff(b02(x),x) = -a0(x)*b00(x) -a1(x)*b01(x) - a2(x)*b02(x) + b12(x),
diff(b11(x),x) = 2*b12(x),
diff(b12(x),x) = -a0(x)*b01(x) -a1(x)*b11(x) - a2(x)*b12(x) + b22(x),
diff(b22(x),x) = -2*a0(x)*b02(x) -2*a1(x)*b12(x) -2*a2(x)*b22(x)
}:

# Remove denominator den1 from solution, only needed to work around
# bugs in LinearFunctionalSystems[RationalSolution].  Once those bugs
# are fixed, this line shouldn't be necessary anymore:
if assigned(den1) then
	sys := {seq(lhs(i) = rhs(i) + int(lhs(i),x) * normal(diff(den1,x)/den1), i = sys)}
fi:

sysS := subs(a2(x)=coeff(L,Dx,2),a1(x)=coeff(L,Dx,1),a0(x)=coeff(L,Dx,0), sys):

with(LinearFunctionalSystems):
RationalSolution(sysS, [b00(x),b01(x),b02(x),b11(x),b12(x),b22(x)]):
C1 := %[1]*b0^2 + 2*%[2]*b0*b1 + 2*%[3]*b0*b2 + %[4]*b1^2 + 2*%[5]*b1*b2 + %[6]*b2^2:
# This would be the conic in the Singer approach:
C1 := primpart(C1, {b0,b1,b2}):
C1 := collect(C1,{b0,b1,b2},distributed,factor);


# This is the system one would get following [v. Hoeij,v.d. Put 2006].  The
# approach is more or less dual to Singer's approach; in Singer's approach,
# one computes a first order right-hand factor of L6 which corresponds to a
# 1-dimensional quotient module, whereas in [v. Hoeij,v.d. Put 2006]
# you compute a 1-dimensional sub-module of the symmetric square module.
sys := {
diff(a00(x),x) = a0(x)*a02(x),
diff(a01(x),x) = -2*a00(x) + a1(x)*a02(x) + a0(x)*a12(x),
diff(a02(x),x) = -a01(x) + a2(x)*a02(x) + 2*a0(x)*a22(x),
diff(a11(x),x) = -a01(x) + a1(x)*a12(x),
diff(a12(x),x) = -a02(x) -2*a11(x) + a2(x)*a12(x) + 2*a1(x)*a22(x),
diff(a22(x),x) = -a12(x) + 2*a2(x)*a22(x)
}:

# Again, once the bug in LinearFunctionalSystems[RationalSolution] is fixed
# this line can be removed:
if assigned(den2) then
        sys := {seq(lhs(i) = rhs(i) + int(lhs(i),x) * normal(diff(den2,x)/den2), i = sys)}
fi:

vars := [a00(x),a01(x),a02(x),a11(x),a12(x),a22(x)]:

sysS := subs(a2(x)=coeff(L,Dx,2),a1(x)=coeff(L,Dx,1),a0(x)=coeff(L,Dx,0), sys):

RationalSolution(sysS, vars):
C2 := %[1]*a0^2 + %[2]*a0*a1 + %[3]*a0*a2 + %[4]*a1^2 + %[5]*a1*a2 + %[6]*a2^2: 
# This would be the conic in the vH.vdP approach:
C2 := primpart(C2, {a0,a1,a2}):
C2 := collect(C2,{a0,a1,a2},distributed,factor);

"length of C1 is:", length(C1);
"length of C2 is:", length(C2);