...one of the most highly
regarded and expertly designed C++ library projects in the
world.

— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards

#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

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

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:

- 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. On exist this is set to actual number of invocations performed.

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.

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.

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:

- 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. On exit is set to the actual number of invocations performed.

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.

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.

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:

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

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 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^{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.

The implementation of the bisection algorithm is extremely straightforward and not detailed here. TOMS Algorithm 748: enclosing zeros of continuous functions 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++.