Boost C++ Libraries of the most highly regarded and expertly designed C++ library projects in the world. Herb Sutter and Andrei Alexandrescu, C++ Coding Standards


Generalizing to Compute the nth root

If desired, we can now further generalize to compute the nth root by computing the derivatives at compile-time using the rules for differentiation and boost::math::pow<N> where template parameter N is an integer and a compile time constant. Our functor and function now have an additional template parameter N, for the root required.

[Note] Note

Since the powers and derivatives are fixed at compile time, the resulting code is as efficient as as if hand-coded as the cube and fifth-root examples above. A good compiler should also optimise any repeated multiplications.

Our nth root functor is

template <int N, class T = double>
struct nth_functor_2deriv
{ // Functor returning both 1st and 2nd derivatives.
  static_assert(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
  static_assert((N > 0) == true, "root N must be > 0!");

  nth_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of)
  { /* Constructor stores value a to find root of, for example: */ }

  // using boost::math::tuple; // to return three values.
  std::tuple<T, T, T> operator()(T const& x)
    // Return f(x), f'(x) and f''(x).
    using boost::math::pow;
    T fx = pow<N>(x) - a;                  // Difference (estimate x^n - a).
    T dx = N * pow<N - 1>(x);              // 1st derivative f'(x).
    T d2x = N * (N - 1) * pow<N - 2 >(x);  // 2nd derivative f''(x).

    return std::make_tuple(fx, dx, d2x);   // 'return' fx, dx and d2x.
  T a;                                     // to be 'nth_rooted'.

and our nth root function is

template <int N, class T = double>
T nth_2deriv(T x)
{ // return nth 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.

  static_assert(boost::is_integral<T>::value == false, "Only floating-point type types can be used!");
  static_assert((N > 0) == true, "root N must be > 0!");
  static_assert((N > 1000) == false, "root N is too big!");

  typedef double guess_type; // double may restrict (exponent) range for a multiprecision T?

  int exponent;
  frexp(static_cast<guess_type>(x), &exponent);                 // Get exponent of z (ignore mantissa).
  T guess = ldexp(static_cast<guess_type>(1.), exponent / N);   // Rough guess is to divide the exponent by n.
  T min = ldexp(static_cast<guess_type>(1.) / 2, exponent / N); // Minimum possible value is half our guess.
  T max = ldexp(static_cast<guess_type>(2.), exponent / N);     // Maximum possible value is twice our guess.

  int digits = std::numeric_limits<T>::digits * 0.4;            // Accuracy triples with each step, so stop when
                                                                // slightly more than one third of the digits are correct.
  const std::uintmax_t maxit = 20;
  std::uintmax_t it = maxit;
  T result = halley_iterate(nth_functor_2deriv<N, T>(x), guess, min, max, digits, it);
  return result;
    show_nth_root<5, double>(2.);
    show_nth_root<5, long double>(2.);
#ifndef _MSC_VER  // float128 is not supported by Microsoft compiler 2013.
    show_nth_root<5, float128>(2);
    show_nth_root<5, cpp_dec_float_50>(2); // dec
    show_nth_root<5, cpp_bin_float_50>(2); // bin

produces an output similar to this

 Using MSVC 2013

nth Root finding Example.
Type double value = 2, 5th root = 1.14869835499704
Type long double value = 2, 5th root = 1.14869835499704
Type class boost::multiprecision::number<class boost::multiprecision::backends::cpp_dec_float<50,int,void>,1> value = 2,
  5th root = 1.1486983549970350067986269467779275894438508890978
Type class boost::multiprecision::number<class boost::multiprecision::backends::cpp_bin_float<50,10,void,int,0,0>,0> value = 2,
  5th root = 1.1486983549970350067986269467779275894438508890978
[Tip] Tip

Take care with the type passed to the function. It is best to pass a double or greater-precision floating-point type.

Passing an integer value, for example, nth_2deriv<5>(2) will be rejected, while nth_2deriv<5, double>(2) converts the integer to double.

Avoid passing a float value that will provoke warnings (actually spurious) from the compiler about potential loss of data, as noted above.

[Warning] Warning

Asking for unreasonable roots, for example, show_nth_root<1000000>(2.); may lead to Loss of significance like Type double value = 2, 1000000th root = 1.00000069314783. Use of the the pow function is more sensible for this unusual need.

Full code of this example is at root_finding_n_example.cpp.