function lambda = bisection(T,i,tol)
%bisection Compute eigenvalue i of the unreduced symmetric tridiagonal
%matrix T using the bisection method.
%
%   lambda = bisection(T,i,tol) computes eigenvalue lambda_i,
%   where lambda_1 < lamabda_2 < ... < lambda_i < ... < lambda_n.

if ~isTriDiag(T)
	error('The matrix is not tridiagonal.');
end
% remove any small entries off the three diagonals
T = diag(diag(T))+ diag(diag(T,1),1)+diag(diag(T,-1),-1);
if ~issym(T)
	error('The matrix is tridiagonal but not symmetric.');
end

n = size(T,1);
if i > n
   error('The matrix only has %d eigenvalues.',n);
end
if nargin == 2
	tol = 1.0e-12;
end

left = -norm(T,'inf');
right = -left;
mid = (left+right)/2;
% continue until finding the eigenvalue or until the enclosing
% interval has length < tol.
while (right-left) >= tol
	% count the number sign changes between successive entries
	% of p0(mid), p1(mid), ..., pn(mid)
   sign_changes = 0;
   havezero = false;
   pkm2 = 1;
   pkm1 = T(1,1)-mid;
   if pkm1 == 0
      pkm1 = -pkm2;
      havezero = true;
   end
   if pkm2*pkm1 < 0
      sign_changes = sign_changes + 1;
   end
   if havezero
      pkm1 = 0;
   end
   for j = 2:n
      havezero = false;
      pk = (T(j,j)-mid)*pkm1 - T(j-1,j)^2*pkm2;
      if pk == 0
         pk = -pkm1;
         havezero = true;
      end
      if pk*pkm1 < 0
         sign_changes = sign_changes + 1;
      end
      if havezero
         pk = 0;
      end
      pkm2 = pkm1;
      pkm1 = pk;
   end      % move to the the appropriate 1/2 length interval
   if sign_changes >= i
      right = mid;
   else
      left = mid;
   end
   mid = (left+right)/2;
end

lambda = mid;
