with(Groebner):

sf_ord_upper := proc()
	local eqs, nterms, basis, realroots, posroots, i;
	eqs := [
	  -S + A + B + D + x + x,
	  -A + A_B*x^2 + B*A_B*x + D*A_D*x,
	  -A_t + B*A_B*x + D*A_D*x,
	  -A_B + B + B*A_D*x + A_D,
	  -A_D + D + D*A_D*x,
	  -B + A*B_A*x^2 + D*B_D,
	  -B_A + A*x^2 + A*B_A*x^2 + B_D,
	  -B_D + D + D*B_D,
	  -D + x
	];

	nterms := [ D, B_D, B_A, B, A_D, A_B, A_t, A, S ];

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

	print("sf ord upper 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_upper := proc()
        local eqs, nterms, basis, realroots, posroots, i;

	eqs := [
	  -S + A + B + D + x + x,
	  -A + A_B*x^2 + B*A_B*x + D*A_D*x,
	  -A_t + B*A_B*x + D*A_D*x,
	  -A_B + B + B*A_D*x + A_D,
	  -A_D + D + D*A_D*x,
	  -B + A*B_A*x + D*B_D*x,
	  -B_A + A + A*B_A*x + B_D,
	  -B_D + D + D*B_D*x,
	  -D + x
	];

	nterms := [ D, B_D, B_A, B, A_D, A_B, A_t, A, S ];

        basis := gbasis(eqs, lexdeg(nterms[1..-2], [S]));
	print("sf rpn upper 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_upper := proc()
        local eqs, nterms, basis, realroots, posroots, i;
	eqs := [
	  -S + A + B + D,
	  -A + A_B + B*A_B + D*A_D,
	  -A_t + B*A_B + D*A_D,
	  -A_B + B + B*A_D + A_D,
	  -A_D + D + D*A_D,
	  -B + A*B_A + D*B_D,
	  -B_A + A + A*B_A + B_D,
	  -B_D + D + D*B_D,
	  -D + x
	];

	nterms := [ D, B_D, B_A, B, A_D, A_B, A_t, A, S ];

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