read "DeltaHyperSols" ;
read "DeltaSuperForms" ;
_Env_LRE_tau := tau;
_Env_LRE_x   := x;

##########################  Application #1:  ProjectiveHom  ####################################

# Given L1, L2 in C(x)[tau],  compute all G of the form:
#    Solution(operator of order 1) * "an element of C(x)[tau]"
# for which G applied to any solution of L1 gives a solution of L2.
#
ProjectiveHom := proc(L1, L2)
	local m,S,i,j,T;
	global x, tau;
	m := degree(L1,tau);
	S := deltaHS(TensorSystem(adjointOperator(L1), L2), x);
	S := [seq(seq([i[3] * i[1]/i[2], adjustR(j[1..m],L1)], j=i[4]), i=[S])];
	T := [seq(content(i[2],tau), i=S)];
	seq( SolOf(tau-factor(S[i][1] * subs(x=x+1, T[i])/T[i])) * collect(S[i][2]/T[i],tau,factor), i=1..nops(S))
end:
adjointOperator := proc(L) global x,tau; local n, i;
	n := degree(L, tau); 
	add(normal(subs(x = x+i-1, coeff(L, tau, n-i)/lcoeff(L,tau))) * tau^i, i=0..n)
end:
subsDual := proc(L, n, bd) 
	local j;
	global x, tau;
	{seq(bd[j] = -coeff(L, tau, j+1)*bd[0]/coeff(L, tau, 0)+`if`(j = n-1, 0, bd[j+1]), j = 0 .. n-1)} 
end:

# change of basis from v=d[n-1] to d[i]'s
adjustR:=proc(mat,L)
	local action, V, b, i, n, bd, v;
	global x,tau;
	n := degree(L, tau);
	action := subsDual(L, n, bd);
	V[0] := bd[n-1];
	for i to n-1 do 
		V[i] := normal(subs(action, subs(x = x+1, V[i-1]))) 
	od;
	V := SolveTools:-Linear({seq(V[i]-v[i], i = 0 .. n-1)}, {seq(bd[i], i = 0 .. n-1)});
	b := subs(V, add(bd[i]*tau^i, i = 0 .. n-1));
	collect(subs({seq(v[i-1] = mat[i], i=1..n )}, b), tau, normal)
end:




# L1 and L2 are two operators in C(x)[tau].
# This program produces the matrix A for which tau(Y) = AY is
# the system for the tensor product of L1 and L2
#
TensorSystem := proc(L1, L2)
	global x, tau;
	local B,m,n,b,i,j;
	if coeff(L1,tau,0) = 0 or coeff(L2,tau,0) = 0 then
		error "expect tau^0 have non-zero coefficients"
	fi;
	m, n := degree(L1,tau), degree(L2,tau);
	B := [seq(seq([i,j], i=1..m), j=1..n)];
	matrix([seq(TensorInTermsOfBasis(L1,L2, b, B, m,n), b=B)])
end:
TensorInTermsOfBasis := proc(L1,L2, b, B, m,n)
	global x, tau;
	local tau_b, i, j;
	tau_b := b + [1,1];
	if tau_b = [m+1, n+1] then
		add(add(WriteInTermsOf([i,j], B, normal(
			coeff(L1, tau, i-1)/lcoeff(L1, tau) *
			coeff(L2, tau, j-1)/lcoeff(L2, tau) )), i=1..m), j=1..n)
	elif tau_b[1] = m+1 then
		add(WriteInTermsOf([i,tau_b[2]], B, normal(
			-coeff(L1, tau, i-1)/lcoeff(L1, tau) )), i=1..m)
	elif tau_b[2] = n+1 then
		add(WriteInTermsOf([tau_b[1], i], B, normal(
			-coeff(L2, tau, i-1)/lcoeff(L2, tau) )), i=1..n)
	else
		WriteInTermsOf(tau_b, B, 1)
	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:

if run_tests = true then
	L1 := (x-1/2)^2*(x-1)*tau^3 + tau + (x+1/2)^2*(x-3);
	G := (x-5)*tau^2 + x*tau - 5*(x-1);
	L2 := collect(primpart(`LREtools/rightdivision`(`LREtools/LCLM`(L1, G), G)[1],tau),tau,factor);
	ProjectiveHom(L1, L2);
	# OK, matches G.

	L := x*(x+3)*(x+1)*(x^5+7*x^4+15*x^3-32*x-22)*tau^4-(x+2)*(4*x^4+25*x^3+43*x^2+7*x-24)*x*tau^3
	-(x+1)*(x^6+7*x^5+19*x^4+21*x^3-21*x^2-87*x-62)*tau^2+x*(x+2)*(11*x^3+61*x^2+103*x+62)*tau
	+(x-1)*(x+2)*(x+1)*(x^5+12*x^4+53*x^3+97*x^2+46*x-31) ;

	_Env_ratsols_only := true;	# Compute only rational solutions instead of hypergeometric solutions:
	ProjectiveHom(L, subs(tau = -tau,L));
	_Env_ratsols_only := false;
fi;


##########################  Application #2:  Factoring  ####################################


# Input: L in C(x)[tau]
# Output: matrix A such that order-d factors of L correspond to hypergemetric solutions of AY = tau(Y).
ExteriorPowerSystem := proc(L,d)
	global x, tau;
	local n, B, b;
	n := degree(L,tau);
	B := Subsets(n,d);  # if b = [1,3,4] in B then view it as tau^0 wedge tau^2 wedge tau^3
	# Now compute the action of tau on B.
	matrix([seq(WriteInTermsOfBasis(L, b, B, n), b=B)])
end:
# Subsets of {1..n} given as ordered lists
Subsets := proc(n::nonnegint, d::nonnegint) local i;
	options remember;
	if d=0 then
		[ [ ] ]  # 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, b, B, n) global x, tau; local i, tau_b, S, j;
	tau_b := [seq(i+1, i=b)];
	if tau_b[-1] <= n then
		WriteInTermsOf(tau_b, B, 1)
	else
		S := 0;
		tau_b := tau_b[1..-2];
		for i to n do if not member(i, tau_b) then
			S := S + WriteInTermsOf(sort([op(tau_b), i]), B,
				normal(-coeff(L, tau, i-1)/coeff(L, tau, n))
				* mul(`if`(j>i,-1,1), j=tau_b)	)
		fi od;
		map(normal,S)
	fi
end:


BekeBronstein := proc(L, d::posint)
	global x, tau;
	local n,i,j,S;
	n := degree(L,tau);
	if n <= d then
		error "d should be less than the order of the input"
	fi;
	S := deltaHS(ExteriorPowerSystem(L,d), x);
	S := {seq( add(_C[j] * ConvertToOperator(i[-1][j], d), j=1..nops(i[-1])), i = [S])};
	map(SimplifyR, S)

	# After that, [Bronstein 1994] then uses a right-division to find the correct values of _C[1], _C[2], ...
	# However, that right-division produces equations of degree n-d+1, while the Plucker relations have degree 2.
	# So instead of a right-division, I'll implement something else that computes equations of degree 2.
end:
ConvertToOperator := proc(v, d)
	global x,tau;
	local i;
	add( (-1)^i * v[i+1] * tau^(d-i), i=0..d)
end:
SimplifyR := proc(R) global tau;
	sort(collect(primpart(R,tau), tau,factor), tau)
end:


if run_tests = true then
	L6 := `LREtools/LCLM`( (x+3) * (x-1/2) * tau^2 + 3*(x-1/2)*tau - 5 * x^2,  7*(x+2)*(x+1/2) * tau^4 + x*tau^2 - x*(x-5/2));
	BekeBronstein(L6, 2);
	BekeBronstein(L6, 4);
	L4 := `LREtools/LCLM`(tau^2-x, x*tau^2-tau-x^2-2*x-1);
	S := BekeBronstein(L4,2);
fi:


# This is a version of BekeBronstein that, if there are more than two homogeneous unknown _C[i]'s
# then instead of using right-division to compute equations for them (resulting in equations of degree n-d+1)
# we'll use Plucker relations (resulting in equations of degree 2).
#
BBP := proc(L, d::posint)
	global x, tau;
	local n,i,S;
	n := degree(L,tau);
	if n <= d then
		error "d should be less than the order of the input"
	fi;
	S := deltaHS(ExteriorPowerSystem(L,d), x);
	{seq(PluckerRelations(L, n, d, i[-1]), i=[S])}
end:

# Input: L, n, d, V = list-of-poly-sols

PluckerRelations := proc(L, n, d, V)
	global x, tau;
	local R, i, vars, r, P, B, iR, EQ, l, s, Y, sb;
	R := add(_C[i] * ConvertToOperator(V[i], d), i=1..nops(V));
	vars := {seq(_C[i],i=1..nops(V))};
	sb := RowReduceR(R, d, vars); # Change-of-basis to make R smaller.
	R := SimplifyR( subs(sb,R) );
	if R = 0 then
		NULL
	elif d = 1 or n-d = 1 then
		R
	elif nops(V) <= 2 then
		# With 1 or 2 homogeneous variables we only need 1 equation:
		r := `LREtools/rightdivision`(lcoeff(R,tau)^(n-d+1) * L, R)[2];
		if r = 0 then
			R
		elif nops(V)=1 then
			NULL
		else
			r := content(numer(normal(r)), {x,tau});
			seq(SimplifyR(eval(R, i)), i = [solve(r, vars)])
		fi
	else
		# Convert lists to polynomials so we can divide by the gcd:
		P := add(add(subs(sb,_C[i]) * V[i][l] * Y^l, l=1..nops(V[1])), i=1..nops(V));
		P := primpart(P, Y);
		B := Subsets(n,d);
		iR := R;
		EQ := {};
		# Now construct binomial(n-1, d)-1 quadratic equations for the _C[i]
		for i to n-1-d do
			iR := primpart(subs(x=x+1, iR*tau), tau); # = tau^i * R
			EQ := EQ union {seq( add(
			  coeff(iR, tau, l) * CoeffOfP(s, l, B, P, Y),
			  l=i..d+i), s = Subsets(i+d-1, d-2) )}		# s wedge iR = 0
		od:
		EQ := {seq(coeffs(collect(primpart(i,vars),x,primpart),x), i = EQ)} minus {0};
		if EQ = {} then
			R
		else
			l := Groebner:-HilbertDimension(EQ, vars);
			if l = 0 then
				NULL
			elif l = 1 then
				seq(SimplifyR(eval(R, i)), i = [solve(EQ, vars)])
			else
				[R, cat("contains a P^", l-1, "-family. Equations for the _C[i] :"), EQ]
			fi
		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:

# Used to reduce the size.
RowReduceR := proc(R, d, vars)
	global x, tau;
	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], solve({cf = tmp}, v[1]));
		treated := treated union {v[1]}
	od od;
	sb
end:

if run_tests = true then
	L4 := `LREtools/LCLM`(tau^2-x, x*tau^2-tau-x^2-2*x-1);
	BBP(L4,2);
	L6 := `LREtools/LCLM`(tau^2-x+2, tau^2-x, x*tau^2-tau-x^2-2*x-1);
	S := BBP(L6,2);
fi;