Boost C++ Libraries

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
	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.
  }
private:
  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);
#endif
    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
	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
	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.