Algorithm TOMS 748: Alefeld, Potra and Shi: Enclosing zeros of continuous functions

template <class F, class T, class Tol>
std::pair<T, T>
   toms748_solve(
      F f,
      const T& a,
      const T& b,
      Tol tol,
      boost::uintmax_t& max_iter);

template <class F, class T, class Tol, class Policy>
std::pair<T, T>
   toms748_solve(
      F f,
      const T& a,
      const T& b,
      Tol tol,
      boost::uintmax_t& max_iter,
      const Policy&);

template <class F, class T, class Tol>
std::pair<T, T>
   toms748_solve(
      F f,
      const T& a,
      const T& b,
      const T& fa,
      const T& fb,
      Tol tol,
      boost::uintmax_t& max_iter);

template <class F, class T, class Tol, class Policy>
std::pair<T, T>
   toms748_solve(
      F f,
      const T& a,
      const T& b,
      const T& fa,
      const T& fb,
      Tol tol,
      boost::uintmax_t& max_iter,
      const Policy&);

These functions implement TOMS Algorithm 748: it uses a mixture of cubic, quadratic and linear (secant) interpolation to locate the root of f(x). The two pairs of functions differ only by whether values for f(a) and f(b) are already available.

Generally speaking it is easier (and often more efficient) to use bracket and solve rather than trying to bracket the root yourself as this function requires.

This function is provided rather than Brent's method as it is known to be more effient in many cases (it is asymptotically the most efficient known, and has been shown to be optimal for a certain classes of smooth functions). It also has the useful property of decreasing the bracket size with each step, unlike Brent's method which only shrinks the enclosing interval in the final step. This makes it particularly useful when you need a result where the ends of the interval round to the same integer: as often happens in statistical applications for example. In this situation the function is able to exit after a much smaller number of iterations than would otherwise be possible.

The TOMS 748 algorithm parameters are:

f: A unary functor that is the function whose root is to be solved. f(x) need not be uniformly increasing or decreasing on x and may have multiple roots. However, the bounds given must bracket a single root.
a: The lower bound for the initial bracket of the root.
b: The upper bound for the initial bracket of the root. It is a precondition that a < b and that a and b bracket the root to find so that f(a) * f(b) < 0.
fa: Optional: the value of f(a).
fb: Optional: the value of f(b).
tol: A binary functor that determines the termination condition for the search for the root. tol is passed the current brackets at each step, when it returns true, then the current brackets are returned as the result. See also predefined termination functors.
max_iter: The maximum number of function invocations to perform in the search for the root. On exit, max_iter is set to actual number of function invocations used.

The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.

toms748_solve returns: a pair of values r that bracket the root so that:

f(r.first) * f(r.second) <= 0

and either

tol(r.first, r.second) == true

max_iter >= m

where m is the initial value of max_iter passed to the function.

In other words, it's up to the caller to verify whether termination occurred as a result of exceeding max_iter function invocations (easily done by checking the updated value of max_iter against its previous value passed as parameter), rather than because the termination condition tol was satisfied.