Let's now suppose we want to find the fifth root of a number a.

The equation we want to solve is :

  f(x) = x⁵ -a

If your differentiation is a little rusty (or you are faced with an function whose complexity makes differentiation daunting), then you can get help, for example, from the invaluable WolframAlpha site.

For example, entering the commmand: differentiate x ^ 5

or the Wolfram Language command: D[x ^ 5, x]

gives the output: d/dx(x ^ 5) = 5 x ^ 4

and to get the second differential, enter: second differentiate x ^ 5

or the Wolfram Language command: D[x ^ 5, { x, 2 }]

to get the output: d ^ 2 / dx ^ 2(x ^ 5) = 20 x ^ 3

To get a reference value, we can enter: fifth root 3126

or: N[3126 ^ (1 / 5), 50]

to get a result with a precision of 50 decimal digits:

5.0003199590478625588206333405631053401128722314376

(We could also get a reference value using multiprecision root).

The 1st and 2nd derivatives of x⁵ are:

  f'(x) = 5x⁴

  f''(x) = 20x³

Using these expressions for the derivatives, the functor is:

template <class T>
struct fifth_functor_2deriv
{
  // Functor returning both 1st and 2nd derivatives.
  fifth_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of)
  { /* Constructor stores value a to find root of, for example: */ }

  std::tuple<T, T, T> operator()(T const& x)
  {
    // Return both f(x) and f'(x) and f''(x).
    T fx = boost::math::pow<5>(x) - a;    // Difference (estimate x^3 - value).
    T dx = 5 * boost::math::pow<4>(x);    // 1st derivative = 5x^4.
    T d2x = 20 * boost::math::pow<3>(x);  // 2nd derivative = 20 x^3
    return std::make_tuple(fx, dx, d2x);  // 'return' fx, dx and d2x.
  }
private:
  T a;                                    // to be 'fifth_rooted'.
}; // struct fifth_functor_2deriv

Our fifth-root function is now:

template <class T>
T fifth_2deriv(T x)
{
  // return fifth root of x using 1st and 2nd derivatives and Halley.
  using namespace std;                  // Help ADL of std functions.
  using namespace boost::math::tools;   // for halley_iterate.

  int exponent;
  frexp(x, &exponent);                  // Get exponent of z (ignore mantissa).
  T guess = ldexp(1., exponent / 5);    // Rough guess is to divide the exponent by five.
  T min = ldexp(0.5, exponent / 5);     // Minimum possible value is half our guess.
  T max = ldexp(2., exponent / 5);      // Maximum possible value is twice our guess.
  // Stop when slightly more than one of the digits are correct:
  const int digits = static_cast<int>(std::numeric_limits<T>::digits * 0.4);
  const boost::uintmax_t maxit = 50;
  boost::uintmax_t it = maxit;
  T result = halley_iterate(fifth_functor_2deriv<T>(x), guess, min, max, digits, it);
  return result;
}

Full code of this example is at root_finding_example.cpp and root_finding_n_example.cpp.