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, using statements are provided to list what functions are being used in this example: you can, of course, partly or fully qualify the names in other ways. (For your application, you may wish to extract some parts into header files, but you should never use using statements globally in header files).

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.

Cube root function without derivatives

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;                  // Initally 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 boost::uintmax_t maxit = (std::numeric_limits<boost::uintmax_t>::max)(); which is 18446744073709551615 and is more than anyone would wish to wait for! 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 policies::get_max_root_iterations<Policy>();.

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

Cube root function with 1st derivative (slope)

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ⁿ for positive n (actually all nonzero n) is n x^n-1, allows us to get the 1st differential as 3x².

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 chosing the value of digits. It is tempting to simply chose std::numeric_limits<T>::digits, but this may mean some inefficient and unnecessary iterations as the function thrashes around trying to locate the last bit. In theory, since the precision doubles with each step it is sufficient to stop when half the bits are correct: as the last step will have doubled that to full precision. Of course the function has no way to tell if that is actually the case unless it does one more step to be sure. In practice setting the precision to slightly more than std::numeric_limits<T>::digits / 2 is a good choice.

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.

Cube root with 1st & 2nd derivative (slope & curvature)

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,