`LREtools/RightFactors` := module () description
"Factoring Recurrence Operators and Applications, by Mark van Hoeij, www.math.fsu.edu/~hoeij/algorithms (2020).";

  export RightFactors, MinimalRecurrence, GuessRecurrence,
	 toOperator, toRecurrence, LTSing, GenExp;

  local	RFactors, SimplifyR, degreeBound, ExteriorPowerSystem, Subsets, WriteInTermsOfBasis,
	WriteInTermsOf, CommonFactor, ProductShifts, ValueTopCoeffNotSorted, ValueTopCoeff,
	coeffRemember, PolySols, ConvertToOperator, RowReduceR, PluckerRelations, CoeffOfP,
	FindIntShifts, IntShift, PairsAB, MatchAB, CombinationsAB,

	RFactorsHeuristic, NthTerm, DivisionCheck, MinRec, MinRecChar0, MinRecChRem,
	FirstSolution, BuildMatrix, BasSol, RFactorsLCLM, Complete, BasisSolutions,
	SolSubSpace, ApplyR, SelectFactors,

	VarTR, VarT, evalaLcoeff, GenExp_sc, TrimT0N, TrimN, EvalMonomial, EvalMonRem,
	GenExpDelta, ShortForm, ValTerm, TermWithVal, SymprodTn, Tntauk, LowestValuationTerm,
	IntDiff, eNorm, AdjustmentTerm, eTopCoeffs, DetFactors, SubLists, MultiplyLocalDeterminants;


RightFactors := proc(L, u) local v;
	v := indets([_passed], {radical, nonreal});
	if v <> {} then
		v := radfield(v);
		return op(eval( [procname(op(eval([_passed], v[1])))],  v[2]))
	fi;
	v := indets(L,function) minus indets(L,RootOf);
	if type(u, function) then
		toRecurrence(procname(toOperator(L, u), _passed[3..-1]), u)
	elif not has(L, _Env_LRE_tau) and v<>{} and nops(indets(op(0,v[1]),name))=1 then
		procname(L, op(0,v[1])(indets(op(1,v[1]),name)[1]), _passed[2..-1])
	elif member(u, {'DegreeBounds', 'DegreeBound'}) then
		RFactors(L, 0, _passed[3..-1])
	elif u = 'Heuristic' then
		RFactorsHeuristic(L), forget(NthTerm)
	elif u = 'LCLM' then
		RFactorsLCLM(L), forget(NthTerm),forget(BasisSolutions),forget(SolSubSpace)
	elif u = 0 then
		{1}
	else
		RFactors(L, _passed[2..-1])
	fi
end:

# Convert a recurrence to a difference operator in C(x)[tau]
toOperator := proc(L, un, x,tau) local v,i,m,u,n;
	if type(L,`=`) then
		procname(lhs(L) - rhs(L), _passed[2..-1])
	elif _npassed=1 or not type(un,function) then
		v := indets(L,function)[1];
		procname(L, op(0,v)(indets(op(1,v),name)[1]), _passed[2..-1])
	elif _npassed=2 then
		procname(_passed, _Env_LRE_x, _Env_LRE_tau)
	else
		u, n := op(0,un), op(1,un);
		v := {seq(`if`(op(0,i)=u, op(1,i)-n, NULL), i=indets(L,function))}:
		if (not type({u,n},set(name))) or has(v,n) or eval(L, {seq(u(n+i)=0, i=v)}) <> 0 then
			error "Input should be a homogeneous linear recurrence"
		fi;
		m := min(v);
		sort(eval(add(coeff(L,u(n+i)) * tau^(i-m), i=v), n = x-m), tau)
	fi
end:
toRecurrence := proc(L, un::function(name), x, tau) local u,n,i,v;
	if _npassed=2 then
		procname(_passed, _Env_LRE_x, _Env_LRE_tau)
	elif type(L,set) then
		{seq(procname(i,_passed[2..-1]), i=L)}
	elif type(L,list) then	# --> process list output from PluckerRelations:
		[seq(`if`(has(i,tau), procname(i,_passed[2..-1]), i), i=L)]
	else
		u, n := op(0,un), op(1,un);
		v := [seq(i, i = ldegree(L,tau)..degree(L,tau))];
		add(eval(coeff(L,tau,v[-i]), x=n)*u(n+v[-i]), i=1..nops(v)) = 0
	fi
end:

###########################################################################
######         Compute Right-Factors of Recurrence Operators         ######
###########################################################################

# The algorithm starts by computing generalized exponents, defined in
#      www.math.fsu.edu/~hoeij/papers/issac10/J.pdf
# These are used to efficiently compute hypergeometric solutions
# of the d'th exterior power system described in
#      www.math.fsu.edu/~hoeij/slides/bronstein/slides.pdf
# Solutions of this system that satisfy Plucker relations correspond to
# factors of L.
#
# Mark van Hoeij, 2020.

# Input: A recurrence operator L in C(x)[tau].
#        An integer d.
#        (optional): a set ext giving the algebraic extensions in C.
#
# Output: If d > 0:  { Right-Factors of order d }.
#         If d = 0:  [ Degree-bounds for Right-Factors of order 1..n ].
#
RFactors := proc(L, d::nonnegint, ext::set)
	local n, i, M, v, x, tau, pL, P;
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	n := degree(L, tau);
	if L = 0 or n = FAIL then
		error "First argument should be a non-zero polynomial in the shift operator", tau
	elif _npassed < 3 then
		return procname(_passed, indets(L, {'RootOf','radical','nonreal'}))
	elif ldegree(L,tau)<>0 then
		WARNING(cat("Operators are factored in C(",x,")[",tau,",",tau,"^(-1)]. The factor ",tau^ldegree(L,tau)," is ignored as it is a unit in this ring."));
		return procname(`LREtools/mult`(tau^(-ldegree(L,tau)),L), _passed[2..-1])
	elif not type(L, polynom(ratpoly(anything,x), tau)) then
		error "First argument should be polynomial in", tau, "and rational in", x
	elif d > n then
		return {}
	elif d = n then
		return `if`(d=0, [], {SimplifyR(L, x, tau)})
	elif ext = {} then
		Normalizer := normal
	else
		Normalizer := evala
	fi;
	pL := collect(primpart(L,tau), tau);
	P := PairsAB(pL, x, tau, ext);
	if d=0 then
		# compute degree-bounds only but no actual factors
		M := {seq([v[3,2], degreeBound(pL, x, n, v[3,2], v, v[1])], v=P)};
		return [seq(max(seq(`if`(v[1]=i, v[2], NULL), v=M)), i=1..n)]
	fi;
	P := [seq(`if`(v[3,2]=d, v, NULL), v=P)];
	if P = [] then return {} fi:
	# We computed all pairs, then discarded pairs of order <> d. We could have
	# avoided computing pairs of order > d but it makes little CPU difference.
	M := collect(pL/factor(lcoeff(pL,tau)),tau,Normalizer);
	M := ExteriorPowerSystem(M, tau, n, d);
	i := evalb(ext union indets(L) minus {x,tau} = {}); # input over Q?
	if _Env_print_number_cases = true then
		lprint("Number of cases", nops(P))
	fi;
	M := {seq(PolySols(M, x, tau, n, d, v, v[1] * lcoeff(v[3,1],x), i), v=P)};
	forget(coeffRemember); forget(ValueTopCoeff);
	map(SimplifyR, {seq(PluckerRelations(n, d, x, tau, i, ext), i=M)}, x, tau)
end:
SimplifyR := proc(R, x, tau)
	if type(R,list) then return [procname(R[1], x, tau),op(2..-1,R)] fi:
	sort(map(sort, collect(primpart(R*`tools/sign`(R,[tau,x]),tau),tau,factor),x), tau)
end:
degreeBound := proc(L, x, n, d, v, r) local i;
	if v[2] - degree(CommonFactor(n-d, x, r), x) < 0 then # implies lcoeff(R,tau)=0
		-infinity
	elif d < n then
		v[2] - degree(CommonFactor(n-d-1, x, r), x) # Bounds deg(leading coeff)
		     + add(degree(i[2],x), i=v[3,3])        # adds NewtonPolygon
	else
		degree(L,x) # Handle d=n separately because CommonFactor
	fi                  # does not work for n-d-1 = -1
end:


###############################################################################
######                       Exterior Power System                       ######
###############################################################################

# Input: L in C(x)[tau]
# Output: matrix M such that order-d factors of L are hypergeometric solutions of MY = tau(Y).
ExteriorPowerSystem := proc(L, tau, n, d)
	local B, b;
	B := Subsets(n,d);  # if b = [1,3,4] in B then view it as tau^0 wedge tau^2 wedge tau^3
	Matrix([seq(WriteInTermsOfBasis(L, tau, b, B, n), b=B)]) # = the action of tau on B.
end:
# Subsets of {1..n} given as ordered lists
Subsets := proc(n::nonnegint, d::nonnegint) local i;
	options remember;
	if d=0 then
		[ [ ] ]  # the only list with 0 elements is [ ]
	elif n=0 then
		[ ]
	else
		[ op(procname(n-1, d)), seq( [op(i),n], i=procname(n-1, d-1)) ]
	fi
end:
WriteInTermsOfBasis := proc(L, tau, b, B, n)
	local i, tau_b, j;
	tau_b := [seq(i+1, i=b)];	# b = tau^i1 wedge ... So tau_b = tau^(i1+1) wedge ...
	if tau_b[-1] <= n then
		WriteInTermsOf(tau_b, B, 1)
	else
		tau_b := tau_b[1..-2];
		add(`if`(member(i, tau_b), 0, WriteInTermsOf(sort([op(tau_b), i]), B,
		 -coeff(L, tau, i-1)/coeff(L, tau, n) * mul(`if`(j>i,-1,1), j=tau_b))), i=1..n)
	fi
end:
WriteInTermsOf := proc(b, B, c) local i;
	for i do if B[i]=b then return [0$(i-1), c, 0$(nops(B)-i)] fi od
end:


###############################################################################
######    Solving M Y = r tau(Y) where M is the Exterior Power System    ######
###############################################################################


CommonFactor := proc(s, x, r) options remember, system;
	denom( ProductShifts(eval(r, x=x-s), x, s) )
	# lcm(%, procname(s-1,x,r)) if s>0 is also a valid bound,
	# but it turns out to be the same.
end:
ProductShifts := proc(r, x, n) options remember, system; local i;
	Normalizer( mul(eval(r, x=x+i), i=0..n-1) )
end:

ValueTopCoeffNotSorted := proc(S, x, tau, d, G, g)
	local k;
	if member(S[-1],S[1..-2]) then return 0 fi;
	k := 1;
	while S[k] < S[-1] do k := k+1 od;
	(-1)^(d-k) * ValueTopCoeff([op(S[1..k-1]),S[-1],op(S[k..-2])], x, tau, d, G, g)
end:
ValueTopCoeff := proc(S, x, tau, d, G, g) options remember; local i,j,DB;
	if max(S) = d then
		return 1
	fi;
	DB := floor( -add((S[j]-j)*g[j], j=1..d) );
	coeff(-add(coeffRemember(G,x,tau,i) * ValueTopCoeffNotSorted([op(S[1..-2]), S[-1]-d+i],
		x, tau, d, G, g), i=0..d-1), x, DB) * x^DB
end:
coeffRemember := proc(G,x,tau,i) options remember; # (table cleared at the end of RFactors)
	local j,L;
	if _npassed=3 then
		L := mul(j[2], j=G);
		collect(L/lcoeff(L,tau), tau, normal)
	else
		L := coeff(procname(G,x,tau),tau,i);
		lcoeff(L,x) * x^degree(L,x)
	fi
end:

# Solve M Y = r tau(Y) where M is an exterior power matrix.
PolySols := proc(M, x, tau, n, d, v, r, overQ)
	local a, j, Ansatz, CF, DB, G, S, S1, Sn, vars, g, i, isz, l, s, sb,
	 top1, SB, SBN, ANS, vl, Y, MAT, e, c, varY, npiv, piv, sol;

	G := sort([op(v[3,3])], (i1,i2) -> evalb(i1[1] > i2[1]));
	# Example: [[-1/2, tau^2-x], [-2/3, tau^3-x^2]]
	g := [seq(i[1]$degree(i[2],tau), i=G)];

	i := 0; SB := {};
	Sn := Subsets(n,d);
	for S in Sn do if member(1,S) then
		i := i+1;
		CF[S] := CommonFactor(n-max(S), x, r);
		DB[i] := v[2] - add((S[j]-j)*g[j], j=1..d) - degree(CF[S],x);
		sb := {};
		if i = 1 then
			if DB[1] < 0 then
				return NULL # DegreeBound DB[1] is negative
			fi;
			userinfo(3,'LREtools', "Building equations for type", tau-r,
			  "with top coefficients", mul(j[2],j=G),
			  "and component-wise Content Bound", [seq(CommonFactor(n-j, x, r), j=d..n)]);
			top1 := a[i, DB[i]] * lcoeff(CF[S], x);
		else
			DB[i] := floor(DB[i]);
			sb := a[i,DB[i]] = lcoeff(ValueTopCoeff(S,x,tau,d,G,g), x) * top1 / lcoeff(CF[S],x);
			SB := SB union {lhs(sb)-rhs(sb)}
		fi;
		Ansatz[S] := subs(sb, add(a[i,j] * x^j, j=0..DB[i])) # Part of Ansatz without common factor
	fi od;
	for S in Subsets(n,d) do if not member(1,S) then
		s := min(S)-1;
		S1 := [seq(j-s,j=S)];	# Ansatz[S] is a shift of Ansatz[S1]
		CF[S] := Normalizer(ProductShifts(r, x, s) * eval(CF[S1], x = x+s));
		Ansatz[S] := subs(x=x+s, Ansatz[S1]) # don't Normalize this
	fi od;
	Ansatz := [seq(CF[S] * Ansatz[S], S = Sn)];
	# Further improvements to "content bound" CF[S] may not be needed as it
	# tends to be close to content( ConvertToOperator(Ansatz,tau,d), tau).

	# Put primary variables a[1..d, *] at the end (ordering affects CPU time)
	g := [seq(i-j, j=0..i-1)];
	vars := [seq(seq(a[i,DB[i]-j], j=`if`(i=1,0,1)..DB[i]),i=g)];
	c := nops(vars);

	# Replacing many variables a[i,j] by 1 variable Y speeds polynomial arithmetic:
	varY := {seq(vars[i] = Y^i, i=1..c)}:
	l := Vector( eval(Ansatz, varY) );
	isz := M . l - r * subs(x=x+1, l); # should be zero

	g := nops(Sn); # Reverse ordering of isz and take coefficients:
	isz := map(coeffs,[seq(`if`(member(n, Sn[g-i]),
		sort(collect(numer(isz[g-i]),x),x),
		NULL), # Discard entries that are zero by construction
	i=0..g-1)], x);

    if overQ then
	MAT := Matrix(nops(isz),c);
	i := 0;
	for g in isz do e := expand(g); if e<>0 then
		i := i+1;
		for j in `if`(type(e,`+`),e,[e]) do
			MAT[i,degree(j,Y)] := eval(j, Y=1)
		od;
	fi od:
	MAT, piv, npiv := SolveTools:-LinearSolvers:-RationalDenseNullSpace(MAT, ':-skiptest'=true);
	if npiv = [] then return NULL fi;
	sol := {seq(vars[i] = add(MAT[i,j] * vars[npiv[j]], j=1..nops(npiv)), i = piv)};
	l := eval(Ansatz, sol);
	i := convert(l, Vector);
	# Skipped test in NullSpaceOverQ, to prove correctness check that this is zero:
	MAT := M . i - r * subs(x=x+1, i);
	MAT := {seq(`if`(member(n, Sn[i]), normal(MAT[i]), NULL), i=1..nops(Sn))} minus {0};
    fi;
	if MAT <> {} then
		userinfo(1,'LREtools',`if`(overQ, "Check failed, redoing solve", "Input is not over Q"));
		sol := subs([seq(Y^(c-i) = vars[c-i], i=0..c-1)], map(collect,{op(isz)},Y));
		sol := SolveTools:-Linear(sol, vars);
		l := eval(Ansatz, sol);
	fi;
	ANS := NULL;
    do
	while v[2] - degree(l[1],x) > 0 do
		if l[1] = 0 then
			return ANS # stop because order(l) < d
		fi;
		SBN := eval(SB, {seq(i = a[op(1,i),op(2,i)-v[2]+degree(l[1],x)],
		  i = indets(SB) minus indets([_passed]))});  # Not indets intersect vars
		if Normalizer(eval(SBN, sol)) minus {0} = {} then
			break # SBN holds
		fi;
		# Require SBN to hold so that top coeffs match GenExp data.
		# That reduces dim(output), reducing EQ in PluckerRelations.
		sol := SolveTools:-Linear(sol union SBN, vars);
		l := eval(l, sol)
	od;
	vl := {op(vars)} intersect indets(l);
	# lprint("Number of solutions", nops(vl));
	if vl = {} then
		return ANS
	else
		ANS := ANS, [seq( map(coeff, l, i), i=vl)];
		if nops(vl)=1 then return ANS fi:
		sol := SolveTools:-Linear(sol union {lcoeff(l[1],x)}, vars);
		l := eval(l, sol);
		# In case the leading coefficient is not normalized:
		while l[1]<>0 and Normalizer(lcoeff(l[1],x))=0 do
			l := subs( lcoeff(l[1],x) = 0, l)
		od:
	fi
    od
end:


################################################################################
######    Use Plucker Relations to select solutions that give RFactors    ######
################################################################################


# Converts a solution v of the exterior power system to an operator.
ConvertToOperator := proc(v, tau, d) local i;
	add( (-1)^i * v[i+1] * tau^(d-i), i=0..d)
end:
# Input: R is a linear combination of vars.
# Output: change-of-basis that typically makes R smaller.
RowReduceR := proc(R, x, tau, d, vars)
	local cc,i,treated,sb,c,cf,v,tmp;
	if has(R, tau) then
		cc := [lcoeff(R,tau), seq(coeff(R,tau,i), i=0..d-1)]
	else
		cc := R
	fi;
	treated := {};
	sb := NULL;
	for c in cc do for i from degree(c,x) to 0 by -1 while treated<>vars do
		cf := subs(sb, coeff(c,x,i));
		v := (indets(cf) intersect vars) minus treated;
		if v={} then next fi;
		sb := sb, subs(tmp=v[1], SolveTools:-Linear({cf = tmp}, {v[1]}));
		treated := treated union {v[1]}
	od od;
	sb
end:

# Input: n, d, V = list-of-poly-sols
PluckerRelations := proc(n, d, x, tau, V, ext)
	local R, i, vars, P, B, EQ, l, s, Y, sb, N, RI;
	N := nops(V);
	R := add(_Z[i] * ConvertToOperator(V[i],tau,d), i=1..N);
	R := primpart(R, {tau,seq(_Z[i],i=1..N)});
	EQ := {};
	if d > 1 and d < n-1 then
		# Convert lists to polynomials so we can divide by the gcd:
		P := add(collect(add(_Z[i] * V[i][l], i=1..N), x) * Y^l, l=1..nops(V[1]));
		P := primpart(P, Y); # In ex2 this is faster
		B := Subsets(n,d);
		# Now construct binomial(n-1, d)-1 quadratic equations for the _Z[i]
		EQ := {seq(seq(add(
			  eval(coeff(R, tau, l-i), x=x+i) * CoeffOfP(s, l, B, P, Y),
			  l=i..d+i), s = Subsets(i+d-1, d-2) )	# s wedge tau^i R = 0
			,i=1..n-1-d)};
		EQ := {seq(coeffs(collect(i,x),x), i=EQ)}
	fi;
	vars := {seq(_Z[i],i=1..N)};
	sb := RowReduceR(R, x, tau, d, vars); # Change-of-basis to make R smaller.
	l := indets(subs(sb,lcoeff(lcoeff(R,tau),x))) intersect vars;
	if nops(l) <> 1 then error "Due to RowReduceR this should have only one variable" fi;
	sb := sb, {l[1] = _Z[N], _Z[N] = l[1]};
	R := collect(subs(sb,R), tau, Normalizer);
	RI := subs(_Z[N] = 1, R);    # Inhomogeneous R and
	vars := vars minus {_Z[N]};  # its variables
EQ := map(numer,map(evala, subs(sb, _Z[N]=1, EQ))) minus {0};
	if EQ = {} then
		return `if`(N=1, RI, [RI, "a family of dimension", N-1, "with equations", {} ])
	fi;
	EQ := Groebner['Basis'](EQ, tdeg(op(vars)));
	if member(1,EQ) then
		NULL
	else
		l := {seq(`if`(degree(i,vars) = 1, i, NULL), i=EQ)};
		if l <> {} then
			l := SolveTools:-Linear(l, vars);
			RI, EQ := subs(l, RI), subs(l, EQ);
			vars := indets(map(rhs,l)) intersect vars;
			EQ := Groebner['Basis'](EQ, tdeg(op(vars)))
		fi;
		N := nops(vars);
		RI, EQ, vars := op(subs( {seq(vars[i] = _Z[i], i=1..N)}, [RI, EQ, vars]));
		if N = 0 then
			RI
		elif N = 1 and has(EQ,_Z[1]) then # Solving a univariate quadratic equation:
			seq(eval(RI, i), i=solve(evala(map(Factor,EQ,ext)), vars))
		else
			[RI, "a family of dimension", Groebner:-HilbertDimension(EQ,vars),
			 "with equations", {op(EQ)}]
		fi
	fi
end:
# s wedge l corresponds to an element of basis B.
# Find that element and return the corresponding coefficient of P.
#
CoeffOfP := proc(s, l , B, P, Y) local i,b,sg;
	if member(l, s) then
		0
	else
		sg := mul(`if`(i > l, -1, 1), i=s); # = sign of sort-reordering:
		b := sort( [op(s),l] );
		b := [1, seq(i+1, i=b)]; # Using only wedge products that include tau^0
		for i do if B[i]=b then return sg * coeff(P, Y, i) fi od
	fi
end:


##################################################################################
###     Construct r = c * A/B for which M Y = r tau(Y) could have PolySols.    ###
###     Pairs A,B come from trailing/leading coefficient while c and degree    ###
###     bound come from Generalized Exponents.                                 ###
##################################################################################


# Example:
#	A = [[x+1/2,1], [x,2], [x-1,1]];
#	P = x
# Then Output = [[x + 1/2, 1]],
#		[[x, 0], [x, 0], [x - 1, 1]]
#
FindIntShifts := proc(A, P, x) local res, resA, i, d;
	resA := NULL; res := NULL;
	for i in A do
		if IntShift(P, i[1], x, 'd') then
			res := res, [i[1],d]$i[2]
		else
			resA := resA, i
		fi
	od;
	[resA], sort([res],proc(i,j) evalb(i[2]<j[2]) end)
end:
# Check if two monic irreducible poly's are an integer shift from each other.
IntShift := proc(A, P, x, d) local n,e;
	n := degree(A,x);
	if n=degree(P,x) then
		e := Normalizer(coeff(A,x,n-1)-coeff(P,x,n-1))/n;
		if type(e,integer) and Normalizer(A-subs(x=x+e,P))=0 then
			if _npassed=4 then
				d := e
			fi;
			return true
		fi
	fi;
	false
end:

# Leading and Trailing Singularities after some desingularization.
#
LTSing := proc(L, x, tau) local i, l, d, M;
	l, d := ldegree(L,tau), degree(L,tau);
	i := iquo(l+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);
	[gcd(coeff(L,tau,d), subs(x=x-1, coeff(M,tau,d+1))), gcd(coeff(L,tau,l), coeff(M,tau,l))];
end:
# Faster version
#
LTSing := proc(L, x, tau, ext) local M,l,d,m,c,A,i,r, c1,c0;
	l, d := ldegree(L,tau), degree(L,tau);
	A := seq(`LREtools/g_p/factors`(coeff(L,tau,c),x,ext), c=[d,l]);
	if ext union (indets(L) minus {x,tau}) <> {} then return A fi:
	m := iquo(l+d, 2);
	c1 := coeff(L, tau, m+1);  c0 := subs(x=x+1, coeff(L, tau, m));
	M := c1 * subs(x=x+1, collect(L*tau,tau)) - c0 * L;	# (c1 tau - c0) * L
	for c to 2 do for i to nops(A[c]) do
	    m[c,i] := A[c,i,2];
	    if degree(A[c,i,1],x)=1 then
		r := solve(A[c,i,1],x) - `if`(c=1,1,0);
		if eval(M,x=r) = 0 and Normalizer(eval([c1,c0][c], x=r))<>0 then
		    m[c,i] := `if`(m[c,i]=1, 0, m[c,i] - ldegree(convert(subs(x=x+r,
				series(M, x=r, m[c,i])),polynom), x));
		    userinfo(2,'LREtools',cat("Reduced ",
			`if`(c=1,"lead","trail"),"ing singularity"),A[c,i,1])
		fi
	    fi
	od od:
	seq([seq(`if`(m[c,i]>0, [A[c,i,1], m[c,i]], NULL), i=1..nops(A[c]))], c=1..2)
end:


# Look for pairs A,B such that we have to find solutions of M the form
#  	Sol(tau - c * A/B) * Polynomial
#
PairsAB := proc(L, x, tau, ext) local Af, Bf, i, Pseq, P, vA, vB, w, dP, resP, s;
	options remember, system;
	Pseq := NULL;
	Bf, Af := LTSing(_passed);	# Example:
					#	Af = [[x+1/2,1], [x,2], [x-1,1]]
	do				#	Bf = [[x-5/2,2], [x,1], [x-2,1], [x+2,1]]
		# Select the next factor P:
		i := [op(Af), op(Bf)];
		if i = [] then break fi;
		P := i[1,1];
		Af, vA := FindIntShifts(Af, P, x);
		Bf, vB := FindIntShifts(Bf, P, x);

		# Example: If P was x, then we now have
		# 	Af = [[x + 1/2, 1]]
		#	Bf = [[x - 5/2, 2]]
		#	vA = [[x, 0], [x, 0], [x - 1, 1]]
		#	vB = [[x + 2, -2], [x, 0], [x - 2, 2]]
		#
		# For A, use the last ones from vA in resP.
		# For B, use the first ones from vB in resP (these are moved from vB to w)

		w := []; resP := [vA,w];

		# Now loop over the nops(vA)+nops(vB)+1 combinations:
		do
			if [vA,vB]=[[],[]] then
				break
			elif vA = [] or (vB<>[] and vA[1,2] > vB[1,2]) then
				w := [op(w), vB[1]];
				vB := vB[2..-1]
			else
				vA := vA[2..-1]
			fi;
			resP := resP, [vA,w]
		od;

		dP := degree(P,x);
		# Give list of all [A,B,s,d] and add this list to Pseq
		# Here s corresponds to a slope of the NewtonPolygon
		# So s = -v where v is as in AAECC, Vol. 17, 83-115 (2006).
		# The d is the same as in that paper (e.g. a polynomial solution
		# of degree 5 would correspond to some [c,s,d] with c=1, s=0, d=5.
		# Here [c,s,d] corresponds to a factor tau - c*x^s*(1+d*t+O(t^2))
		# where t = 1/x.  In other words:   tau - c*t^(-s)*(1+d*t+O(t^2))
		Pseq := Pseq, [seq([mul(i[1],i=w[1]), mul(i[1],i=w[2]),
			(nops(w[1])-nops(w[2])) * dP, Normalizer(
				 add(coeff(i[1], x, dP-1), i=w[1])
				-add(coeff(i[1], x, dP-1), i=w[2]))], w=[resP])]
	od;
	# Check p-curvature? It'd have to be fast (p = 2 or 3) for it to make sense.
	w := DetFactors(L,x,tau,ext); forget(EvalMonRem);
	if assigned(_Env_RightFactors_data) then
		w := {seq(`if`(DetFactorsSelect(i, x, tau), i, NULL), i=w)}
	fi;
	s := {seq(degree(i[1],x),i=w)};
	{seq(`if`(member(P[3],s), seq(MatchAB(op(P),i,x),i=w), NULL), P=CombinationsAB(Pseq))}
end:
# Check if an A/B can be matched.
#
MatchAB := proc(A,B,s,d,L,x) local db;
	if s=degree(L[1],x) then
		db := Normalizer( coeff(L[1],x,s-1)/lcoeff(L[1],x)-d);
		if type(db, nonnegint) then
		    if not assigned(_Env_RightFactors_data) or DeterminantSelect(A/B, x) then
			return [A/B, db, L]  # db = degree-bound for lcoeff(R, tau)
		    fi
		fi
	fi;
	NULL
end:
CombinationsAB := proc(L) local i,j;
	if _npassed=0 then
		[ [1,1,0,0] ]
	elif _npassed=1 then
		L
	else
		[seq(seq([i[1]*j[1], i[2]*j[2], i[3]+j[3], i[4]+j[4]],
		 j=L),i=procname(_passed[2..-1]))]
	fi
end:

#############################################################################
######      Fast heuristic, produces a subset of the right-factors     ######
#############################################################################

# Mark van Hoeij, 2020.

RFactorsHeuristic := proc(L0) local x,tau,  P, A, ANS, L, Ord, db, i, lc, r;
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	P := collect(primpart(L0, tau), tau);
	Normalizer := `if`(indets(P, {'RootOf','radical','nonreal'}) = {}, normal, evala);
	Ord := degree(P, tau);
	db := RFactors(P, 0);
	ANS := {};
		# P -> subs(x=-x,tau=tau^(-1), P) is an automorphism of C(x)[tau, tau^(-1)]
  for L in [P, collect(subs(x=-x,tau=tau^(-1), P) * tau^Ord,tau)] do
	lc := {seq(i[1]+Ord, i=roots(lcoeff(L,tau),x))};
	while lc <> {} do
		r := lc[1]; lc := lc minus {r};
		if member(true,[seq(type(Normalizer(i-r),posint),i=lc)]) then
			next
		fi;						# Heuristic Factoring:
		A := [L,x,tau,Ord,r,0];				# NthTerm plus list A
		NthTerm(0, A) := 1;				# define a solution u of L.
		for i to Ord do					# MinRec below computes the
			NthTerm(-i, A) := 0			# minimal operator r of u.
		od;						# If r <> L then we found
		r, i := MinRec(L, Ord, x, tau, A,r, db);	# a non-trivial factor.
		if r<>L then
			i := degree(r,tau);
			ANS := ANS union {SimplifyR(`if`(P=L, r,
				collect(subs(x=Ord-i-x,tau=tau^(-1),r)*tau^i, tau)),x,tau)}
		fi			# reversing the -x,tau^(-1) automorphism
	od
  od;
	ANS
end:


###############################################################################
####     MinRec: Minimal Operator of a sequence given by NthTerm(n, A)     ####
###############################################################################

# Compute the n'th term of a recurrence given by a table A.
#
NthTerm := proc(n::nonnegint,A)	# A = [L, x, tau, Ord, offset]
	options remember;	# remember table is cleared when done.
	local m,i;
	m := n-A[4];
	Normalizer(-add(eval(coeff(A[1],A[3],i   ), A[2] = A[5]+m) * procname(m+i,A), i = 0..A[4]-1)
		      / eval(coeff(A[1],A[3],A[4]), A[2] = A[5]+m) )
end:

# To see if R divides L check if `LREtools/rightdivision`(L,R)[2] = 0.
# This is a fraction-free mod-p version of that.
#
DivisionCheck := proc(L,R,x,tau,p) local l,r,d,rd,q;
	l := traperror( `mod/Primfield/ReduceField`([L,R], {x,tau}, p) );
	if l = lasterror or map([degree,ldegree],l,tau) <> map([degree,ldegree],[L,R],tau) then
		userinfo(1,'LREtools',"Full division check");
		return evalb(Normalizer(`LREtools/rightdivision`(L,R)[2]) = 0)
	fi;
	l,r := op(l);
	d := degree(l,tau) - degree(r,tau);
	while d >= 0 do
		rd := subs(x=x+d, r);
		q := Normal(lcoeff(l, tau)/lcoeff(rd, tau)) mod p;
		l := Primpart(denom(q)*l - numer(q) * rd * tau^d, tau) mod p;
		d := degree(l,tau) - degree(r,tau);
	od;
	userinfo(`if`(l=0,2,1),'LREtools',cat("Division check ",`if`(l=0,"passed","FAILED")));
	evalb(l = 0)	# Note: Primpart mod p already reduces algebraic extensions.
end:

#  Input is L in C(x)[tau], its order, degree-bound for factors, and
#  a table A with data to compute the NthTerm of the sequence.
#  Output is the minimal order operator that the sequence satisfies
#  for all  n >= min_n  for some finite min_n.
#
#  To ensure termination, only use a proven degree-bound (key = generalized exponents).
#
MinRec := proc(L, Ord, x, tau, A,offset, DegBound) local tc, i, R, min_n, ord, overQ;
	userinfo(2,'LREtools',"Degree bounds", DegBound);
	tc := coeff(L,tau,0);
	if offset <> 0 then
		tc := collect(subs(x = x + offset, tc), x, Normalizer)
	fi;
	min_n := max(0, seq(`if`(type(i[1],integer), i[1]+1, NULL), i=roots(tc, x)) );
	overQ := evalb({} = indets(L,{'RootOf','radical','nonreal'}) union indets(L) minus {x,tau});
	for ord from 0 to Ord-1 do
		if _Env_LRE_GuessR = true and min_n > 2500 then break fi;
		tc := _passed[3..6], min_n,ord, `if`(ord=0,0,DegBound[ord]);
		if tc[-1] < 0 then
			next
		elif overQ then
			R := traperror(MinRecChRem(tc)); # trapping(modular inverse does not exist)
			if R=lasterror then
				userinfo(1,'LREtools', R, "Switching to Char0 code.")
			elif R = 0 then
				next
			elif DivisionCheck(L,R,x,tau, modp1(Prime(1))) then
				return R, min_n
			elif _Env_LRE_GuessR = true then
				next	# Don't call slow Char0 code under GuessRecurrence
			fi
		fi;
		R := MinRecChar0(tc);
		if R <> 0 and DivisionCheck(L,R,x,tau, modp1(Prime(1))) then
			return R, min_n
		fi
	od;
	L, min_n	# L already had minimal order.
end:
MinRecChar0 := proc(x,tau, A,offset, min_n,ord,db)
	local a, i, j, m, R, eq, ne, s, var;
	# Make an ansatz for the right-factor R based on the degree bound:
	R := add(add(a[i,j] * x^i, i=0..db) * tau^j, j=0..ord);
	var := indets(R) minus {x,tau};
	s := min_n;	# The part of the sequence we use is: u(n+offset) for n >= min_n
	while var <> {} do
		ne := nops(var) + 4; # Taking 4 more equations than the minimum.
		eq := evala({seq(add(eval(coeff(R,tau,m), x = offset+i) * NthTerm(i+m,A),
			m = 0..ord), i = s .. s + ne-1)});
		if s > min_n and eq = {0} and _Env_LRE_GuessR = true then
			# GuessRecurrence uses _Env_LRE_GuessR to avoid bad recurrences
			# but that can interfere with termination. Fix:
			eq := {lcoeff(lcoeff(R,tau),x)}
		fi;
		eq := SolveTools:-Linear(eq, var);
		R := evala(Primpart( subs(eq, R), tau));
		s := s + ne;
		var := indets(R) intersect var
	od;
	R
end:
MinRecChRem := proc(x,tau,A,offset,min_n,ord,db)
	local c,r,i,j;
	c := (db+1)*(ord+1);	# c unknowns and c+8 equations.
	r := c+8;
	j := [seq(NthTerm(i,A), i = min_n .. min_n + ord+r-1)];
	r := FirstSolution(j, r,c, _passed[4..7]);
	if r=0 then
		0	
	else
		add(add(r[i*(ord+1)+j+1] * x^i * tau^j, i=0..db), j=0..ord)
	fi
end:
# Compute the First element of a NullSpace basis:
#
FirstSolution := proc(terms, r,c, offset,min_n, ord,db)
	local m,k,d,p,A,rnk,bp,i,s,f,N,lc,rc,x;
	userinfo(3,'LREtools',cat("Trying order ",ord," and degree ",db));
	m := nops(terms);
	rc := `if`(_Env_LRE_GuessR = true, ceil(m/3) .. m-1, {});
	m := 1;
	k := kernelopts(wordsize);
	d := 'integer'[k/8];
	p := 2^(k/2 - 4);	# Replacing NthTerm with NthTermModP did not make it faster.
	do
		p := nextprime(p);
		A := LinearAlgebra:-Modular:-Mod(p,BuildMatrix(_passed[2..7],terms mod p, p),d);
		# Only do a full Row Reduction if m > 1 (then rank < c is likely)
		LinearAlgebra:-Modular:-RowReduce(p,A,r,c,c,0,0,'rnk',0,0, evalb(m > 1));
		if rnk = c then
			userinfo(4,'LREtools',"No solution");
			return 0
		elif m=1 then
			# NullSpace non-trivial so we need a full Row Reduction:
			LinearAlgebra:-Modular:-RowReduce(p,A,r,c,c,0,0,0,0,0,true)
		else
			bp := false;
			for i to f-1 do if A[i,i]=0 then bp := true fi od;
			if bp then
				next	# discard p
			elif A[f,f] <> 0 then
				m := 1	# discard all previous primes
			fi
		fi:
		for f while A[f,f]<>0 do od;		# First free variable
		A := LinearAlgebra:-Column(A, f);
		A[f] := -1;
		if m=1 then
			N := A; lc := rc
		else
			N := chrem([N,A], [m, p]);
		fi;
		if lc <> {} then
			s := 0; k := ord+1;
			while s=0 do
				s := add(A[i*(ord+1)+k]*x^i, i=0..db); k := k-1
			od;
			lc := {seq(`if`(eval(s, x=i+offset-k) mod p = 0, i, NULL), i=lc)}
		fi;
		m := m * p;
		if lc = {} then		# MinRec/DivisionCheck will check correctness:
			A := iratrecon(N, m, 'scaled');
			if A <> FAIL and length(m) > 3 + max(map(length,A)) then
				userinfo(3,'LREtools',"Solved"); return A
			fi
		elif m > p then		 # _Env_LRE_GuessR ==> don't return a bad recurrence.
			userinfo(4,'LREtools',
			  "Discarding a recurrence that does not determine higher terms");
			return 0
		fi;
		if assigned(_Env_LRE_GuessM) and length(m) > _Env_LRE_GuessM then
			return 0	# _Env_LRE_GuessM => no large coeff size.
		fi
	od				
end:
BuildMatrix := proc(r,c, offset, min_n, ord, db, vals, p) local e,a,a1,i,j,M;
	a[0] := 1;
	M := Matrix(r,c);
	for e to r do
		a1 := min_n + offset + e - 1;
		if _npassed = 8 then a1 := a1 mod p fi:
		for i to db do
			a[i] := a1 * a[i-1];
			if _npassed = 8 then a[i] := a[i] mod p fi:
		od;
		for i from 0 to db do for j from 0 to ord do
			M[e, i*(ord+1)+j+1] := a[i] * vals[e+j]
		od od
	od;
	M
end:


#######################################################################
######         Guess Recurrence for a list of terms              ######
#######################################################################

# Similar to gfun[listtorec] but this can find larger recurrences.
#
GuessRecurrence := proc(vals::list, un::function, offset)
	local x,tau, i,j,  A, c, db, ex, n, ord, overQ, r;
	if _npassed = 2 or offset = 'Minimize' then return
		procname(_passed[1..2], 0, _passed[3..-1])
	fi;
	ord := 0; n := 0;
	for i in vals do
		if i=0 then
			n := n+1; ord := max(n, ord)	# Initial order should be more than
		else					# the number of consecutive 0's
			n := 0				
		fi
	od;
	n := nops(vals);
	if n = ord then return un = 0 fi;
	_Env_LRE_GuessR := true;
	overQ := evalb(type(vals, list(rational)));
	do
		ord := ord + 1;
		ex := max(2, floor(n/15));
		db := floor( (n-2*ord-1-ex)/(ord+1) );
		if db < 0 then return FAIL fi;
		ord := floor( (n-db-ex-1)/(db+2) );
		c := (db+1)*(ord+1);
		r := n - ord;
		if r - c < ex then error "unexpected" fi;
		if overQ then
			A := FirstSolution(vals, r,c,offset,0,ord,db)
		else
			A := LinearAlgebra:-NullSpace(BuildMatrix(r,c,offset,0,ord,db,vals));
			if A={} then A := 0 else A := A[-1] fi
		fi;
		if A <> 0 then
			r := add(add(A[i*(ord+1)+j+1] * x^i * tau^j, i=0..db), j=0..ord);
			if _passed[-1]='Minimize' and has(r,x)	# has(r,x): Don't "simplify"
			and ldegree(r,tau) = 0 then		# a constant recurrence.
				r := evala(Primpart(r,tau));	# ldegree > 0 ==> bad initial
				j := degree(r,tau);		# conditions were given.
				A := [r,x,tau,j,offset,1];
				for i to nops(vals) do
					NthTerm(i-1, A) := vals[i]
				od;
				_Env_LRE_x, _Env_LRE_tau := x,tau;
				_Env_LRE_GuessM := 30 + 2 * max(map(length,[coeffs(r,{x,tau})]));
				i := RFactors(r,0);		# don't let coeffs grow too much
				if max(i) < n then
					r := traperror(MinRec(r,j,x,tau,A,offset,i)[1])
				fi;
				if r = lasterror then
					if r[20..27] = "division" then
						next	# We found a bad recurrence.
					else		# Maybe "return FAIL" is better?
						error r
					fi
				fi
			fi;
			return toRecurrence(SimplifyR(r,x,tau), un, x,tau)
		fi
	od
end:

#####################################################################
######       Reduce Recurrence to proven Minimal Order         ######
#####################################################################

MinimalRecurrence := proc(R, un, init, V)
	local u,n, x, tau, A, L, Ord, ext, i, j, iv, lc, m, offset, r, B,N, Range;
	if _npassed=1 or not type(un, function(name)) then
		u, n := {seq(op(0,i),i = indets(R,function))}, indets(R,name);
		if nops(u)*nops(n) <> 1 then
			error "Second argument should be of the form u(n) for some names u,n."
		else
			return procname(R, u[1](n[1]), _passed[2..-1])
		fi
	fi;
	u, n := op(0,un), op(1,un);
	if not (type(u,name) and type(n,name)) then
		error "Second argument should be of the form u(n) for some names u,n."
	elif type(R,set) then
		L := [seq(`if`(has(i,n),i,NULL), i=R)];
		if nops(L)<>1 then
			error "First argument should contain one linear homogeneous recurrence"
		fi;
		return procname(L[1], un, R minus {L[1]}, _passed[3..-1])
	elif not has(R,u) then
		error "dependent variable does not appear in the recurrence"
	fi;
	# Now check if the sequence is defined for all n
	iv := {seq(`if`(op(0,i)=u, i, NULL), i=indets(init,function))};
	if not type(init,set(equation)) or iv <> map(lhs,init) or has(init,n) then
		error "Third argument should be a set of initial terms {u(0)=.., u(1)=.., ...}."
	fi;
	offset := min(seq(op(1,i), i=iv));
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	ext := indets([_passed], {'RootOf','radical','nonreal'});
	Normalizer := `if`(ext = {}, normal, evala);
	L := toOperator(R, u(n), x,tau);
	Ord := degree(L, tau);
	m := [seq([i, Normalizer(subs(x=offset+i, init, add(coeff(L,tau,j)*u(offset+i+j),
		j=0..Ord)))], i=0..nops(init)-Ord-1)];
	m := [seq(`if`(i[2]<>0 and not has(i[2],u), i[1]+offset, NULL),i=m)];
	if m <> [] then WARNING(cat(
		"Discarding some initial conditions; Recurrence does not hold for ", n=op(m)));
		m := max(m);
		return procname(R, un, {seq(`if`(op(lhs(i))>m,i,NULL),i=init)}, _passed[4..-1])
	fi;
	if {seq(u(offset+i), i=0..Ord-1)} minus iv <> {} then
		error "Need at least as many initial terms as the order of the recurrence."
	fi;
	lc := `if`(offset=0, lcoeff(L,tau), evala(Expand(subs(x = x + offset, lcoeff(L,tau)))));
	m := min(seq(`if`(type(i[1],nonnegint) and not member(u(i[1]+Ord+offset), iv),
		i[1], NULL), i=roots(lc,x)));
	if m<>infinity then
		error "Recurrence + initial terms do not determine", u(m+Ord+offset)
	fi;
	A := [L,x,tau,Ord,offset,2];
	for i in init do
		NthTerm(op(1,lhs(i))-offset, A) := rhs(i)
	od:
	r, m := MinRec(L, Ord, x, tau, A,offset, RFactors(L, 0, ext));
	N := degree(r,tau);
	if _npassed > 3 and _passed[-1] = `+` then

	####################################################################################
	###      SumDecompose.  Write u(n) as u1(n) + u2(n) + .... of minimal order      ###
	####################################################################################

		# Key step: Write the Minimal Order Operator r as an LCLM:
		L := [op(RFactorsLCLM(r))];
		forget(BasisSolutions); forget(SolSubSpace);
		m := max(m,seq(`if`(type(i[1],integer),i[1]+m+1,NULL), i = roots(
		 mul(evala(Expand(eval(lcoeff(i,tau),x=x+offset) )),i=L),x)));
		B := [seq(BasSol(i,x,tau, m+offset,N), i=L)];
		j := convert(B,`+`) - [seq(NthTerm(m+i,A), i=0..2*N)];
		forget(NthTerm);
		iv := indets(B,name) minus indets(L,name);
		B := eval(B, SolveTools:-Linear(j, iv));
		iv := iv intersect indets(B);
		B := subs({seq(iv[i] = _Z[i], i=1..nops(iv))},B);
		if _npassed > 4 and type(V,name) then
		    V :=  proc(c,n) options remember;
			local L, x, tau, ord, offset;
			L, x, tau, ord, offset := procname(c,infinity);
			if not type(n - offset, posint) then
				'procname(_passed)'
			else
				L := eval(L, x=n-ord);
				evala(-add(coeff(L,tau,i)*procname(c,n+i-ord),
					i = 0..ord-1) / lcoeff(L,tau))
			fi
		    end:
		    for i to nops(L) do
			V(i,infinity) := L[i], x, tau, degree(L[i],tau), m+offset;
			for j from 0 to 2*N do
				V(i,j+m+offset) := B[i][j+1]
			od
		    od
		fi;
		return un = add((u || i)(n), i=1..nops(L)),
		  seq(toRecurrence(SimplifyR(L[i],x,tau), (u||i)(n), x, tau), i=1..nops(B)),
		  "Relation holds for", n >= offset + m
	else
	  if r = L then
		Range := "Relation holds for", n >= offset, "Initial terms", init
	  else
		while m > 0 and Normalizer(add(eval(coeff(r,tau,i),x = offset+m-1) *
			NthTerm(i+m-1,A),i = 0..N)) = 0 do	m := m-1
		od;
		lc := {seq(m+i, i=0..N-1),seq(`if`(type(i[1]-m,nonnegint),i[1]+N, NULL),
			i=roots(evala(Expand(subs(x = x + offset, lcoeff(r,tau)))),x))};
		Range := "Relation holds for", n >= offset + m,
			  "Initial terms", {seq(u(i+offset) = NthTerm(i, A), i=lc)}
	  fi
	fi;
	forget(NthTerm);
	`if`(r=L, R, toRecurrence(SimplifyR(r,x,tau), u(n), x, tau)), Range
end:

# A variation on "BasisSolutions" described in Section RFactorsLCLM below.
BasSol := proc(L,x,tau, offset,n) local A, i, a;
	A := [L,x,tau,degree(L,tau),offset,a];
	for i to degree(L,tau) do
		NthTerm(i-1, A) := a[i]
	od;
	[seq(NthTerm(i, A), i=0..2*n)]
end:


###########################################################################
###      Write L as an LCLM of lower order operators when possible      ###
###########################################################################

# Write L as a Complete LCLM of operators, meaning that the output
# operators can not be further decomposed as an LCLM of lower order
# operators (that does not imply that they are irreducible!) (instead,
# it means that they have a {unique} irreducible left-factor).
#
# This allows a sequence u(n) to be decomposed as u(n) = u1(n) + ...
# Completeness means that each term does not allow a further sum
# decomposition. We need the recurrence for u(n) to have minimal order
# which is why "SumDecompose" sits inside program MinimalRecurrence.
#
RFactorsLCLM := proc(L) local S, ds, d, v, n, red, x,tau, id;
	options remember;
	x, tau := _Env_LRE_x, _Env_LRE_tau;
	n := degree(L,tau);
	if n < 2 then return {L} fi;
	S := RFactorsHeuristic(L);
	if Complete(S,L,x,tau,'red') then return red fi;
	ds := [seq(`if`(n=2*d,d,op([d,n-d])),d=1..floor(n/2))];
	id := indets(L) union {x,tau} union indets(L, {'RootOf','radical','nonreal'});
	for d in ds do
		v := SelectFactors(RFactors(L,d),id) minus S;
		if v <> {} then
			S := S union v;
			if Complete(S,L,x,tau,'red') then return red fi;
		fi
	od;
	{L}
end:

# At any stage in RFactorsLCLM,  S is a subset of the set of right-factors.
# This program checks if L = LCLM(S).
# It does this by taking bases of solutions of the operators in S and
# checking if their combined rank is order(L).
# These BasisSolutions are represented by their values in some range.
#
# If L = LCLM(S), the program checks if still holds after removing elements
# of S, starting with the largest ones. After that, to ensure a complete
# decomposition, it applies recursion (RFactorsLCLM) to S.
#
Complete := proc(S, L,x,tau, red) local R, A, B, Ord, i;
	Ord := degree(L,tau);
	if add(degree(i,tau),i=S) < Ord then return false fi;
	B := BasisSolutions(L,x,tau,0) ;
	B := {seq(SolSubSpace(R, B,x,tau), R=S)};
	if LinearAlgebra:-Rank(Matrix([op(B)])) < Ord then
		return false
	elif _npassed < 5 then
		return true
	fi;
	A := [seq([i, degree(i,tau), degree(i,x)], i=S)];
	A := sort(A, (i,j) -> (i[2]>j[2] or (i[2]=j[2] and i[3]>j[3])));
	A := map(i -> i[1], A);
	i := 1;
	while i <= nops(A) do
		B := subsop(i = NULL, A);
		if Complete(B, L,x,tau) then A := B else i := i+1 fi
	od;
	A := {op(A)}; i := nops(A);
	if _passed[-1] <> 0 then
		A := map(op,map(RFactorsLCLM,A));
		if nops(A) > i then
			procname(A, L,x,tau, 'A', 0)
		fi
	fi;
	red := A; true
end:
BasisSolutions := proc(L,x,tau, offset) local A, i, j, m, n, v, a;
	options remember;		# remember table is cleared when done.
	m := 0;
	n := degree(L,tau);
	v := [tcoeff(L,tau), lcoeff(L,tau)];
	if offset <> 0 then
		v := evala(Expand( subs(x=x+offset, v) ))
	fi;
	for j to 2 do for i in roots(v[j],x) do if type(i[1],integer) then
		m := max(m, i[1]+1)
	fi od od;
	m := m + offset;
	A := [L,x,tau,n,m,a];
	for i to n do
		NthTerm(i-1, A) := a[i]
	od;
	[m, {seq(a[i], i=1..n)}, [seq(NthTerm(i, A), i=0..2*n)]]
end:

# The subspace of SPAN(B) that satisfies R.
#
SolSubSpace := proc(R, B, x,tau) local S, m, n, var, i;
	options remember;		# remember table is cleared when done.
	m, n:= B[1], nops(B[2]);
	S := SolveTools:-Linear({ApplyR(m, B[3], R, x,tau,degree(R,tau))}, B[2]);
	S := eval(B[3][1..n], S);
	var := indets(S) intersect B[2];
	if nops(var) <> degree(R,tau) then error "dimension mismatch" fi;
	seq(map(coeff,S,i), i=var);
end:
ApplyR := proc(m,b,R, x,tau,ord) local i,j;
	seq(add(eval( coeff(R,tau,i), x=m+j) * b[1+i+j], i=0..ord), j=0..nops(b)-ord-1)
end:

# In case there are families with infinitely many factors, see procedure PluckerRelations
# for the format, then we have to select finitely many of them. This works best when
# these families have low dimension d and few equations eq, which tends to be the case
# when when using MinimalRecurrence. 
#
SelectFactors := proc(S,id) local i, R, a, d, eq, s, v;
  if type(S,list) then
	R, d, eq := S[1], S[3], S[5];		# R is d-dimensional family with equations eq.
	v := indets(R) minus id;			# Find a value for the unknowns in v.
	for s from 0 to nops(v)-d do
		a := solve(eq union v[s+1..s+d], v);	# Want 1 member of a d-dimensional family,
		if a<>NULL then				# so we equated d unknowns from v to 0.
			a := sort([a],length)[1];
			if indets(map(rhs,a), {'RootOf','radical','nonreal'}) minus id = {} then
				# Avoiding new algebraic numbers
				eq := {seq(i=0, i=indets(map(rhs,a)) intersect v)};
				return eval(eval(R, a), eq)
			fi
		fi
	od;
	NULL
  elif type(S,set) then
	map(procname,S,id)
  else
	S
  fi
end:

###########################################################################
######                   Generalized Exponents                       ######
###########################################################################

# Mark van Hoeij, 2020.

# Input:
#
#  L must be in C[x, x^(-1), tau, tau^(-1)], use primpart(L,tau) to end up there.
#
#  x, tau are names. Multiplication in C[x,tau] is tau * x = (x+1) * tau
#
#  ext (optional) is a set that describes the field of constants, if not
#      specified then ext = the field generated by the coefficients in L.
#
GenExp := proc(L, x, tau::name, ext::set)  local l,d,i,c,v,RES,s,r,m,P;

    if nargs=1 or type(x,set) then
	return procname(L, _Env_LRE_x, _Env_LRE_tau, _passed[2..-1])
    elif not has(L,tau) then
	error "expecting an operator in", tau
    elif not type(L, ratpoly(anything,x)) then
	error "expecting coeffients to be rational functions in", x
    elif not type(L, polynom(anything,x)) then
	return procname(primpart(L,tau), _passed[2..-1])
    elif _npassed=3 then
	return procname(_passed, indets(L, {'RootOf','radical','nonreal'}))
    fi;

    # Compute degrees of the coefficients and use evala to make sure
    # the degrees are correct. Example 1: degree(1+(1+1/t-(t+1)/t)*x, x);
    #
    l, d := ldegree(L,tau), degree(L,tau);
    for i from l to d do
	c := coeff(L, tau, i);
	v[i] := degree(c, x);
	if c<>0 and evala(coeff(c,x,v[i])) = 0 then
		return procname(collect(L, tau, i -> evalaLcoeff(i,x)), x, tau, ext)
	fi
    od;

    # Now use the degrees to compute the tau-NewtonPolygon + Polynomials.
    # Factor these polynomials over Q(ext) and then compute the rest of
    # the generalized exponents using GenExp_sc:
    #
    RES := NULL; m := d;
    while m > l do
	# Next slope in the tau-polygon:
	s := min(seq((v[m]-v[i])/(m-i), i = l .. m-1));
	# Points on slope s:
	r := {seq(`if`( v[m]-v[i] = s*(m-i), i, NULL), i = l .. m)};
	# First point on slope s:
	m := min(op(r));
	# Newton polynomial:
	P := add(coeff(coeff(L,tau,i),x,v[i])*x^((i-m)/denom(s)), i=r);
	P := select(has, evala(Factors(P, ext))[2], x);
	RES := RES, GenExp_sc(L,x,tau,ext,l,d, s,P,max(seq(i[2],i=P)),VarTR())
    od;
    RES := {seq(ShortForm(op(i),evalb(_passed[-1]='LongForm')), i=[RES])};
    userinfo(3,'LREtools',"Generalized Exponents", RES);
    RES
end:

# Produce variables to represent Generalized Exponents.
#
VarTR := proc() local N,T0,T1,T2,T3,T4,T5,T6,T7,T8; options remember;
	N,T0,T1,T2,T3,T4,T5,T6,T7,T8
end:
VarT := proc(n) local T; options remember;
	if n<9 then VarTR()[n+2] else T[n] fi
end:

# Use evala to make sure that lcoeff(nonzero poly) really is not zero.
# See Example 1 above why this can be necessary.
evalaLcoeff:=proc(a,x) local d,i;
	d := degree(a,x);
	if a = 0 or evala(coeff(a,x,d)) <> 0 then return a fi;
	procname(add(coeff(a,x,i)*x^i, i=ldegree(a,x)..d-1),x)
end:

# If K is a difference field and r is in K then tau -> r*tau is an
# endomorphism of K[tau]. The image of L is denoted as Symprod(L, r*tau-1).
# The reason for this name is because the solutions of that operator
# are products of solutions of L and the solution of r*tau-1.
#
# How to compute a Generalized Exponent?  The idea is: Repeatedly make use
# of SymProd until the operator has a generalized exponent e = 1, which
# is equivalent to the indicial equation having 0 as a root.
# For that to happen, the first necessary condition is that the tau-polygon has
# a slope equal to 0, and its polynomial has a root 1.
# This first necessary condition is equivalent to saying that the Delta polygon
# has a slope > 0. As long as that slope is < 1 we can make another SymProd to
# make that slope higher, until we have a slope >= 1. That corresponds to a
# non-constant indicial equation, whose roots then give whats left of e.

# Send L = sum(a_k tau^k) to Ls = SymProd(L, T1*tau - 1) = sum(a_k (T1 tau)^k).
# This sends slope s to 0. Then compute c = RootOf(poly for slope s) and
# compute SymProd(Ls, c*tau-1) which sends c to 1. Then call GenExpDelta.
#
GenExp_sc := proc(L,x,tau,ext,l,d, s,P,Mult, N,T0,T1)
    local m, mm, k, a_k, dg, v, c, ORD, i, Ls, RES, ds, ns, j, e;

    m := min(seq(s*k-degree(coeff(L,tau,k),x), k=l..d)); # min(val(a_k T1^k))
    mm := `if`(s=0, 0, mods(m/s, denom(s))); # Makes v-(k-mm)*s below in Z.
    for k from l to d do
	a_k := coeff(L,tau,k);	   # L = sum(a_k tau^k, k=l..d)
	dg := degree(a_k, x);	   # degree(a_k,x) = -val(a_k)
	v := s*k-dg-m;		   # v = val(a_k T^s) - min(val). Hence v >= 0.
	if v > Mult then	   # v is used to trim a_k (v > Mult means trim(a_k) = 0)
	    c[k] := NULL
	else
	    ORD := Mult-ceil(v)+1; # Use only terms from a_k T1^k/T0^v with val <= Mult
	    # This represents the Symmetric Product sum(a_k (T1 tau)^k)) applied to x^N:
	    c[k] := [convert(series(
		  # Rewrite a_k in terms of T0 = 1/x and scale it to T0-valuation 0:
		series(add(coeff(a_k,x,i)*T0^(dg-i), i = dg-Mult+ceil(v) .. dg), T0=0,ORD)*
		  # (T1 * tau)^k is the following series times T1^k * tau^k :
		series(`if`(k<0,mul(1+i*T0, i=k..-1),1/mul(1+i*T0, i=1..k-1))^s, T0=0,ORD)*
		  # tau^k(x^N)/x^N where N is a variable is:
		series((1+k*T0)^N, T0=0, ORD),    T0=0,ORD),polynom),
		  # The factor T0^v = T0^(v-(k-const)*s) * T1^(k-const) belongs to tau^k:
		T0^(v-(k-mm)*s) * T1^(k-mm), v, k]
	fi
    od;
    # Symmetric product of L that sent the slope s to 0 and scaled m = min(val) to 0:
    Ls := [seq(c[k], k=l..d)];  RES := NULL;
    for i in P do
	c := RootOf(i[1],x);          # Goal: SymProd(Ls, c*tau - 1) to send c to 1
	ds, ns := denom(s), numer(s); #       Insert c so that T1^ds = c * T0^ns
	e := [[T0, T0    = 1/x      , 1      ,   1    ,   1           ],
	      [T1, T1^ds = c * T0^ns, s      ,   ds   , degree(i[1],x)]];
	#     [Tn, equation(Tn)     , val(Tn), Ram(Tn), AlgDeg(Tm)    ]
	v := [seq(`if`(j[3]>i[2],NULL, [
		TrimT0N(j[1], i[2], i[2]-j[3], N,T0), # In C[T0,N]. Valuation zero.
		EvalMonomial(j[2],e), j[3], # Monomial and its valuation j[3] >= 0
		j[4] # power of tau         # Call to EvalMonomial is optional here.
	]), j=Ls)];
	RES := RES, GenExpDelta(ext,v,x,e,N, [i[2],i[2]], 1)
    od;
    RES
end:
# f is C[N, T0]. Remove terms with N-degree > Ncut or valuation > T0cut
#
TrimT0N := proc(f, Ncut, T0cut, N,T0) local i;
    add(TrimN(coeff(f, T0, i), Ncut, N) * T0^i, i=0..floor(T0cut))
end:
TrimN := proc(f, Ncut, N) local i;
    add(coeff(f, N, i) * N^i, i=0..Ncut)
end:
EvalMonomial := proc(M, e) local p,q,i,T,k,d,m;
    # e = [ [T0,    T0 = 1/x     ,   1    ,  1,       1   ],
    #       [T1, T1^r1 = c1*T0^n1, val(T1), r1, algDeg(c1)],  ...]
    p := M; q := 1;
    for i from nops(e) to 2 by -1 do T := e[i,1]; if has(p, T) then
	k := degree(p, T); d := e[i,4]; m := modp(k,d); # use mods?
	q := q * T^m;
	p := coeff(p, T, k) * EvalMonRem(T, k-m, e[1..i], d)
    fi od;
    p * q
end:
# Compute T^k where k is a multiple of d and T in {T1, T2, ...}
#
EvalMonRem := proc(T,k,e,d) options remember;
	if k = 0 then return 1 fi;
	evala(EvalMonomial(procname(T,k-sign(k)*d,e,d),e) * rhs(e[-1,2])^sign(k))
end:

# The Delta Polygon is the Newton Polynomial of L written as an element of
#	C(( T0^(1/r) ))[Delta]
# where Delta = tau - 1.  We do not explicitly write L this way, instead,
# we apply L to x^N where N is a variable. That produces an element Npoly in
#	C(( T0^(1/r) ))[[ N ]] 
# and we compute the Newton Polygon of that. To ensure that we only compute
# slopes of interest, we cut Npoly down to N-degree Top[1] (we throw away
# everything of degree > Top[1]). We also throw away terms with valuation
# greater than Top[2] since those can not affect the Newton Polygon/Polynomials.
#
GenExpDelta := proc(ext,Ls,x, e, N, Top,SlopeBound,Npoly)
	local s,P,   k,cf,v,r,m,R,oldR,   i,d,T,Tn;
    if _npassed=7 then return procname(_passed,
	collect(add(i[1]*i[2], i=Ls),{N,seq(i[1],i=e)},'distributed'))
    fi:
    for k from Top[1] to 0 by -1 do
	# Evaluate only terms needed for the maximal slope.
	# This way in recursion, no terms will be recomputed.     max(val):
	cf[k], v[k] := LowestValuationTerm(coeff(Npoly, N, k), e, Top[2]-(Top[1]-k)*s,
		Top[2]-(Top[1]-k)*SlopeBound -`if`(k=Top[1],1/1000,0) ); # < min(val)
	s := `if`(k=Top[1], 0, max(s, (v[Top[1]]-v[k])/(Top[1]-k))) # max slope found so far
    od:	                                                            # 0 <= s < SlopeBound <= 1
	r := {seq(`if`(v[Top[1]]-v[k] = s*(Top[1]-k), k, NULL), k=0..Top[1])};
	if s=0 then m := 0 else m := min(op(r)) fi;
	oldR := mul(k[4],k=e); R := denom(s*oldR); # old and new ramification

	# Newton Polynomial for the maximal slope:
	P := add(cf[k]*x^((k-m)/R), k=r);
	d := degree(P,x);
	T := TermWithVal(R*s, e);
	P := add(EvalMonomial(coeff(P,x,i)*numer(T)^(d-i)*denom(T)^i,e)*x^i, i=0..d);
	P := select(has,evala(Factors(evala(Primpart(P,x)),ext union indets(e,RootOf)))[2],x);
	Tn := VarT(nops(e));
	if s = 0 then
	    P := IntDiff([seq(collect(i[1]/lcoeff(i[1],x),x,evala)$i[2],i=P)],x,oldR);
	    return seq([e,i,x,N], i=P)
	fi;                                                      # val, ram, algdeg
	seq(SymprodTn([op(e),[Tn,Tn^R = RootOf(i[1],x)*e[1,1]^R/T, 1-s,  R , degree(i[1],x)]],
	    [i[2], Top[2]-s*(Top[1]-i[2])], s, ext,Ls,x,N), i=P),
	# Treated slope s with SymprodTn. Use recursion for slopes < s
	`if`(m=0, NULL, procname(_passed[1..5], [m, v[m]], s, Npoly))
end:
ShortForm := proc(e,f,x,N,LongForm) local i;	# Convert to a shorter format.
	if LongForm then
		[op(e),[op(f), 1, 1, degree(f[1],x)]]
	else
		[e[2,1]*mul(1+e[i,1],i=3..nops(e))*(1+N/x),
		 seq(eval(e[i,2],e[1,2]),i=2..nops(e)), N = [seq(RootOf(f[1],x)+i, i=f[2])]]
	fi
end:
ValTerm := proc(t, e) local i;  # Each entry of e:  [Tn, equation, val(Tn), ram, algdeg]
	add(degree(t,i[1])*i[3], i=e)
end:
TermWithVal := proc(v,e) # Product of T0,T1,... with valuation v
	local w, i, q, f, m;
	w := v; q := 1;
	for i from nops(e) to 1 by -1 do f := e[i]; if f[4] > 1 then
		m := mods(w/f[3], f[4]);
		q := q * f[1]^m;  w := w - m*f[3]
	fi od;
	q * f[1]^w
end:
SymprodTn := proc(e, Top,s, ext,Ls,x,N)
	local L,i,j,A;
	L := NULL;     # Compute a symmetric product:
	for j in Ls do if j[3] <= Top[2] then
		A := Tntauk(j[4], e[1,1], e[-1,1], e[-1,3], Top[2]-j[3]);
		A := collect(A*j[2],{seq(i[1],i=e)},'distributed');
		A := [`if`(type(A,`+`),op(A),A)];
		A := add(`if`(ValTerm(i,e) <= Top[2], EvalMonomial(i,e), 0), i=A);
		L := L,[TrimT0N(j[1], Top[1], Top[2]-j[3], N, e[1,1]), A, j[3], j[4]]
	fi od;
	GenExpDelta(ext,[L],x, e, N, Top, s)
end:
# Compute ((1+Tn)*tau)^k = output * tau^k. Compute terms up to valuation V.
# 0 < val(Tn) = v < 1
#
Tntauk := proc(k, T0, Tn, v, V)
	local p,s,i;
	if   V <   v or k=0 then 1
	elif V < 2*v or k=1 then 1 + k*Tn
	else
		p := Tn^floor(1/v);
		s := `if`(k<0,1/mul(1+Tn/(1+i*T0*p)^v, i=k..-1),
				mul(1+Tn/(1+i*T0*p)^v, i=0..k-1));
		s := convert(series(s, Tn=0, floor(V/v)+1), polynom);
		collect(eval(s, T0=T0/p),{T0,Tn},'distributed')
	fi
end:
# Output: only the term of lowest valuation, if its
# valuation is in (BelowMinVal, MaxVal].
#
LowestValuationTerm := proc(F, e, MaxVal, BelowMinVal)
	local B,v,m,M,i,w;
	B := BelowMinVal;
	if F = 0 then return 0, infinity fi:
    do
	v := MaxVal; m := 0; M := 0;
	for i in [`if`(type(F,`+`),op(F),F)] do
		w := ValTerm(i, e);
		if  w <= B then M := M + i
		elif w = v then m := m + i
		elif w < v then m :=     i; v := w
		fi;
	od; # Uncomment this for bugtesting:
	    # if evala(M) <> 0 then error "bug detected" fi;
	# Either apply EvalMonomial to terms of m or to 2nd argument for GenExpDelta
	m := evala(m);
	if m <> 0     then return m, v        fi;
	if v = MaxVal then return 0, infinity fi;
	B := v;	    # lprint("Terms cancelling out")
    od:
end:
# v is a list of irreducible polynomials. Compute which of them
# are integer/r shifts of each other, determine the "minimal" ones, and
# the non-negative integers you have to shift them with to get the other ones.
#
IntDiff := proc(v::list, x::name, r::posint)
    local i,j,d,e,large,dif,n;
    large := {}; n := nops(v); for i to n do dif[i] := 0; d[i] := degree(v[i],x) od;
    for i from 1 to n-1 do for j from i+1 to n do if d[i]=d[j] then
	e := evala(coeff(v[i],x,d[i]-1)-coeff(v[j],x,d[i]-1));
	if type(r*e/d[i],integer) and evala(v[i]-subs(x=x+e/d[i],v[j]))=0 then
	    if e >= 0 then
		dif[i] := dif[i], e; if e > 0 then large := large union {j} fi
	    fi;
	    if e <= 0 then
		dif[j] := dif[j], -e; if e < 0 then large := large union {i} fi
	    fi
	fi
    fi od od;
	[seq(`if`(member(i,large),NULL,[v[i],sort([dif[i]])]),i=1..n)]
end:


###############################################################################
###     To get a degree bound we have to multiply generalized exponents     ###
###     (taking a Norm when defined over an extension), and an AdjustTerm   ###
###     which is the sum of the valuations of the differences.              ###
###############################################################################

eNorm := proc(e, x, tau, ext) local ext2, p, S, i, s, R, r, a;
	ext2 := ext union indets(e[2], RootOf);
	# p := evala(Norm(tau-rhs(e[2,2]), ext2, ext))^mul(i[4]*i[5], i=e[3..-1])
	# has NewtonPolygon terms, but this p has more terms at non-integer slopes:
	p := eTopCoeffs(e, x, tau, ext);
	S := 0;
	for i in e[3..-1] do 
		ext2 := ext2 union indets(i,RootOf);
		S := S * i[5];
		if type(i[2],list) then  # This is the last entry of e
			s := -coeff(i[1],x,i[5]-1)/coeff(i[1],x,i[5])
		elif numer(i[3])=1 then
			s := EvalMonomial( -(-i[1])^denom(i[3]) / e[1,1], e) / denom(i[3])
		else
			next
		fi;
		S := S + evala(Trace(s, ext2, ext));
	od;
	r := mul(i[4],i=e);     # total ramification in e
	R := mul(i[5],i=e) * r; # total degree [C((t))[e]:C((t))]
	a := ( AdjustmentTerm(e,e) - R*(e[2,3]+1) )/2;
	S := collect(coeff(p,tau,0) * (1 + (r * S - a)/ x), x, evala);
	(-1)^R * S, R, R * e[-1, 2], e[2, 3], p
end:
AdjustmentTerm := proc(e1, e2) local e, N, S, d, i, m, mn, n, v1, v2, vm;
	v1 := e1[2, 3];
	v2 := e2[2, 3];
	vm := min(v1, v2);
	d := 2;
	while d < nops(e1) and d < nops(e2) and e1[d] = e2[d] do
		d := d+1
	od;
	# Now d is the first entry where e1 and e2 disagree.

	S := vm; N := 1; m := 0;
	for i from 2 to d-1 do
		n := e1[i,4] * e1[i,5];
		N := n * N;
		mn := min( e1[i+1,3],e2[i+1,3]);
		S := n^2 * S + N * (mn-m);
		m := mn
	od;
	S * mul(mul(i[4]*i[5], i=e[d..-1]), e=[e1,e2])
end:

# In addition to computing the "Norm + AdjustmentTerm", which gives the dominant
# part of the determinant, we also compute the top coefficients. That includes
# the Newton polygon, Newton polynomial, as well as (in case of non-integer slope)
# the coefficient of the highest integer point in the polygon.
#
eTopCoeffs := proc(e, x, tau, ext)
	local S, dS, ext2, i, id, t;
	S := tau - expand(mul(1+e[i,1], i=3..nops(e)-1));
	id := {tau, seq(e[i,1],i=1..nops(e)-1)};
	S := add(`if`(ValTerm(t,e)<1, EvalMonomial(t,e), 0), t=S);
	for i from nops(e)-1 to 2 by -1 do if i=2 or has(S, e[i,1]) then
		if i=2 then S := subs(tau=tau/e[2,1], S) fi;
		if lhs(e[i,2])=e[i,1] then
			S := subs(e[i,2],S)
		else
			S := resultant(numer(S), numer(lhs(e[i,2])-rhs(e[i,2])), e[i,1])
		fi;
		S := collect(S/lcoeff(S,tau), id, 'distributed');
		if not type(S,`+`) then error "unexpected format" fi;
		dS := degree(S,tau);
		S := add(`if`( ValTerm(t,e)+`if`(i=2,e[2,3]*(degree(t,tau)-dS),0)<1,
			EvalMonomial(t,e), 0), t=S)
	fi od;
	S := eval(S, e[1,2]);
	ext2 := ext union indets(e,RootOf);
	if ext2 <> ext or e[-1,5]>1 then
		S := collect(evala(Norm(S, ext2, ext))^e[-1,5], {x,tau}, 'distributed');
		dS := degree(S,tau);
		S := add(`if`(-degree(t,x)+e[2,3]*(degree(t,tau)-dS)<1, t, 0), t=S)
	fi;
	collect(S,tau,Normalizer)
end:

# Faster version, does the same.
#
eTopCoeffs := proc(e, x, tau, ext)
	local S, dS, ext2, i,j, id, t, z,pz,S0;
	S := tau - mul(1+EvalMonomial(e[i,1],e), i=3..nops(e)-1);
	id := {tau, seq(e[i,1],i=1..nops(e)-1)};
	for i from nops(e)-1 to 2 by -1 do
		if i > 2 and not has(S, e[i,1]) then next fi;
		dS := degree(S,tau);
		if i=2 then S := e[2,1]^dS * subs(tau=tau/e[2,1], S) fi;
		S := collect(S, id, 'distributed', Normalizer);
		S := add(`if`(ValTerm(t,e)+`if`(i=2,e[2,3]*(degree(t,tau)-dS),0)<1,
			EvalMonomial(t,e), 0), t=S);
		S0 := S;
		pz := NumberTheory['CyclotomicPolynomial'](e[i,4],z);
		for j to e[i,4]-1 do
			S := rem(S * rem(subs(e[i,1] = z^j * e[i,1], S0),pz,z), pz,z);
			S := collect(S, id, 'distributed', Normalizer);
			dS := degree(S,tau);
			S := add(`if`(ValTerm(t,e)+`if`(i<3,e[2,3]*(degree(t,tau)-dS),0)<1,
				EvalMonomial(t,e), 0), t=S);
		od;
		S := collect(S, id, 'distributed', Normalizer);
		if has(S, {z, e[i,1]}) then error "unexpected" fi
	od;
	S := eval(S, e[1,2]);
	ext2 := ext union indets(e,RootOf);
	if ext2 <> ext or e[-1,5]>1 then
		S := collect(evala(Norm(S, ext2, ext))^e[-1,5], {x,tau}, 'distributed');
		dS := degree(S,tau);
		S := add(`if`(-degree(t,x)+e[2,3]*(degree(t,tau)-dS)<1, t, 0), t=S)
	fi;
	collect(S,tau,Normalizer)
end:



###############################################################################
###     Now form all possible products, which gives the dominant part of    ###
###     the determinant and top coefficients for all possible RFactors.     ###
###############################################################################

DetFactors := proc(L, x, tau, ext) local i,j,k,S;
	S := SubLists(GenExp(_passed,'LongForm'), x, tau, ext);
	eval({seq([i[1], i[2], {seq( [k, mul( `if`(j[2] = k, j[1]^j[3], 1), j=i[3])],
		k = map(l -> l[2], i[3]))}], i = S)}, VarT(1)=tau)
end:
SubLists := proc(E::set, x, tau, ext)
	local j, k, ANS1, ANS2, L, P, R, a, d, le, i, p, ve;
	P[0] := 1;
	le := E[1];	# A generalized exponent in LongForm
	d, R, L, ve, p := eNorm(le, x, tau, ext);
	if nops(L)>1 then
		a := AdjustmentTerm(le, le)
	fi;
	for i from 1 to nops(L) do
		# Maybe these can be better written as cnd lists:
		P[i] := MultiplyLocalDeterminants(P[i-1], d * (1+(L[-i] - (i-1)*a)/x), x)
	od;
	# OUTPUT IF nops(E) = 1:   { [P[i], R*i, {[p^i, ve]}, {[i, le]} ] , i > 0 }
	# Otherwise, we have to combine these with MultiplyLocalDeterminant, +, union
	# and now use the last one for further adjustment-terms.
	#
	ANS1 := {seq( [P[i], R*i, {[p,ve,i,R,le]}], i=1..nops(L))};
	if nops(E)>1 then
		ANS2 := procname(E minus {le}, x, tau, ext);
		ANS1 := ANS1 union ANS2 union {seq(seq( [MultiplyLocalDeterminants(
		  P[i], j[1]*(1 - add(i*k[3]*AdjustmentTerm(le, k[5]), k=j[3])/x), x),
		  R*i + j[2], j[3] union {[p,ve,i,R,le]}], i=1..nops(L)), j=ANS2)}
	fi;
	ANS1;
end:
# Local determinants are of the form (...)*x^d + (...)*x^(d-1)
# where lower-degree terms are irrelevant. So to multiply them,
# we can simply discard the lower degree terms.
#
MultiplyLocalDeterminants := proc(d1, d2, x) local p,d;
	p := d1*d2; d := degree(p,x);
	Normalizer(coeff(p,x,d))*x^d + Normalizer(coeff(p,x,d-1))*x^(d-1)
end:

end module:
#savelib(`LREtools/RightFactors`):