[section:roots2 Root Finding Without Derivatives: Bisection, Bracket and TOMS748]

[h4 Synopsis]

``
#include <boost/math/tools/roots.hpp>
``
		
   namespace boost{ namespace math{
   namespace tools{
   // Bisection
   template <class F, class T, class Tol>
   std::pair<T, T> 
      bisect(
         F f, 
         T min, 
         T max, 
         Tol tol, 
         boost::uintmax_t& max_iter);

   template <class F, class T, class Tol>
   std::pair<T, T> 
      bisect(
         F f, 
         T min, 
         T max, 
         Tol tol);

   template <class F, class T, class Tol, class ``__Policy``>
   std::pair<T, T> 
      bisect(
         F f, 
         T min, 
         T max, 
         Tol tol, 
         boost::uintmax_t& max_iter,
         const ``__Policy``&);

		// Bracket and Solve Root
   template <class F, class T, class Tol>
   std::pair<T, T> 
      bracket_and_solve_root(
         F f, 
         const T& guess, 
         const T& factor, 
         bool rising, 
         Tol tol, 
         boost::uintmax_t& max_iter);

   template <class F, class T, class Tol, class ``__Policy``>
   std::pair<T, T> 
      bracket_and_solve_root(
         F f, 
         const T& guess, 
         const T& factor, 
         bool rising, 
         Tol tol, 
         boost::uintmax_t& max_iter,
         const ``__Policy``&);

	// TOMS 748 algorithm
   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``&);

   // Termination conditions:
   template <class T>
   struct eps_tolerance;

   struct equal_floor;
   struct equal_ceil;
   struct equal_nearest_integer;
   
   }}} // namespaces
   
[h4 Description]

These functions solve the root of some function /f(x)/ ['without the 
need for any derivatives of /f(x)/].

The `bracket_and_solve_root` functions use TOMS Algorithm 748 that is asymptotically the most efficient known,
and have been shown to be optimal for a certain classes of smooth functions.
Variants with and without __policies are provided.

Alternatively, there is a simple bisection routine which can be useful
in its own right in some situations, or alternatively for narrowing 
down the range containing the root, prior to calling a more advanced
algorithm.

All the algorithms in this section reduce the diameter of the enclosing
interval with the same asymptotic efficiency with which they locate the 
root.  This is in contrast to the derivative based methods which may /never/
significantly reduce the enclosing interval, even though they rapidly approach
the root.  This is also in contrast to some other derivative-free methods
(for example the methods of [@http://en.wikipedia.org/wiki/Brent%27s_method Brent or Dekker)]
which only reduce the enclosing interval on the final step.
Therefore these methods return a std::pair containing the enclosing interval found,
and accept a function object specifying the termination condition.
Three function objects are provided for ready-made termination conditions:
/eps_tolerance/ causes termination when the relative error in the enclosing
interval is below a certain threshold, while /equal_floor/ and /equal_ceil/ are
useful for certain statistical applications where the result is known to be
an integer.  Other user-defined termination conditions are likely to be used
only rarely, but may be useful in some specific circumstances.

[h6 Bisection]

   template <class F, class T, class Tol>
   std::pair<T, T> 
      bisect(
         F f, 
         T min, 
         T max, 
         Tol tol, 
         boost::uintmax_t& max_iter);

   template <class F, class T, class Tol>
   std::pair<T, T> 
      bisect(
         F f, 
         T min, 
         T max, 
         Tol tol);

   template <class F, class T, class Tol, class ``__Policy``>
   std::pair<T, T> 
      bisect(
         F f, 
         T min, 
         T max, 
         Tol tol, 
         boost::uintmax_t& max_iter,
         const ``__Policy``&);

These functions locate the root using bisection: function arguments are:

[variablelist
[[f]  [A unary functor which is the function whose root is to be found.]]
[[min] [The left bracket of the interval known to contain the root.]]
[[max] [The right bracket of the interval known to contain the root.
        It is a precondition that /min < max/ and /f(min)*f(max) <= 0/,
        the function signals evaluation error if these preconditions are violated.
        The action taken is controlled by the evaluation error policy.
        A best guess may be returned, perhaps significantly wrong.]]
[[tol]  [A binary functor that specifies the termination condition: the function
        will return the current brackets enclosing the root when /tol(min,max)/ becomes true.]]
[[max_iter][The maximum number of invocations of /f(x)/ to make while searching for the root.]]
]

[optional_policy]

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

or

   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/ when the function returns), rather than 
because the termination condition /tol/ was satisfied.

[h6 Bracket and solve]

   template <class F, class T, class Tol>
   std::pair<T, T> 
      bracket_and_solve_root(
         F f, 
         const T& guess, 
         const T& factor, 
         bool rising, 
         Tol tol, 
         boost::uintmax_t& max_iter);

   template <class F, class T, class Tol, class ``__Policy``>
   std::pair<T, T> 
      bracket_and_solve_root(
         F f, 
         const T& guess, 
         const T& factor, 
         bool rising, 
         Tol tol, 
         boost::uintmax_t& max_iter,
         const ``__Policy``&);

This is a convenience function that calls /toms748_solve/ internally
to find the root of /f(x)/.  It's usable only when /f(x)/ is a monotonic
function, and the location of the root is known approximately, and in
particular it is known whether the root is occurs for positive or negative
/x/.  The parameters are:

[variablelist
[[f][A unary functor that is the function whose root is to be solved.
    f(x) must be uniformly increasing or decreasing on /x/.]]
[[guess][An initial approximation to the root]]
[[factor][A scaling factor that is used to bracket the root: the value
         /guess/ is multiplied (or divided as appropriate) by /factor/
         until two values are found that bracket the root.  A value
         such as 2 is a typical choice for /factor/.]]
[[rising][Set to /true/ if /f(x)/ is rising on /x/ and /false/ if /f(x)/
         is falling on /x/.  This value is used along with the result
         of /f(guess)/ to determine if /guess/ is
         above or below the root.]]
[[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.]]
[[max_iter] [The maximum number of function invocations to perform in the search 
            for the root.]]
]

[optional_policy]
         
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

or

   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 value of /max_iter/ when the function returns), rather than 
because the termination condition /tol/ was satisfied.

[h6 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 two 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 functions differ only by whether values for /f(a)/ and
/f(b)/ are already available.  The toms748_solve parameters are:

[variablelist
[[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.]]
[[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.]]
[[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.]]
]

[optional_policy]

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
or
   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.

   template <class T>
   struct eps_tolerance
   {
      eps_tolerance(int bits);
      bool operator()(const T& a, const T& b)const;
   };
   
`eps_tolerance` is the usual termination condition used with these root finding functions.
Its operator() will return true when the relative distance between /a/ and /b/
is less than twice the machine epsilon for T, or 2[super 1-bits], whichever is
the larger.  In other words, you set /bits/ to the number of bits of precision you
want in the result.  The minimal tolerance of twice the machine epsilon of T is
required to ensure that we get back a bracketing interval: since this must clearly
be at least 1 epsilon in size.
   
   struct equal_floor
   {
      equal_floor();
      template <class T> bool operator()(const T& a, const T& b)const;
   };
   
This termination condition is used when you want to find an integer result
that is the /floor/ of the true root.  It will terminate as soon as both ends
of the interval have the same /floor/.
   
   struct equal_ceil
   {
      equal_ceil();
      template <class T> bool operator()(const T& a, const T& b)const;
   };

This termination condition is used when you want to find an integer result
that is the /ceil/ of the true root.  It will terminate as soon as both ends
of the interval have the same /ceil/.

   struct equal_nearest_integer
   {
      equal_nearest_integer();
      template <class T> bool operator()(const T& a, const T& b)const;
   };

This termination condition is used when you want to find an integer result
that is the /closest/ to the true root.  It will terminate as soon as both ends
of the interval round to the same nearest integer.

[h4 Implementation]

The implementation of the bisection algorithm is extremely straightforward
and not detailed here.  __TOMS748 is described in detail in:

['Algorithm 748: Enclosing Zeros of Continuous Functions, 
G. E. Alefeld, F. A. Potra and Yixun Shi, 
ACM Transactions on Mathematica1 Software, Vol. 21. No. 3. September 1995. 
Pages 327-344.]

The implementation here is a faithful translation of this paper into C++.

[endsect] [/section:roots2 Root Finding Without Derivatives]

[/ 
  Copyright 2006, 2010 John Maddock and Paul A. Bristow.
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or copy at
  http://www.boost.org/LICENSE_1_0.txt).
]