# 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) have valuations >= 0 at x = i+R with i <= 0
# and b(right-sols) have valuations >= 0 at x = i+R with i > 0
# 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).


# To figure out if b = 1 + DL in D/DL is integral at x = i+R, we need the
# valuations of left-sols (if i <= 0) or right-sols (if i > 0) at x = i+R.
#
# Notation:
#	L is in Q(x)[tau] where tau is the shift operator
#	R = RootOf(a normalized irreducible polynomial)
#	gmin, gmax, a0, an, left, right = data coming from `LREtools/g_p`(L, R)
#
# Input: integer i and data coming from `LREtools/g_p`(L, R)
# Output of the next two procedures: valuation(left-sols) and valuation(right-sols) at x = i+R.
LeftValuation := proc(i, a0, left, gmin, n) local c;
	# left-sols  have valuation 0 at left-n  .. left-1  (the first n 0's in a0)
	c := n+1+i-left;
	if c <= 0 then 0 elif c > nops(a0) then gmin else a0[c] fi;
end:
RightValuation := proc(i, an, right, gmax, n) local d;
	# right-sols have valuation 0 at right+1 .. right+n (the last n 0's in an)
	d := i-right-n+nops(an);
	if  d <= 0 then -gmax elif d > nops(an) then 0 else an[d] fi
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, a0,an,left,right, gmin, gmax;
	P := 1;
	for p in S do
		g := gp[p];
		a0,an,left,right := g[-4..-1];
		if n=1 then an := [seq(i-an[-1], i=an)] fi; # (fixes a quirk in the implementation)
		gmin := min(g[1]); gmax := max(g[1]);
		P := P / mul(collect(subs(x=x-i, p),x)^`if`(i>0,
			RightValuation(k+i, an, right, gmax, n),
			 LeftValuation(k+i, a0, left,  gmin, n)),
		i = min(1,left-k) .. max(0,right-k))
	od;        #    ^ cut-off-points ^   (need to match the "i>0" cut-off).
	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));
		if has(i,x) and degree(R,tau)=2 then
			R := normal(lcoeff(R,tau)/subs(x=x+1,i))*tau^2 + coeff(R,tau,1)*tau + coeff(R,tau,0)*i
		fi;
		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: