Boost C++ Libraries

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

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

 Copyright © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle Walker and Xiaogang Zhang Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)