##############################
#                            #
#   Absolute Factorization   #
#                            #
##############################

# Let D = C(x)[tau] where tau is the shift operator.
# Let D_p = C(x)[tau^p] = a subring of D that is isomorphic to D.
# 
# IsomD( _ , p) is an isomorphism from D to D_p given by:  tau, x |--> tau^p, x/p.
#
# The inverse automorphism is IsomD( _ , 1/p).
#
IsomD := proc(L, p) 
	subs(_Env_LRE_tau = _Env_LRE_tau^p, _Env_LRE_x = _Env_LRE_x / p, L)
end:

# Let M = D/DL, this is a D-module and also a D_p-module.
#
# L is irreducible in D if and only if M is an irreducible D-module.
#
# Definition:  We say that L is absolutely irreducible if M is an irreducible D_p-module
# for every positive integer p.
#
# Theorem: this is equivalent to saying that M is an irreducible D_p-module for every
# prime number p that divides the order of L.
#
# Assuming that L is irreducible, our goal is to either (a) use the theorem to prove that L
# is absolutely irreducible, or (b) explicitly show that L is not absolutely irreducible by
# returning a generator of a non-trivial D_p-submodule of M.


# section_op(L, p, 1) returns the smallest operator in D_p that is divisible by L.
# This is the minimal annihilator in D_p of 1+DL in M.
# The notation in the paper is "L downarrow p".
#
# section_op(L, p)returns the image of "L downarrow p" in D under IsomD( _ , 1/p).
# The notation in the paper is "L^(p)".
#
section_op := proc(L, p)
	local a, ord, Lp, i, eqns, solutions, nf, tau;
	if nargs = 2 then return procname(L, p, p) fi;
	tau := _Env_LRE_tau;
  	ord := degree(L, tau); 
	Lp := add(a[i]*tau^(p*i), i = 0 .. ord); 
	eqns := {coeffs(LREtools[RightDivision](Lp, L)[2], tau)}; 
	solutions := SolveTools:-Linear(eqns, {seq(a[i], i = 0 .. ord)}); 
	nf := add(`if`(lhs(i) = rhs(i), 1, 0), i = solutions); 
	if nf > 1 then
		# The order is less than expected, the highest nf-1 terms are zero:
		eqns := {op(eqns), seq(a[i] = 0, i = ord - nf + 2 .. ord)};
		solutions := solve(eqns, {seq(a[i], i = 0 .. ord)})
	fi;
	Simpl_op(IsomD(eval(Lp, solutions), 1/args[3]), tau)
end:

# Simplify the operator L:
Simpl_op := proc(L, tau)
	sort(collect(primpart(L,tau),tau,factor),tau)
end:

ReducedDet := proc(v::list,x) local r,s,i;
	# Loop over v, search for same type, then combine.
	if v=[] then return v fi;
	r := NULL;
	s := v[1,2];
	for i from 2 to nops(v) do
		if IntShift(v[1,1],v[i,1],x) then
			s := s + v[i,2]
		else
			r := r, v[i]
		fi
	od;
	[`if`(s=0, NULL, [v[1,1],s]), op(procname([r],x))]
end:

# Check if two monic irreducible poly's are an integer shift from each other.
IntShift := proc(A, P, x) local n,e;
	n := degree(A,x);
	if n <> degree(P,x) then
		false
	elif lcoeff(A,x)<>1 or lcoeff(P,x)<>1 then
		procname(A/lcoeff(A,x), P/lcoeff(P,x), x)
	else
		e := Normalizer(coeff(A,x,n-1)-coeff(P,x,n-1))/n;
		type(e,integer) and Normalizer(A-subs(x=x+e,P))=0
	fi;
end:

# Input: L in C[x,tau], collected w.r.t. tau
# Output: The highest-degree terms of L that are invariant under gauge-transformations
# These are the terms whose x-degree is greater than "point on Newton polygon" - 1.
# For example, if L = a2*tau^2 + a1*tau + a0 and degree(a2,x) = 3 and degree(a1,x) = 3
# and degree(a0,x) = 4, then the slope is 1/2, and the points on the polygon are
# (0,4), (1, 3+1/2), (2,3) and then the coefficients we take are those of x^4*tau^0,
# and x^3*tau^1, and x^3*tau^2 as those are invariant under gauge-transformations.
#
InvariantTerms := proc(L,x,tau)
local n, d, degs, powd, v, L0, m, i, j, c;
	n := degree(L,tau);
	d := degree(lcoeff(L,tau),x);
	degs := [seq(degree(coeff(L,tau,i),x) - d,i=0..n)];
	powd := 0;
	v := degs[1];
	L0 := lcoeff(coeff(L,tau,0),x) * x^v;
	while powd<n do
		m := max(seq((degs[powd+j+1]-degs[powd+1])/j,j=1..n-powd));
		i := n;
		while degs[i+1]-degs[powd+1]<>(i-powd)*m do i := i-1 od;
		for j to i-powd do
			c := floor(v + j * m);
			L0 := L0 + coeff(coeff(L, tau, powd + j), x, c+d) * x^c * tau^(powd+j)
		od;
		powd := i;
		v := degs[powd+1]
	od;
	L0
end:
DetectSymSquare := proc(L,x,tau) local A;
	A := InvariantTerms(collect(L,tau),x,tau);
	testeq(Normalizer(coeff(A,tau,3) * coeff(A,tau,1)^3 - coeff(A,tau,0) * coeff(A,tau,2)^3))
end:



# Input: L in D, assumed irreducible. This assumption is not checked here (use RightFactors for that).
#
# Output:  true if L is absolutely irredubile, otherwise [p, {R_1,...}] where p is a prime number
# and where each IsomD_Dm(R_i, p) generates a non-trivial D_p-submodule(s) of D/DL, explicitly showing
# that L is not absolutely irreducible.
#
AbsFactorization := proc(L)
	local i,ord,p,c,G,L_tilde,L_p,S,x,tau,N,t;
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	N := InvariantTerms(L, x, tau);
	coeffs(N, tau, 't');
	# Only primes p for which N in Q(x)[tau^p] can occur:
	for p in map(i -> i[1], ifactors(igcd(op(map(degree,[t],tau))))[2]) do
	L_p := section_op(L, p); 
	if degree(L_p, tau) < degree(L, tau) then 
		return [p, {1}]
	elif p=2 and args[-1] = `Hom method` then
		# Step 1(a) in Section 3.3:
		L_tilde := subs(tau = -tau, L);
		# Step 1(b) in Section 3.3: (this step needs the file "Hom.txt")
		G := Hom(L, L_tilde);
		if G <> [] then
			if nops(G) > 1 then error "Input was not irreducible" fi;
			G := G[1];
			c := LREtools[MultiplyOperators](subs(tau=-tau,G), G);
			c := LREtools[RightDivision](c, L)[2];
			if has(c, {x,tau}) then error "should not happen" fi;
			G := collect(G/sqrt(c), tau, normal);
			return [p, {seq(Simpl_op(IsomD(LREtools[RightDivision](LREtools[LCLM](L, 1 + i*G), 1 + i*G)[1], 1/2), tau), i = {1,-1})}]
		fi
	else
		# Pass along data to RightFactors to only compute factors that are consistent with this data.
		# Need to load an update for RightFactors for this to make a difference.
		#
		_Env_RightFactors_data := [p, IsomD(N, 1/p), ReducedDet(select(has,factors(tcoeff(L,tau)/lcoeff(L,tau))[2],x), x)];

		S := LREtools[RightFactors](L_p, degree(L,tau)/p);
		if S <> {} then
			return [p,S]
		fi
	fi
	od;
	# Step 2
	true
end:


# LREtools[RightFactors] makes a list of "candidate determinants" for
# the factors that it is trying to construct. The larger this list is,
# the longer the computation will take.
# However, some of these candidate determinants can be discarded if
# we know that the input of RightFactors comes from AbsFactorization.
# Given a "candidate determinant" v, the following program tells a
# modification of RightFactors if we should keep v or not.
# 
# To measure if this improvement makes a difference, you can turn off
# this improvement by inserting "return true" at the beginning.
DetFactorsSelect := proc(v, x, tau)
	local N,f,i;
	N := _Env_RightFactors_data[2];
	f := mul(i[2],i=v[-1]);
	if N=`3-2` then
		evalb(degree(f,tau)=3 and DetectSymSquare(f,x,tau))
	else
		testeq(rem(InvariantTerms(collect(f,tau),x,tau),N,tau))
	fi
end:


# LREtools[RightFactors] also makes a list of "candidate Invariant terms". But like
# the previous program, some of these may be discarded if we know that the input comes
# from AbsFactorization.
#
DeterminantSelect := proc(v, x)
	local N;
	N := _Env_RightFactors_data;
	evalb([] = ReducedDet( [op(N[3]), op(select(has,factors(subs(x=x/N[1], 1/v))[2],x))], x))
end:




##############################
#                            #
#   Reduce Order 3 --> 2     #
#                            #
##############################


# Input: Order 3 element L in Q[x][tau] whose SymSquare has order 5.
# Output: b,r for which L = SymSquare(tau^2+tau+b) (s) tau-r.
#
Theorem51 := proc(L, x, tau) local i, c, dA, b, r, g;
	if degree(L,tau) <> 3 then
		error "expecting an operator of order 3"
	fi;
	for i from 0 to 2 do
		c[i] := normal( coeff(L,tau,i)/coeff(L,tau,3) )
	od:
	if c[1] = 0 and c[2] = 0 then
		return "Case 3(a)"
	elif c[1] = 0 or c[2] = 0 then
		return false  # (I in Theorem 5.1(1) is not zero)
	elif normal(c[0] - subs(x=x-1, c[1]) * c[2]) = 0 then
		return "Case 3(b)"
	fi;
	dA := 1 - subs(x=x-1,c[1]) * c[2]/c[0];
	b := normal( (1 - subs(x=x-1,c[2]) * c[2]/c[1])/dA );
	b := normal(b);
	r := normal( subs(x=x-2, c[2])/ (subs(x=x-1, b)-1) );
	g := normal( ( 1-c[1]/c[0] * subs(x=x-1, c[1]/c[2]) )/dA );
	if normal(subs(x=x+1, b) - g) <> 0 then
		return false
	fi;
	[b, r]
end:


# The extended Euclidean Algorithm for D = C(x)[tau]
#
GCDEX := proc(a, b) 
	local r1, q1, s, t; 
	if b = 0 then 
		1/content(a, _Env_LRE_tau), 0
	else
		q1, r1 := LREtools[RightDivision](a, b);
		t, s := procname(b, r1);
		s, collect(t-s*q1, _Env_LRE_tau, normal)
	fi
end:

# Compute the m'th SYMMetric power of L1:
symmpower := proc(L1, m)  global x,tau ;
	local c, n1, sb1, i, P, var, j, IsZero, EQN, s, NumberFree, k, N;
	n1 := degree(L1, tau); 
	sb1 := u(x + n1) = -add(coeff(L1, tau, i)*u(x + i)/lcoeff(L1, tau), i = 0 .. n1 - 1); 
	P[0] := u(x)^m; 
	N := binomial(m + n1 - 1, n1 - 1); 
	for i to N do 
		P[i] := subs(sb1, subs(sb1, subs(x = x + 1, P[i - 1])))
	od;
	var := {seq(u(x + i), i = 0 .. n1 - 1), seq(u(x + j), j = 0 .. n1 - 1)}; 
	IsZero := collect(add(c[i]*P[i], i = 0 .. N), var, distributed); 
	EQN := {coeffs(IsZero, var)}; 
	var := {seq(c[i], i = 0 .. N)}; 
	s := solve(EQN, var); 
	NumberFree := add(`if`(lhs(i) = rhs(i), 1, 0), i = s); 
	if 1 < NumberFree then
		s := solve(s union {seq(c[k] = 0, k = N + 2 - NumberFree .. N)}, var)
	fi;
	sort(collect(primpart(subs(s, add(c[i]*tau^i, i = 0 .. N)), tau), tau, factor), tau)
end:

# Input: Order 3 element of Q(x)[tau].
# Output: If L gauge-equivalent to some SymSquare(tau^2+tau+b) (s) tau-r
# then the transformation, its inverse, tau^2+tau+b and r.
# Otherwise:  FAIL
# With the optional argument `projective` we just return SimpIB(tau^2+tau+b, `projective`)
# which tries to find a smaller operator L2 that is projectively equivalent to it.
#
ReduceOrder32 := proc(L3)
	local x,tau, L6, L1, ux, i, G2, G, S, vars, rel, rel1, rel2, C, pt, LG, rel3, X, a, b, d, l;
	tau := _Env_LRE_tau;
	x := _Env_LRE_x;

    # Run some checks to see if L3 could be a symmetric square:
	if not DetectSymSquare(L3,x,tau) then return FAIL fi;
	d := tcoeff(L3,tau)/lcoeff(L3,tau);
	l := surd(lcoeff(d,x),3) * x^(degree(d,x)/3); # Can also replace tau-l^2 with 1.
	d := ReducedDet(select(has,factors(d)[2],x),x);
	if convert(map(i -> i[2] mod 3, d),`+`) <> 0 then return FAIL fi;
	_Env_RightFactors_data := [1, tau-l^2, [seq([i[1], 2*i[2]/3], i = d)] ];

    # Step 1:
	L6 := symmpower(L3,2);

    # Step 2:
	if degree(L6,tau) < 6 then
		i := Theorem51(L3, x, tau);
		if not type(i,list) then return i fi;
		if args[-1]=`projective` then return SimpIB(tau^2+tau+i[1], `projective`) fi;
		return 1, 1, tau^2+tau+i[1], i[2]
	fi;

    # Step 3:
	L1 := LREtools[RightFactors](L6, 1);
	if nops(L1)<>1 then return FAIL fi;
	L1 := L1[1];

    # Step 4:
	ux[3] := solve(LREtools[OperatorToRecurrence](L3, u(x)), u(x + 3));
	ux[0] := u(x);
	for i to 5 do ux[i] := normal(subs(x = x + 1, u(x + 3) = ux[3], ux[i - 1])) od;

    # Step 5:
	G2 := add(a[i]*ux[i]^2, i = 0 .. 5);
	G := add(b[i]*ux[i], i = 0 .. 2);

    # Step 6:
	S := {coeffs(collect(G2 - G^2, {seq(ux[i], i = 0 .. 2)}, distributed), {seq(ux[i], i = 0 .. 2)})};
	G2 := add(a[i]*tau^i, i = 0 .. 5);
	vars := indets(G2) minus {tau};

    # Steps 7 and 8:
	rel := collect(LREtools[RightDivision](subs(solve(S, vars), G2), L1)[2], {b[0], b[1], b[2]});
	rel1 := b[0] = X[1] - 1/2*coeff(rel, b[0], 1);
	rel := collect(subs(rel1, rel), {X[1], b[1], b[2]}, Normalizer);
	if degree(rel,b[1]) = 2 then
		rel2 := b[1] = X[2] - 1/2*coeff(rel, b[1], 1)/coeff(rel, b[1], 2);
		rel3 := b[2] = X[3];
    # Step 9:
		C := collect(primpart(subs(rel2,rel3, rel), {X[1], X[2], X[3]}), {X[1], X[2], X[3]});
    # Step 10, code from www.math.fsu.edu/~hoeij/files/ConicProgram
		pt := Conic(expand(coeff(C, X[1], 2)), expand(coeff(C, X[2], 2)), expand(coeff(C, X[3], 2))); 
		if pt = FAIL then return FAIL fi;
		S := solve({pt[1] - X[1], pt[2] - X[2], pt[3] - X[3], rel1, rel2, rel3}, {X[1], X[2], X[3], b[0], b[1], b[2]})
	else
		S := rel1, X[1] = 0, b[1]=1, b[2]=0
	fi;

    # Step 11:
	G := subs(S, add(b[i]*tau^i, i = 0 .. 2)); 

    # Step 12:
	LG := LREtools[RightDivision](LREtools[LCLM](L3, G), G)[1];
	LG := collect(LG/lcoeff(LG, tau), tau, factor);
	i := Theorem51(LG, x, tau);
	if type(i, list) then
		if args[-1]=`projective` then return SimpIB(tau^2+tau+i[1], `projective`) fi;
		collect(G, tau, factor), collect(GCDEX(G,L3)[1], tau, factor), tau^2+tau+i[1], i[2]
	else
		i
	fi 
end:


###################################
#                                 #
#   Reduce Orders {3,4} --> 2     #
#                                 #
###################################


# Input: an operator L4 of order 3 or 4
# This program checks if L4 is 2-solvable, and if so, reduces L4 to order 2.
#
ReduceOrder := proc(L4)
	local A,i,x,tau,c0,c1,c2,c3,c4,d,b1,b2;
	x, tau := _Env_LRE_x, _Env_LRE_tau;

	lprint("Factorization...");	# Add optional flag to skip this?
	A := [seq(LREtools[RightFactors](L4, i), i=1..3)];
	if degree(L4,tau) = 3 then
		A := A[1..2];
		if A<>[{},{}] then
			return ["Reducible, right factors are", op(A)]
		else
			lprint("Calling ReduceOrder32...");
			return ReduceOrder32(args)
		fi
	elif A <> [{},{},{}] then
		if nops(A[3])=1 and A[2]={} and nops(A[1]) <= 1 then
			i := ReduceOrder32(A[3][1],`projective`)
		fi;
		return ["Reducible, right factors are", op(A), `if`(has(i,tau),op(["which is solvable in terms of", i]),NULL)]
	fi;

	A := InvariantTerms(L4,x,tau);
	c0,c1,c2,c3,c4 := seq(coeff(A,tau,i),i=0..4);
	if c1=0 and c3=0 then
		if coeff(L4,tau,1)=0 and coeff(L4,tau,3)=0 then
			return ["Input is the image of", IsomD(L4,1/2), "under isomorphism C(x)[tau] -> C(x)[tau^2],  tau,x |-> tau^2, x/2."]
		fi;
		lprint("Absolute factorization...");
		A := AbsFactorization(L4);
		if type(A,list) then
			return ["Absolute factorization, right factors of operator for u(2*n) are", A[2]]
		fi
	fi;
	d := ReducedDet(select(has,factors(tcoeff(L4,tau)/lcoeff(L4,tau))[2],x),x);
	if member(1, map(i -> i[2] mod 2, d)) then return FAIL fi;
	_Env_RightFactors_data := [1, `3-2`, [seq([i[1], 3*i[2]/2], i = d)] ];
	b1,b2 := TestSum(c4*c1^2, -c0*c3^2), TestSum(c4*c1^2, c0*c3^2);
	if b1 or b2 then
		lprint("Checking Symmetric Product...");
		A := CaseSymProd(L4,b1,b2);
		if not has(A, FAIL) then
			return [cat(`if`(A[-1]=2,"Operator for u(2*n)","Input"), " is projectively equivalent to SymProd of"), A[1]]
		fi
	fi;
	if b1 and TestSum(c3*c1^2*c4, - 2*c1*c2^2*c4, c3^2*c1*c2, c2^3*c3, - c1*c3^4/c4) then
		lprint("Checking Symmetric Cube... (can be time consuming...)");
		A := CaseSymCube(L4);
		if A <> FAIL then
			return ["input is projectively equivalent to the symmetric cube of", A]
		fi
	fi;
	FAIL
end:
TestSum := proc()
	local x,d;
	x := _Env_LRE_x;
	d := max(map(degree,[args],x));
	testeq(coeff(convert([args],`+`),x,d))
end:

# L4 is an irreducible order 4 operator
# If L4 is projectively equivalent to some SymProd(L2a,L2b) then the output is [{L2a,L2b},1]
# If section_op(L4,2) is projectively equivalent to some SymProd(L2a,L2b) then [{L2a,L2b},2]
# Otherwise FAIL.

CaseSymProd := proc(L4,b1,b2)
	local L6,R,b,i;
	L6 := ExtSquare(L4);
	R := `if`(b1, LREtools[RightFactors](L6,3), {});
	b := 1;
	if R={} then
		lprint("Checking operator for u(2*n)...");
		R := IsomD(L6, 1/2);
		if not (b2 and DetectSymSquare(add(coeff(L6,tau,2*i)*tau^i,i=0..3),x,tau)) then
			return FAIL
		elif type(R,polynom(anything,tau)) then
			R := {R, subs(x=x+1/2,R)}
		else
			R := AbsFactorization(L6);
			if R=true then return FAIL fi;
			if R[1]=3 then error "unexpected" fi;
			R := R[2]
		fi;
		b := 2;
	fi;
	[map(ReduceOrder32, R, `projective`), b]
end:
# Returns the exterior square of an operator of order 4.
ExtSquare := proc(L)
	local u1, u2,sb,i,P,IsZ,S, x,tau;
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	sb := u1(x+4) = -add( coeff(L,tau,i)/lcoeff(L, tau) * u1(x+i), i=0..3),
	u2(x+4) = -add( coeff(L,tau,i)/lcoeff(L, tau) * u2(x+i), i=0..3);
	P[0] := u1(x)*u2(x+1) - u1(x+1)*u2(x);
	for i to 6 do
		P[i] := subs(x=x+1, sb, P[i-1]);
		P[i] := collect(P[i], indets(P[i],function), normal)
	od;
	IsZ := add(a[i] * P[i], i=0..6);
	IsZ := collect(IsZ, indets(IsZ,function),distributed);
	S := solve({coeffs(IsZ, indets(IsZ,function))}, {seq(a[i], i=0..6)});
	Simpl_op( subs(S, add(a[i]*tau^i, i=0..6)), tau)
end:


# Input: L4 is an irreducible order 4 operator
# Output = FAIL if L4 is not proj-equivalent to some L2 symm^3.
# Otherwise output = [L2, transformation]

CaseSymCube:=proc(L4,d)
	local L10,L3;
	L10:=symmpower(L4,2);
	if degree(L10,tau)=7 then
		return SpecialCase(L4)
	fi;
	_Env_RightFactors_SymCubeCase := true;
	L3 := LREtools[RightFactors](L10,3);
	if L3={} then return FAIL fi;
	ReduceOrder32(L3[1], `projective`);
end;:


# Input: L is an irreducible order 7 operator
# Output: L2 or FAIL
SpecialCase:=proc(L)
	local P,i,s,IsZ,var,S,L3,L2,sb1;
	P[0] := u(x)^2;
	sb1 := u(x+4) = -add(coeff(L, tau, i)*u(x+i)/lcoeff(L, tau), i = 0..3);
	for i from 1 to 6 do
		P[i] := subs(sb1, subs(sb1, subs(x = x + 1, P[i-1])))
	od;
	var := {seq(a[i], i=0..10)};
	for s from 1 to 3 do
		P[7] := u(x)*u(x+s);
		for i from 8 to 10 do
			P[i] := subs(sb1, subs(sb1, subs(x = x + 1, P[i-1])))
		od;
		IsZ := add(a[i] * P[i], i=0..10);
		IsZ := collect(IsZ, indets(IsZ,function), distributed);
		S := solve({coeffs(IsZ, indets(IsZ,function))}, var);
		if nops(indets(map(rhs,S)) intersect var ) > 1 then next fi;
		L3 := primpart(subs(S, add(a[i]*tau^(i-7), i=7..10)),tau);
		return ReduceOrder32(L3,`projective`)
	od;
	lprint("error: u(x)u(x+s) should have worked for some s in {1,2,3}");
	FAIL
end:



########################################
#                                      #
#   Simplify an operator of order 2    #
#                                      #
########################################


# Input: an irreducible polynomial f in Q[x].
# Output: a polynomial, shift-equivalent to f, computed in such a way
# that if two inputs f1,f2 are shift-equivalent, then the output is the same.
#
# For example, if f has a rational root, then the output is x-r where r in [0,1)
# differs from the rational root by an integer.
#
NormalizePoly := proc(f, x) local g, d, p, e; global NormalizePolyTable;
	g := collect(f/lcoeff(f,x), x, Normalizer);
	d := degree(f,x);
	if type(g, polynom(rational,x)) and degree(g,x) < 3 then
		collect(subs(x = x - ceil(coeff(g,x,d-1)/d), g), x, Normalizer)
	else
		# If f is not in Q[x] then we need to do some more work:
		#
		if not assigned(NormalizePolyTable) then NormalizePolyTable := {} fi:
		for p in NormalizePolyTable do if degree(p,x) = d then
			e := Normalizer(coeff(g-p,x,d-1)/d);
			if type(e,integer) then
				return p
			fi
		fi od;
		NormalizePolyTable := NormalizePolyTable union {g};
		g
	fi
end:

# If b is in D/DL, where D = Q(x)[tau, tau^(-1)] and L is in D, and
# if b(left-sols) and b(right-sols) have valuations at x = i+R that are
# at at least as high as the output of the program below, then b
# is considered integral at x = i + R.
#
# Here R is the RootOf of a normalized irreducible polynomial. Common cases:
# if the original polynomial had an integer root, then R = 0,
# if it has a rational root, then R is in [0,1).

# Input: i, gmin, gmax
# Output: [allowed valuation for left-sols, allowed valuation for right-sols]
#
AllowedValuation := proc(i, gmin, gmax)
	if i > 0 then
		[gmin, 0]
	else
		[0, -gmax]
	fi
end:

# To figure out if b = 1 + DL in D/DL is integral at x = i+R, we need the actual
# valuations of left-sols and right-sols at x = i+R. These are computed here.
#
# Notation:
#	L is in Q(x)[tau] where tau is the shift operator
#	R = RootOf(a normalized irreducible polynomial)
#	gp = `LREtools/g_p`(L, R)
#	n = order of L
#	gmin = min(gp[1])
#	gmax = max(gp[1])
#
# Input: integer i, and gp, n, gmin, gmax as above.
# Output: valuation(left-sols) and valuation(right-sols) at x = i+R.
#
ActualValuation := proc(i, gp, n, gmin, gmax) local a0,an,left,right,c,d;
	a0,an,left,right := op(gp[-4..-1]);
	if n=1 then an := [seq(c-an[-1], c=an)] fi; # (fixes a quirk in the implementation)

	# left-sols  have valuation 0 at left-n  .. left-1  (the first n 0's in a0)
	# right-sols have valuation 0 at right+1 .. right+n (the last  n 0's in an)

	c := n+1+i-left;
	d := i-right-n+nops(an);

	[ `if`( c <= 0,     0, `if`(c > nops(a0), gmin, a0[c])) ,
	  `if`( d <= 0, -gmax, `if`(d > nops(an),    0, an[d])) ]
end:

# To make tau^k integral in D/DL, we need to multiply it with the output of this program.
# Input:
#	S = set of irreducible polynomials
#	gp = a table, such that if p in S, then gp[p] = `LREtools/g_p`(L, RootOf(p))
#
MakeIntegral := proc(k, n, S, gp) local P, p, g, i, gmin, gmax;
	P := 1;
	for p in S do
		g := [gp[p]]; gmin := min(g[1]); gmax := max(g[1]);
		P := P * mul(collect(subs(x=x-i, p),x)^max(
			AllowedValuation(i, gmin, gmax) - ActualValuation(k+i, g, n, gmin, gmax)),
		i = min(1,g[-2]-k) .. max(0,g[-1]-k))
	od;        #    ^ cut-off-points  ^  These need to match the cut-off in AllowedValuation.
	P
end:

# The set of singularities of L is defined as:
#
#	{ NormalizePoly(f) | f is an irreducible factor of lcoeff(L,tau) or tcoeff(L,tau) }
#
# To improve the efficiency, we use a gcd-trick to delete some (but not necessarily all) of
# the so-called apparent singularities, defined as singularities where gmin = gmax = 0 and
# moreover ActualValuation(i) >= [0, 0] for all i.
#
Singularities := proc(L, x, tau) local i, d, M, gp, P;
	options remember;
	if not type(L, polynom(anything,x)) then
		return procname(primpart(L,tau), x,tau)
	fi;
	d := degree(L,tau);
	i := iquo(d, 2);
	M := coeff(L, tau, i+1) * subs(x=x+1, L) * tau - subs(x=x+1, coeff(L, tau, i)) * L;
	M := primpart(M, tau);
	M := gcd(subs(x=x+1, lcoeff(L,tau)), coeff(M,tau,d+1)) * gcd(tcoeff(L,tau), coeff(M,tau,0));
	P := {seq(NormalizePoly(i[1],x), i=sort(factors(M)[2],length))};
	for i in P do
		gp[i] := `LREtools/g_p`(L, RootOf(i,x))
	od;
	P, gp
end:

# Compute an Integral Basis of L in triangular form.
IB := proc(L) local x,tau, i,S,n,B,j,unchanged,tj;
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	i := ldegree(L,tau);
	if not has(L,tau) then
		error "Input should have order > 0"
	elif i <> 0 then
		return procname( LREtools[MultiplyOperators](tau^(-i), L) )
	fi;
	S := Singularities(L, x, tau);
	n := degree(L,tau);
	B := [seq(MakeIntegral(i, n, S) * tau^i, i=0..n-1)];
	for i do
		unchanged := true;
		for j in [n-1+i, -i] do
			tj := MakeIntegral(j, n, S) * tau_i(j, n,L,x,tau);
			B := AddToTriangularBasis('unchanged', tj , B, x, tau);
		od;
		if unchanged then return B fi
	od
end:
# Compute tau^i mod L.  Could merge this into program IB.
tau_i := proc(i, n,L,x,tau) local a;
	options remember;
	if i >= 0 and i < n then
		return tau^i
	fi;
	if i < 0 then
		a := tau_i(i+1, n,L,x,tau);
		a := a - coeff(a,tau,0)/tcoeff(L,tau) * L;
		collect(LREtools[MultiplyOperators](tau^(-1), a), tau, normal)
	else
		a := LREtools[MultiplyOperators](tau, tau_i(i-1, n,L,x,tau));
		collect(a - coeff(a,tau,n)/lcoeff(L,tau) * L, tau, normal)
	fi
end:
AddToTriangularBasis := proc(unchanged, v, B, x, tau) local f,d;
	f := ModTriangularBasis(v, B, x, tau);
	if f = 0 then return B fi;
	unchanged := false;
	d := degree(f,tau) + 1;
	procname('unchanged', B[d], subsop(d = f, B), x, tau)
end:
ModTriangularBasis := proc(f, B, x, tau)
	local v,d,q;
	v := f;
	for d from degree(f,tau) to 0 by -1 do
		q := normal( coeff(v,tau,d)/coeff(B[d+1], tau, d) );
		v := collect(v - quo(numer(q), denom(q), x) * B[d+1], tau, Normalizer)
	od;
	`if`(v=0, 0, v / icontent(lcoeff(v,tau)))
end:



# B is in D/DL.  Compute  B wedge tau(B).
DetB := proc(B, x,tau,L) local d;
	d := LREtools[RightDivision](LREtools[MultiplyOperators](tau,B),L)[2];
	d := Normalizer( coeff(B,tau,0)*coeff(d,tau,1) - coeff(B,tau,1)*coeff(d,tau,0) );
	numer(d), degree(denom(d),x)
end:

# The outputs "SoFar" should minimize degree(DetB( b1 - SoFar*b2 )).
#
Reduceb1modb2 := proc(b1,b2,  x,tau,L, d, SoFar) local c,dt,dn,l,v,i,dl;
	dt,dn := DetB(collect(b1 - (SoFar + c*x^d) * b2, tau), x,tau,L);
	l := lcoeff(dt, x);   # Use normal-lcoeff
	v := [seq(`if`(degree(i[1],c)=1,-coeff(i[1],c,0)/coeff(i[1],c,1),NULL), i = factors(l)[2])];
	if v=[] then
		dl := degree(l,c); v := [0]
	else
		dl := 1
	fi;
	if d=0 then		# Use degree-normal-lcoeff
		seq([SoFar + l, [seq(degree(i,x)-dn, i=[coeff(dt,c,2), dt, eval(dt,c=l)])], dl], l=v)
	else
		seq(procname(b1,b2, x,tau,L, `if`(dl=1,d-1,0), SoFar + l*x^d), l=v)
	fi
end:

# Among the outputs of Reduceb1modb2 find the best one.
#
MinReduceb1modb2 := proc(b1,b2,  x,tau,L, ds) local d,M,i,m;
	d := `if`(nargs>5, ds[2]-ds[3], 0);
	M := [Reduceb1modb2(b1,b2,  x,tau,L, d, 0)];
	if nargs>5 and ds[2]>ds[1] and ds[2]>ds[3] then
		# SoFar = 0 must be among the outputs, delete it so that b1 doesn't reduce to b1.
		M := subsop(seq(`if`(M[i,1]=0,i,NULL),i=1..nops(M)) = NULL, M);
		if M = [] then
			return [FAIL, [infinity$3]]
		fi
	fi;
	m := min(map(i -> i[2,3], M));
	M := [seq(`if`(i[2,3]=m, i, NULL), i = M)];	# Only those with minimal deg(b1-SoFar * b2)
	m := min(map(i -> i[2,2], M));
	M := [seq(`if`(i[2,2]=m, i, NULL), i = M)];	# Only those with minimal total degree (of b2, b1-SoFar*b2).
	M := sort(M,length)[1];
	subsop(1 = collect(b1-M[1]*b2, tau, Normalizer), M)
end:


# Input: an integral basis b1,b2.
# Output: a reduced basis, by repeatedly using MinReduceb1modb2, as well as
# a set of the best elements encountered.
#
ReduceBasis := proc(b1,b2, x,tau,L, d1,d,d2) local b,ds,best,D,s,i,m;
	if nargs = 5 then
		b := MinReduceb1modb2(b1,b2,  x,tau,L);
		procname(b2, b[1], x,tau,L, op(b[2]))
	elif d1=-infinity or d2=-infinity then
		[b1,b2], {`if`(d1=-infinity,b1,NULL), `if`(d2=-infinity,b2,NULL)}
	elif d > d1 and d > d2 then
		# Compute a Zsequence
		ds := [d1,d,d2];
		best := convert(ds,`+`), 0;  # best basis so far is b[0],b[1]
		b[0],b[1] := b1,b2;
		D[0],D[1] := d1,d2;
		for i to 4 while b[i]<>FAIL do
			b[i+1], ds := op(1..2, MinReduceb1modb2(b[i-1],b[i], x,tau,L, ds));
			s := convert(ds,`+`);
			if s < best[1] then
				if s = -infinity then return [b[i],b[i+1]], {b[i+1]} fi;
				best := s, i
			fi;
			D[i+1] := ds[3];
		od;
		ds := [d2,d,d1];
		for i to 4 while b[1-i]<>FAIL do
			b[-i], ds := op(1..2, MinReduceb1modb2(b[2-i],b[1-i], x,tau,L, ds));
			s := convert(ds,`+`);
			if s < best[1] then
				if s = -infinity then return [b[1-i],b[i]], {b[-i]} fi;
				best := s, -i
			fi;
			D[-i] := ds[3]
		od;
		m := min(seq(`if`(b[i]=FAIL or not assigned(D[i]),NULL,D[i]),i = -4..5));
		i := best[2];
		[b[i],b[i+1]], {seq(`if`(D[i]=m, b[i], NULL), i=-5..6)}
	else
		b := MinReduceb1modb2(b1,b2, x,tau,L, [d1,d,d2]);
		ds := b[2];
		if b[-1] = 2 then
			[b[1], b2], {}
		elif b[-1] = 0 then
			[b[1], b2], {`if`(ds[3] < ds[1], b[1], b2)}
		else
			procname(b2, b[1], x,tau,L, op(ds));
		fi
	fi
end:
		
		
# Given an element R in D/DL,  find an operator M with MR is 0 in D/DL,
# or return a right-factor of L if R generates a 1-dimensional submodule.
#
MinOp2 := proc(R, x,tau,L) local A, a2,a1,a0;
	A := a2*tau^2 + a1*tau + a0;
	A := subs(solve({coeffs(LREtools[RightDivision](LREtools[MultiplyOperators](A, R), L)[2],tau)},{a2,a1,a0}), A);
	A := primpart(A, tau);
	if has(A,{a2,a1,a0}) then
		lprint("Reducible case, returning right-factor instead:");
		A := primpart(R, tau)
	fi;
	sort(collect(A * sign(lcoeff(A,tau)), tau, factor), tau)
end:

# Input: an operator of order 2.
# Output: if order 1: a right-factor of L.
#         if order 2: a simplified gauge-equivalent operator.
#
SimpIB := proc(L) local x,tau, i, B, c1, c2, R;
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	i := ldegree(L,tau);
	if i <> 0 then
		return procname( LREtools[MultiplyOperators](tau^(-i), L), args[2..-1] )
	elif degree(L,tau)<>2 then
		error "This program only simplifies operators of order 2"
	elif args[-1] = `projective` then
		R := procname(UpdateDet(L,x,tau));
		i := gcd(subs(x=x-1,lcoeff(R,tau)), coeff(R,tau,0));
		while has(i,x) and degree(R,tau)=2 do
			R := normal(lcoeff(R,tau)/subs(x=x+1,i))*tau^2 + coeff(R,tau,1)*tau + coeff(R,tau,0)*i;
			i := procname(R);
			if length(i) < length(R) then
				R := i; i := gcd(subs(x=x-1,lcoeff(R,tau)), coeff(R,tau,0));
			else
				break
			fi
		od;
		return R
	fi;
	B := ReduceBasis(op(IB(L)), x,tau,L);
	# B[1] = reduced basis, and B[2] = {best elements found so far}.
	if B[2]={} then
		B := MinOp2(c1*B[1,1]+c2*B[1,2], x,tau,L);
		if not has(B,{c1,c2}) then
			return B
		fi;
		R := lcoeff(coeff(B,tau,1),x) * resultant(coeff(B,tau,2),coeff(B,tau,0),x);
		R := factors(R)[2];
		R := [seq(`if`(degree(i[1],{c1,c2})=1, solve(i[1],{c1,c2}), NULL), i=R)];
		sort([seq(sort(collect(primpart(eval(B,i),tau),tau,factor),tau), i=R)], length)[1]
	else
		sort([seq(MinOp2(i, x,tau,L), i = B[2])], length)[1]
	fi;
end:
UpdateDet := proc(L, x, tau)
	local a0,a1,a2,r,V,S,i,j;
	a0,a1,a2 := seq(coeff(L,tau,i),i=0..2);
	r := factors(a0/a2);
	V := [seq([NormalizePoly(i[1],x),i[2]], i = r[2])];
	S := {seq(i[1],i=V)};
	S := {seq([primpart(i,x), add(`if`(j[1]=i, j[2], 0), j=V)], i = S)};
	r := ifactors(r[1]/mul(lcoeff(i[1],x)^i[2], i = S));
	r := mul(i[1]^`if`(abs(i[2])>1, floor(i[2]/2), 0), i = [op(r[2]),op(S)]);
	r := r * subs(x=x+1,r) * a2 * tau^2 + r * a1 * tau + a0;
	sort(collect(primpart(r,tau),tau,factor),tau)
end:


################################
#                              #
#   Tools to make examples     #
#                              #
################################


# Calculates the Symmetric product.
# This is needed when (remains to do) we want to compute the transformation.
SymProd := proc(L1, L2)
	local u1, u2,n1,n2,sb,i,P,N,IsZ,f,S,k,x,tau;
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	n1, n2 := degree(L1,tau), degree(L2,tau);
	sb := u1(x+n1) = -add( coeff(L1,tau,i)/lcoeff(L1, tau) * u1(x+i), i=0..n1-1),
	u2(x+n2) = -add( coeff(L2,tau,i)/lcoeff(L2, tau) * u2(x+i), i=0..n2-1);
	P[0] := u1(x)*u2(x);
	N := n1*n2;
	for i to N do
		P[i] := subs(x=x+1, sb, P[i-1]);
		P[i] := collect(P[i], indets(P[i],function), normal)
	od;
	IsZ := add(a[i] * P[i], i=0..N);
	f := indets(IsZ,function) minus indets([seq(subs(x=x+i,[args]),i=0..N)],function);
	IsZ := collect(IsZ, f,distributed);
	S := solve({coeffs(IsZ, f)}, {seq(a[i], i=0..N)});
	k := add(`if`(lhs(i)=rhs(i),1,0), i=S) - 1;
	if k>0 then
		S := solve(S union {seq(a[N-i],i=0..k-1)}, {seq(a[i], i=0..N)})
	fi;
	Simpl_op(subs(S, add(a[i]*tau^i, i=0..N)), tau)
end:

# For this to work you need to load the file "Hom.txt"
ProjectiveHom := proc(L1, L2)
	local x,tau,n,d,R,i;
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	n := degree(L1,tau);
	if degree(L2,tau)<>n then
		error "Expecting inputs of the same order"
	fi;
	d := tcoeff(L2,tau)/tcoeff(L1,tau) * lcoeff(L1,tau)/lcoeff(L2,tau);
	d := factors(d);
	d := surd(d[1],n) * mul(i[1]^(i[2]/n), i=ReducedDet(d[2],x));
	if not type(d,ratpoly(anything,x)) then
		return FAIL
	fi;
	for i in {d, `if`(irem(n,2)=0, -d, NULL)} do
		R := SymProd(L1, tau-i);
		if testeq(rem(InvariantTerms(R,x,tau),InvariantTerms(L2,x,tau),tau)) then
			R := Hom(R, L2);
			if R <> [] then
				return ["SymProd with",tau-i,") and then gauge transformation", op(R)]
			fi
		fi
	od;
	FAIL
end:


##################################################
#                                                #
#   Find a point on a conic over Q(t1..tk)       #
#   www.math.fsu.edu/~hoeij/files/ConicProgram   #
#                                                #
##################################################


# Input: a,b,c in Q(t1,t2,..,tk).
#
# Output: nonzero [x,y,z] for which a*x^2+b*y^2+c*z^2 = 0 if
#         such x,y,z exist in Q(t1,t2,..,tk) and FAIL otherwise.
#
# Author: Mark van Hoeij, Oct 2004.
#
Conic := proc(a,b,c)
	local v,d,i,t,l,s,d2,co,ansatz,c1,j,k,eq;

	# Some preprocessing:
	#
	v := [a,b,c];
	if   a=0 then [1,0,0]
	elif b=0 then [0,1,0]
	elif c=0 then [0,0,1]
	elif not type(v, list(ratpoly(rational))) then
		FAIL # Only Q(t1..tk) is implemented.
	elif args[-1]<>0 then
		d, v := ReduceConic(ReduceList(v));
		v := procname(op(v),0);
		if v = FAIL then return v fi;
		ReduceList([seq(v[i]/d[i], i=1..3)])
	elif indets(v)={} then
		# Solve a conic over Q (requires factoring integers)
		`algcurves/pyth_LLL`(op(v))
	else

	t := SelectVar(v);
	d := map(degree,v,t);
	l := map(lcoeff,v,t);
	s := d[1]+d[2]+d[3];
	d := [seq(s-d[i], i=1..3)];
	d2 := d mod 2;
	# The degrees for the ansatz:
	d := [seq(ceil(i/2)-max(op(d2)), i=d)];
	ansatz := [seq(add(co[i,j]*t^j, j=0..d[i]), i=1..3)];
	eq := {};

	# If the degrees are congruent mod 2, then we need some equations
	# coming from a conic with fewer variables:
	if d2 = [0,0,0] then
		s := procname(op(l));
		if s = FAIL then return s fi;
		# Equate leading coefficients of ansatz to entries of s.
		# Here c1 is a new variable to make the system homogeneous.
		eq := {seq(co[i,d[i]] - s[i]*c1, i=1..3)}
	fi;

	# Compute equations for each factor of a*b*c.
	for k to 3 do
		i,j := op({1,2,3} minus {k});
		for l in factors(primpart(v[k],t))[2] do
			s := msqrt(v[i],v[j],l[1],t);
			if s = FAIL then return s fi;
			s := ansatz[i]*s[1]+ansatz[j]*s[2];
			eq := {op(eq),coeffs(collect(rem(s,l[1],t),t),t)}
		od
	od;

	# Solve linear equations, substitute the solution into ansatz.
	s := SolveTools:-Linear(eq, indets(eq) minus indets(v));
	eval(ansatz, s);
	fi
end:


# Compute the square root of -a/b modulo p(t) and return
# a list with the numerator and denominator of this sqrt.
#
msqrt := proc(a,b,p,t)
	local x,R,v;
	if degree(b,t) > degree(a,t) then
		v := procname(b,a,p,t);
		return `if`(v=FAIL,v,[v[2],v[1]])
	fi;
	R := RootOf(p,t);
	v := subs(t = R, primpart(b*x^2+a,x));
	if degree(p,t)=1 then
		v := factors(v)
	else
		# If the number of vars is > 1 then the modular gcd code 
		# for function fields will really help speed things up here
		v := subs(R=t, evala(Factors(v, {R})))
	fi;
	v := select(has,v[2],x)[1][1];
	if degree(v,x)=2 then
		# Square root is non-rational:
		return FAIL
	fi;
	[coeff(v,x,0), coeff(v,x,1)]
end:

# Clear denominators and remove common factors of a list.
#
ReduceList := proc(v)
	local t,i,f;
	f := primpart(add(v[i]*t^i,i=1..nops(v)),t);
	[seq(coeff(f,t,i),i=1..nops(v))]
end:


# Select a good variable appearing in v.
#
SelectVar := proc(v)
	local vars,i,m,mn,mv;
	vars := indets(v);
	if vars = {} then error fi;
	for i in vars do
		m := max(op(map(degree,v,i)));
		if not assigned(mn) or m<mn then
			mn := m;
			mv := i
		fi;
	od;
	mv
end:

# Make a,b,c squarefree and coprime where v = [a,b,c].
#
ReduceConic := proc(v)
	local i,s,t,j,w,k,g;
	for i to 3 do
		s := sqrfree(v[i]);
		t[i] := mul(j[1]^iquo(j[2],2),j=s[2]);
		w[i] := mul(j[1]^irem(j[2],2),j=s[2])*s[1]
	od;
	for k to 3 do
		i,j := op({1,2,3} minus {k});
		g := gcd(w[i],w[j]);
		if g<>1 then
			w[i] := normal(w[i]/g);
			w[j] := normal(w[j]/g);
			w[k] := w[k]*g;
			t[k] := t[k]/g
		fi
	od;
	[t[1],t[2],t[3]], [w[1],w[2],w[3]]
end:




#################
#               #
#   Examples    #
#               #
#################



_Env_LRE_tau := tau;
_Env_LRE_x := x;
read "/Users/heba/Desktop/ORDER4/Order4-Imp/RightFactors.txt": # LREtools[RightFactors] code with a modification in case _Env_RightFactors_data is assigned.


if false then
 L4 := (4*x-11)*(524160*x^8+9391200*x^7-118179432*x^6-253541284*x^5-339259113*x^4-283416626*x^3-140532705*x^2-35130024*x-2220048)*(x+5)^2*(x+4)^2
 *(2*x+9)^2*(2*x+7)^2*tau^4+16*(524160*x^12+15113280*x^11-364158816*x^10-4278491572*x^9-9186978746*x^8+12166953346*x^7-86741410290*x^6-\
 843333775440*x^5-2144077451746*x^4-3001904754612*x^3-2506144851117*x^2-1178353117620*x-242095406175)*(x+4)^2*(2*x+7)^2*tau^3+(-\
 137438945280*x^17-6482503372800*x^16-90355220358144*x^15-154953056569984*x^14+8139627355615616*x^13+69179680108818000*x^12+
 277321791062784832*x^11+698868352509149328*x^10+1236863662787672992*x^9+1625448731323698944*x^8+1626145247262854144*x^7+
 1235819925815197696*x^6+686291085150978048*x^5+244593652122419200*x^4+24045290042818560*x^3-27607241721839616*x^2-15602879836717056*x-\
 2930851407200256)*tau^2-32768*(1048320*x^12+38613120*x^11-475672512*x^10-11499544808*x^9-68147233556*x^8-184773020492*x^7-262836346620*x^6
 -216526023556*x^5-122659285853*x^4-39783078178*x^3+1029344695*x^2+7429526904*x+2484513522)*(x+2)^2*(2*x+3)^2*tau+4096*(4*x+33)*(524160*x^8
 +13584480*x^7-37764552*x^6-736049716*x^5-3014273813*x^4-6181409598*x^3-7122549901*x^2-4430333096*x-1162363872)*(2*x+3)^2*(2*x+1)^2*(x+2)^2
 *(x+1)^2;

 t0 := time();
 _Env_print_number_cases := true;
 ReduceOrder (L4);
 lprint(time() - t0);

 # 24.702 seconds for 1791 cases (use neither)
 # 11.938 seconds for 597 cases  (use DetFactorsSelect)
 # 6.950 seconds for 363 cases   (use DeterminantSelect)
 # 3.185 seconds for 121 cases   (use both)
fi;



if false then
 # OEIS A219670
 L4 := (x + 3)^2*(x + 4)^2*(x + 5)*(2*x + 3)*(7*x^4 + 56*x^3 + 166*x^2 + 216*x + 105)*tau^4
  - (x + 3)*(x + 4)*(2*x + 3)*(2*x + 7)*(70*x^6 + 1050*x^5 + 6406*x^4 + 20337*x^3 + 35449*x^2 + 32244*x + 12048)*tau^3
  - 3*(x + 3)*(2*x + 5)*(490*x^8 + 9800*x^7 + 84910*x^6 + 416150*x^5 + 1261159*x^4 + 2417840*x^3 + 2860095*x^2 + 1905600*x + 546588)*tau^2
  + 27*(x + 2)^2*(2*x + 3)*(2*x + 7)*(70*x^6 + 1050*x^5 + 6406*x^4 + 20283*x^3 + 35044*x^2 + 31221*x + 11178)*tau
  + 729*(x + 1)^3*(x + 2)^2*(2*x + 7)*(7*x^4 + 84*x^3 + 376*x^2 + 744*x + 550);
 # After the improvements, still 5 minutes of CPU time.

 # OEIS A227845
 L4 := (x + 4)^2*tau^4 - 2*(3*x^2 + 21*x + 37)*tau^3 + 2*(3*x^2 + 15*x + 19)*tau - (x + 2)^2;
 ReduceOrder(L4);
fi;


if false then
 L4 := SymProd(tau^2-tau + 2*x, tau^2-tau-x);
 ReduceOrder (L4);
 L4 := SymProd(tau^2+x*tau-1, tau^2-tau-x);
 ReduceOrder (L4);
 ReduceOrder32( tau^3+(x+3)*tau^2-(x+2)*(x+1)*tau-(x+1)^2*x );
 ReduceOrder( tau^4-(x+1)*(2*x+3)*tau^2+(x+1)*(2*x+3)*tau+1/4*(2*x+3)*(2*x+1)*(x+1)*x );
 ReduceOrder( (x+1)*(2*x+5)*tau^4+2*(2*x+3)*x*tau^3+2*(2*x+1)*x*tau^2+(4*x^2+2*x+2)*tau+(2*x-1)*x );
fi;

if false then
 read "Hom.txt":
 # OEIS A247365
 L4 := (16*x^6 + 96*x^5 + 237*x^4 + 307*x^3 + 222*x^2 + 88*x + 15)*tau^4
  + (x + 1)*(64*x^8 + 800*x^7 + 4244*x^6 + 12430*x^5 + 21920*x^4 + 23837*x^3 + 15726*x^2 +  5872*x + 972)*tau^3
  + (-256*x^10 - 3840*x^9 - 25472*x^8 - 98304*x^7 - 244271*x^6 - 408233*x^5 - 464965*x^4 - 357285*x^3 - 178432*x^2 - 53022*x - 7290)*tau^2
  + (-64*x^9 - 800*x^8 - 4244*x^7 - 12422*x^6 - 21868*x^5 - 23719*x^4 - 15610*x^3 - 5847*x^2 - 1010*x - 6)*tau
  + 16*x^6 + 192*x^5 + 957*x^4 + 2535*x^3 + 3765*x^2 + 2977*x + 981;
 ReduceOrder(L4);
 SymProd(op(%[-1]));
 ProjectiveHom(%, L4);

 L4 := symmpower(tau^2+x*tau+2*x+1, 3);
 ReduceOrder(%);
fi: