function [root,iterations] = bisect(f, a, b, epsilon)
  ## Find a root of the function `f` on the closed interval [a,b].
  ##
  ## It is assumed that there are an odd number of roots within [a,b].
  ## The intermediate value theorem is used to determine whether or not
  ## there is a root on a given interval. So, for example, cos(-2) and
  ## cos(2) are both less than zero; the I.V.T. would not conclude that
  ## cos(x) has a root on -2<=x<=2.
  ##
  ## If `f` has more than one root on [a,b], the smallest root will be
  ## returned. This is an implementation detail.
  ##
  ##
  ## INPUTS:
  ##
  ##   * ``f`` - The function whose root we seek.
  ##
  ##   * ``a`` - The "left" endpoint of the interval in which we search.
  ##
  ##   * ``b`` - The "right" endpoint of the interval in which we search.
  ##
  ##   * ``epsilon`` - How close we should be to the actual root before we
  ##     halt the search and return the current approximation.
  ##
  ## OUTPUTS:
  ##
  ##   * ``root`` - The root that we found.
  ##
  ##   * ``iterations`` - The number of bisections that we performed
  ##     during the search.
  ##

  ## Store these so we don't have to recompute them.
  fa = f(a);
  fb = f(b);

  ## We first check for roots at a,b before bothering to bisect the
  ## interval.  
  if (fa == 0)
    root = a;
    iterations = 0;
    return;
  end

  if (fb == 0)
    root = b;
    iterations = 0;
    return;
  end

  ## Bisect the interval.
  c = (a+b)/2;
  iterations = 1;
  
  ## If we're within the prescribed tolerance, we're done.
  if (b-c < epsilon)
    root = c;
    return;
  end

  if (has_root(fa,f(c)))
    [root,subiterations] = bisect(f,a,c,epsilon);
    iterations = iterations + subiterations;
  else
    [root,subiterations] = bisect(f,c,b,epsilon);
    iterations = iterations + subiterations;
  end
end
