function s=bisection(k,x,a,b)
%input a: main diagonal of symmetric tridiagonal matrix A
%      b: sub(super) diagonal of symmetric tridiagonal matrix A
%      x: the shift in the bisection method.
%      k: degree of the poly. (i.e. the size of the submatrix)
%output s: det(A-xI), as calculated recursively
%
% Note: the number of sign changes in the sturm sequence is equal to the
% number of negative eigenvalues of A.

if k==-1
    s=0;
elseif k==0
    s=1;
else
    if k==1
      s=(a(k)-x)*bisection(k-1,x,a,b);
    else
      s=(a(k)-x)*bisection(k-1,x,a,b)-b(k-1)^2*bisection(k-2,x,a,b);
    end
    out=sprintf('Size of submatrix is %g and its det. is %g',k,s);
    disp(out)
end