with(Groebner):
with(combinat):

gen_ord_lower := proc(k::integer)
	local eqs, nterms, basis, realroots, posroots, i;

	eqs := [
	  -S + Y + Z,
	  -E + Y + Z*x^2 + Y_t*x^4 + Z*x^4,
	  -E_N + Y_N + Z_N*x^2 + Y_Nt*x^4 + Z_N*x^4,
	  -Y + k*P_N,
	  -Y_t + k*P_Nt,
	  -Y_N + (k-1)*P_N,
	  -Y_Nt + (k-1)*P_Nt,
	  -Z + sum('numbcomb(k,i)*P_Nt^i*x^(i-1)', 'i'=2..k),
	  -Z_N + sum('numbcomb(k-1,i)*P_Nt^i*x^(i-1)', 'i'=2..k-1),
	  -P_N + x + E*x + k*x^3 + k*E_N*x^3,
	  -P_Nt + x + E*x + (k-1)*x^3 + (k-1)*E_N*x^3
	];

	nterms := [ P_Nt, P_N, Z_N, Z, Y_Nt, Y_N, Y_t, Y, E_N, E, S ];

	basis := gbasis(eqs, lexdeg(nterms[1..-2], [S]));

	print("gen ord lower basis computed", k);
	realroots := fsolve(discrim(basis[1], S) = 0, x);

	posroots := select(type, [realroots], positive);
	print(posroots);

	1/min(seq(posroots[i], i=1..nops(posroots)));
end:

gen_rpn_lower := proc(k::integer)
	local eqs, nterms, basis, realroots, posroots, i;

	eqs := [
	  -S + Y + Z,
	  -E + Y + Z + Y_t*x^2 + Z*x^2,
	  -E_N + Y_N + Z_N + Y_Nt*x^2 + Z_N*x^2,
	  -Y + k*P_N,
	  -Y_t + k*P_Nt,
	  -Y_N + (k-1)*P_N,
	  -Y_Nt + (k-1)*P_Nt,
	  -Z + sum('numbcomb(k,i)*P_Nt^i*x^(i-1)', 'i'=2..k),
	  -Z_N + sum('numbcomb(k-1,i)*P_Nt^i*x^(i-1)', 'i'=2..k-1),
	  -P_N + x + E*x^2 + k*x^4 + k*E_N*x^5,
	  -P_Nt + x + E*x^2 + (k-1)*x^4 + (k-1)*E_N*x^5
	];

	nterms := [ P_Nt, P_N, Z_N, Z, Y_Nt, Y_N, Y_t, Y, E_N, E, S ];

        basis := gbasis(eqs, lexdeg(nterms[1..-2], [S]));
	print("gen rpn lower basis computed", k);

        realroots := fsolve(discrim(basis[1], S) = 0, x);

        posroots := select(type, [realroots], positive);
	print(posroots);

        1/min(seq(posroots[i], i=1..nops(posroots)));
end:

gen_alph_lower := proc(k::integer)
        local eqs, nterms, basis, realroots, posroots, i;

	eqs := [
	  -S + Y + Z,
	  -E + Y + Z + Y_t + Z,
	  -E_N + Y_N + Z_N + Y_Nt + Z_N,
	  -Y + k*P_N,
	  -Y_t + k*P_Nt,
	  -Y_N + (k-1)*P_N,
	  -Y_Nt + (k-1)*P_Nt,
	  -Z + sum('numbcomb(k,i)*P_Nt^i', 'i'=2..k),
	  -Z_N + sum('numbcomb(k-1,i)*P_Nt^i', 'i'=2..k-1),
	  -P_N + x + E*x + k*x^2 + k*E_N*x^2,
	  -P_Nt + x + E*x + (k-1)*x^2 + (k-1)*E_N*x^2
	];

	nterms := [ P_Nt, P_N, Z_N, Z, Y_Nt, Y_N, Y_t, Y, E_N, E, S ];

        basis := gbasis(eqs, lexdeg(nterms[1..-2], [S]));
	print("gen alph lower basis computed", k);

        realroots := fsolve(discrim(basis[1], S) = 0, x);

        posroots := select(type, [realroots], positive);
	print(posroots);

        1/min(seq(posroots[i], i=1..nops(posroots)));
end: