with(Groebner):

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

	eqs := [
	  -S + Y + Z, # + F0 - G0 + F1 - G1,
	  -E + Y + Z*x^2 + S*x^4,
	  -Y + k*P,
	  -Z + sum('combinat[numbcomb](k, i)*P^i*x^(i-1)', 'i'=2..k),

	  -P + x + E*x
	];

	nterms := [ P, Z, Y, E, S ];

	basis := gbasis(eqs, lexdeg(nterms[1..-2], [S]));
	print("sf 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:

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

        eqs := [
          -S + Y + Z,
          -E + Y + Z + S*x^2,
          -Y + k*P,
          -Z + sum('combinat[numbcomb](k, i)*P^i*x^(i-1)', 'i'=2..k),

          -P + x + E*x^2
        ];

	nterms := [ P, Z, Y, E, S ];

	basis := gbasis(eqs, lexdeg(nterms[1..-2], [S]));
	print("sf 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:

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

        eqs := [
          -S + Y + Z,
          -E + Y + Z + S,
          -Y + k*P,
          -Z + sum('combinat[numbcomb](k, i)*P^i', 'i'=2..k),

          -P + x + E*x
        ];

	nterms := [ P, Z, Y, E, S ];
										
	basis := gbasis(eqs, lexdeg(nterms[1..-2], [S]));
	print("sf 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: