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

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

First some `#includes`

that
will be needed.

#include <boost/math/tools/roots.hpp> //using boost::math::policies::policy; //using boost::math::tools::newton_raphson_iterate; //using boost::math::tools::halley_iterate; // //using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits. //using boost::math::tools::bracket_and_solve_root; //using boost::math::tools::toms748_solve; #include <boost/math/special_functions/next.hpp> // For float_distance. #include <tuple> // for std::tuple and std::make_tuple. #include <boost/math/special_functions/cbrt.hpp> // For boost::math::cbrt.

Tip | |
---|---|

For clarity, |

Let's suppose we want to find the root of a number *a*,
and to start, compute the cube root.

So the equation we want to solve is:

f(x) = x³ -a

We will first solve this without using any information about the slope or curvature of the cube root function.

Fortunately, the cube-root function is 'Really Well Behaved' in that it is monotonic and has only one root (we leave negative values 'as an exercise for the student').

We then show how adding what we can know about this function, first just
the slope or 1st derivative *f'(x)*, will speed homing
in on the solution.

Lastly, we show how adding the curvature *f''(x)* too
will speed convergence even more.

First we define a function object (functor):

template <class T> struct cbrt_functor_noderiv { // cube root of x using only function - no derivatives. cbrt_functor_noderiv(T const& to_find_root_of) : a(to_find_root_of) { /* Constructor just stores value a to find root of. */ } T operator()(T const& x) { T fx = x*x*x - a; // Difference (estimate x^3 - a). return fx; } private: T a; // to be 'cube_rooted'. };

Implementing the cube-root function itself is fairly trivial now: the hardest part is finding a good approximation to begin with. In this case we'll just divide the exponent by three. (There are better but more complex guess algorithms used in 'real life'.)

template <class T> T cbrt_noderiv(T x) { // return cube root of x using bracket_and_solve (no derivatives). using namespace std; // Help ADL of std functions. using namespace boost::math::tools; // For bracket_and_solve_root. int exponent; frexp(x, &exponent); // Get exponent of z (ignore mantissa). T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three. T factor = 2; // How big steps to take when searching. const boost::uintmax_t maxit = 20; // Limit to maximum iterations. boost::uintmax_t it = maxit; // Initially our chosen max iterations, but updated with actual. bool is_rising = true; // So if result if guess^3 is too low, then try increasing guess. int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. // Some fraction of digits is used to control how accurate to try to make the result. int get_digits = digits - 3; // We have to have a non-zero interval at each step, so // maximum accuracy is digits - 1. But we also have to // allow for inaccuracy in f(x), otherwise the last few // iterations just thrash around. eps_tolerance<T> tol(get_digits); // Set the tolerance. std::pair<T, T> r = bracket_and_solve_root(cbrt_functor_noderiv<T>(x), guess, factor, is_rising, tol, it); return r.first + (r.second - r.first)/2; // Midway between brackets is our result, if necessary we could // return the result as an interval here. }

Note | |
---|---|

The final parameter specifying a maximum number of iterations is optional.
However, it defaults to So it may be wise to chose some reasonable estimate of how many iterations may be needed, In this case the function is so well behaved that we can chose a low value of 20.
Internally when Boost.Math uses these functions, it sets the maximum iterations
to |

Should we have wished we can show how many iterations were used in `bracket_and_solve_root`

(this information
is lost outside `cbrt_noderiv`

),
for example with:

if (it >= maxit) { std::cout << "Unable to locate solution in " << maxit << " iterations:" " Current best guess is between " << r.first << " and " << r.second << std::endl; } else { std::cout << "Converged after " << it << " (from maximum of " << maxit << " iterations)." << std::endl; }

for output like

Converged after 11 (from maximum of 20 iterations).

This snippet from `main()`

in root_finding_example.cpp
shows how it can be used.

try { double threecubed = 27.; // Value that has an *exactly representable* integer cube root. double threecubedp1 = 28.; // Value whose cube root is *not* exactly representable. std::cout << "cbrt(28) " << boost::math::cbrt(28.) << std::endl; // boost::math:: version of cbrt. std::cout << "std::cbrt(28) " << std::cbrt(28.) << std::endl; // std:: version of cbrt. std::cout <<" cast double " << static_cast<double>(3.0365889718756625194208095785056696355814539772481111) << std::endl; // Cube root using bracketing: double r = cbrt_noderiv(threecubed); std::cout << "cbrt_noderiv(" << threecubed << ") = " << r << std::endl; r = cbrt_noderiv(threecubedp1); std::cout << "cbrt_noderiv(" << threecubedp1 << ") = " << r << std::endl;

cbrt_noderiv(27) = 3 cbrt_noderiv(28) = 3.0365889718756618

The result of `bracket_and_solve_root`

is a pair
of values that could be displayed.

The number of bits separating them can be found using `float_distance(r.first, r.second)`

. The distance is zero (closest representable)
for 3^{3} = 27 but `float_distance(r.first, r.second) = 3`

for cube root of 28 with this function. The result (avoiding overflow) is
midway between these two values.

We now solve the same problem, but using more information about the function, to show how this can speed up finding the best estimate of the root.

For the root function, the 1st differential (the slope of the tangent to a curve at any point) is known.

This algorithm is similar to this nth root algorithm.

If you need some reminders, then derivatives of elementary functions may help.

Using the rule that the derivative of *x ^{n}* for positive
n (actually all nonzero n) is

To see how this extra information is used to find a root, view Newton-Raphson iterations and the animation.

We define a better functor `cbrt_functor_deriv`

that returns both the evaluation of the function to solve, along with its
first derivative:

To '*return*' two values, we use a std::pair
of floating-point values.

template <class T> struct cbrt_functor_deriv { // Functor also returning 1st derivative. cbrt_functor_deriv(T const& to_find_root_of) : a(to_find_root_of) { // Constructor stores value a to find root of, // for example: calling cbrt_functor_deriv<T>(a) to use to get cube root of a. } std::pair<T, T> operator()(T const& x) { // Return both f(x) and f'(x). T fx = x*x*x - a; // Difference (estimate x^3 - value). T dx = 3 * x*x; // 1st derivative = 3x^2. return std::make_pair(fx, dx); // 'return' both fx and dx. } private: T a; // Store value to be 'cube_rooted'. };

Our cube root function is now:

template <class T> T cbrt_deriv(T x) { // return cube root of x using 1st derivative and Newton_Raphson. using namespace boost::math::tools; int exponent; frexp(x, &exponent); // Get exponent of z (ignore mantissa). T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three. T min = ldexp(0.5, exponent/3); // Minimum possible value is half our guess. T max = ldexp(2., exponent/3); // Maximum possible value is twice our guess. const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. int get_digits = static_cast<int>(digits * 0.6); // Accuracy doubles with each step, so stop when we have // just over half the digits correct. const boost::uintmax_t maxit = 20; boost::uintmax_t it = maxit; T result = newton_raphson_iterate(cbrt_functor_deriv<T>(x), guess, min, max, get_digits, it); return result; }

The result of `newton_raphson_iterate`

function
is a single value.

Tip | |
---|---|

There is a compromise between accuracy and speed when choosing the value
of |

Note that it is up to the caller of the function to check the iteration count after the call to see if iteration stoped as a result of running out of iterations rather than meeting the required precision.

Using the test data in /test/test_cbrt.cpp
this found the cube root exact to the last digit in every case, and in no
more than 6 iterations at double precision. However, you will note that a
high precision was used in this example, exactly what was warned against
earlier on in these docs! In this particular case it is possible to compute
*f(x)* exactly and without undue cancellation error, so
a high limit is not too much of an issue.

However, reducing the limit to `std::numeric_limits<T>::digits * 2 / 3`

gave full precision in all but one of the test cases (and that one was out
by just one bit). The maximum number of iterations remained 6, but in most
cases was reduced by one.

Note also that the above code omits a probable optimization by computing z² and reusing it, omits error handling, and does not handle negative values of z correctly. (These are left as the customary exercise for the reader!)

The `boost::math::cbrt`

function also includes these and other
improvements: most importantly it uses a much better initial guess which
reduces the iteration count to just 1 in almost all cases.

Next we define yet another even better functor `cbrt_functor_2deriv`

that returns both the evaluation of the function to solve, along with its
first **and second** derivative:

f''(x) = 6x

using information about both slope and curvature to speed convergence.

To *'return'* three values, we use a `tuple`

of three floating-point values:

template <class T> struct cbrt_functor_2deriv { // Functor returning both 1st and 2nd derivatives. cbrt_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of) { // Constructor stores value a to find root of, for example: // calling cbrt_functor_2deriv<T>(x) to get cube root of x, } std::tuple<T, T, T> operator()(T const& x) { // Return both f(x) and f'(x) and f''(x). T fx = x*x*x - a; // Difference (estimate x^3 - value). T dx = 3 * x*x; // 1st derivative = 3x^2. T d2x = 6 * x; // 2nd derivative = 6x. return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x. } private: T a; // to be 'cube_rooted'. };

Our cube root function is now:

template <class T> T cbrt_2deriv(T x) { // return cube root of x using 1st and 2nd derivatives and Halley. //using namespace std; // Help ADL of std functions. using namespace boost::math::tools; int exponent; frexp(x, &exponent); // Get exponent of z (ignore mantissa). T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three. T min = ldexp(0.5, exponent/3); // Minimum possible value is half our guess. T max = ldexp(2., exponent/3); // Maximum possible value is twice our guess. const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. // digits used to control how accurate to try to make the result. int get_digits = static_cast<int>(digits * 0.4); // Accuracy triples with each step, so stop when just // over one third of the digits are correct. boost::uintmax_t maxit = 20; T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, get_digits, maxit); return result; }

The function `halley_iterate`

also returns a single value, and the number of iterations will reveal if
it met the convergence criterion set by `get_digits`

.

The no-derivative method gives a result of

cbrt_noderiv(28) = 3.0365889718756618

with a 3 bits distance between the bracketed values, whereas the derivative methods both converge to a single value

cbrt_2deriv(28) = 3.0365889718756627

which we can compare with the boost::math::cbrt

cbrt(28) = 3.0365889718756627

Note that the iterations are set to stop at just one-half of full precision, and yet, even so, not one of the test cases had a single bit wrong. What's more, the maximum number of iterations was now just 4.

Just to complete the picture, we could have called `schroder_iterate`

in the last example:
and in fact it makes no difference to the accuracy or number of iterations
in this particular case. However, the relative performance of these two methods
may vary depending upon the nature of *f(x)*, and the
accuracy to which the initial guess can be computed. There appear to be no
generalisations that can be made except "try them and see".

Finally, had we called `cbrt`

with NTL::RR set to
1000 bit precision (about 300 decimal digits), then full precision can be
obtained with just 7 iterations. To put that in perspective, an increase
in precision by a factor of 20, has less than doubled the number of iterations.
That just goes to emphasise that most of the iterations are used up getting
the first few digits correct: after that these methods can churn out further
digits with remarkable efficiency.

Or to put it another way: *nothing beats a really good initial guess!*

Full code of this example is at root_finding_example.cpp,