function x = EigsOfTridiag(a,b,p,q,tol)
% a is a column n-vector, b is a column (n-1)-vector with nonzero
% entries. Let T be the symmetric tridiagonal with diagonal a and
% subdiagonal b. Denote its eigenvalues by lambda(1)>...>lambda(n).
% tol is a small tolerance.
% Assume p and q satisfy n<=q<=p<=1.
% x = [lambda(q);lambda(q-1);...,lambda(p)];

% Get a bracketing interval [L,R] that includes all the eigenvalues...
OneNorm = max(abs(a)+abs([b; 0])+ abs([0;b])); L = -OneNorm; R =  OneNorm; 
n = length(a);

x = FindEigs(a,b,p,q,tol,L,R,n,0);


function x = FindEigs(a,b,p,q,tol,L,R,iL,iR)
% a is a column n-vector, b is a column (n-1)-vector with nonzero
% entries. Let T be the symmetric tridiagonal with diagonal a and
% subdiagonal b. Denote its eigenvalues by lambda(1)>...>lambda(n).
% tol is a small tolerance.
% Assume p and q satisfy n<=q<=p<=1.
% x = [lambda(q);lambda(q-1);...,lambda(p)];
% iL an integer so  lambda(iL) < L < lambda(iL-1)
% iR an integer so  lambda(iR) < R < lambda(iR-1)

options = optimset('TolX',tol);
if iR == iL-1
    % [L,R] is bracketing and it contains one eigenvalue of interest...
    x = fzero(@(x) CharPoly(x,a,b),[L R],options);
   
else
   mid = (L+R)/2;
   i = InBetween(mid,a,b);
   % lambda(i+1) <= mid < lambda(i)
   if i>=q
       % All the eigenvalues of interest are in [mid,R]
       L = mid; iL = i;
       x = FindEigs(a,b,p,q,tol,L,R,iL,iR);
   elseif i+1<=p
       % All the eigenvalues of interest are in [L,mid]
       R = mid; iR = i;
       x = FindEigs(a,b,p,q,tol,L,R,iL,iR);
   else
    % L <= lambda(q) <...<lambda(i+1) < mid < lambda(i)<...< lambda(p) < R
    % Some eigenvalues are in [L,mid] and some are in [mid,R]...
    x = [FindEigs(a,b,i+1,q,tol,L,mid,iL,i);...
         FindEigs(a,b,p,i,tol,mid,R,i,iR)];
   end
end


   

  

function y = CharPoly(lambda,a,b)
% y = det(T - lambda*I) where T is the n-by-n tridiagonal matrix
% with diagonal a and subdiagonal b.
n = length(a);
b = b.^2;        
pkm1 = 1;
pk = a(1)-lambda;
for k=2:n
    pkm2 = pkm1;
    pkm1 = pk;
    pk = (a(k)-lambda)*pkm1 - b(k-1)*pkm2;
end
y = pk;

function i = InBetween(x,a,b)
% a and b define a tridiagonal matrix T with eigenvalues
%       lambda(n) < lambda(n-1) <...< lambda(1)
% x is a scalar.
% If i = n, then all the eigenvalues of T are to the right of x.
% If i = 0, then all the eigenvalues of T are to the left of x.
% Otherwise, lambda(i+1) < x <= lambda(i)
n = length(a);
b = b.^2;
pkm1 = 1;
pk = a(1)-x;
m = (pkm1*pk<=0);
for k=2:n
    pkm2 = pkm1;
    pkm1 = pk;
    pk = (a(k)-x)*pkm1 - b(k-1)*pkm2;
    m = m + (pk*pkm1<=0);
end
i = n-m;


%  -5 for exponential eval of char poly
% -2 for minor things like returned eigs in reverse order
%  up to -4 for correct but with redundant f-evals