function [ X, fvals, its, Xhistory ] = bisect( f, X0, tol )
%bisect:  Interval bisection root finder.
% 
% Usage:
%
%   >> [ X, fval, its, Xhistory ] = bisect( f, X0, tol )
% 
% Inputs:
% 
%   f  =  function or function handle for f(x) = 0 function ...
%         must evaluate via >> feval( f, x ) ...
%         ex - ... = bisect( @sin, [ 3 4 ] )
%                  =   ""  ( 'sin', [ 3 4 ] )
%                  =   ""  ( @(x) x^2 - 2, [ 1 2 ] )
%                  =   ""  ( inline( 'x^2 - 2' ), [ 1 2 ] )
%         see Help for "inline" and "function_handle"
%   X0  = [ a b ] = initial interval bracketing root
%   tol = (optional) absolute stopping tolerance: | b - a | < tol
%         (default = 10^{-6})
% 
% Outputs:
% 
%   X        = [ a b ] = final bracketing interval
%   fvals    = (optional) final value of [ f(a) f(b) ]
%   its      = (optional) # iterations
%   Xhistory = [ a1 b1 ; ... ; an bn ] = (optional) array of successive
%                                        bisected intervals

% ~/matlab/courses/CompFinance/bisect.m
% 
% C.G. (9-02)
% C.G. (9-03) fixed bug when reaching limits of machine precision
% C.G. (9-05) minor reformatting

a = X0(1) ;  fa = feval( f, a ) ;       % left end point 'a' and f(a)

if ( fa == 0 )                          % test for exceptional case f(a) = 0
  X = [ a a ] ;  Xhistory = X ;
  fvals = [ 0 0 ] ;
  its = 0 ;
  return
end

b = X0(2) ;  fb = feval( f, b ) ;       % right end point 'b' and f(b)

if ( fb == 0 )                          % test for exceptional case f(b) = 0
  X = [ b b ] ;  Xhistory = X ;
  fvals = [ 0 0 ] ;
  its = 0 ;
  return
end

if ( fa*fb > 0 )
  error( '[a,b] not a bracketing interval: f(a) and f(b) same sign' )
end

if ( nargin < 3 || isempty(tol) )
  tol = 1.e-6 ;                                        % default tolerance
end

X = [ a b ] ;  Xhistory = X ;  its = 0 ;               % initialize ...

while ( ( b - a ) > tol )
    
  m = ( a + b ) / 2 ;                                  % mid-point 'm'

  % test for limits of machine precision

  if ( m == a || m == b )
    X = [ a b ] ;  Xhistory = [ Xhistory ; X ] ;
    fvals = [ fa fb ] ;
    warning( 'no further progress possible' )
    return
  end

  fm = feval( f, m ) ;                                 % f( mid-point 'm' )
    
  % test for exceptional case f(m) = 0

  if ( fm == 0 )
    X = [ m m ] ;  Xhistory = [ Xhistory ; X ] ;
    fvals = [ 0 0 ] ;
    return
  end
    
  % update left or right end point as appropriate

  if ( fa*fm > 0 )
    a = m ;  fa = fm ;
  else
    b = m ;  fb = fm ;
  end
    
  % step completed: add to history, increment counter, etc

  X = [ a b ] ;  Xhistory = [ Xhistory ; X ] ;
  fvals = [ fa fb ] ;
  its = its + 1 ;
    
end