Home | Libraries | People | FAQ | More |
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. BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!"); BOOST_STATIC_ASSERT_MSG((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. BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!"); BOOST_STATIC_ASSERT_MSG((N > 0) == true, "root N must be > 0!"); BOOST_STATIC_ASSERT_MSG((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 boost::uintmax_t maxit = 20; boost::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
Passing an integer value, for example,
Avoid passing a |
Warning | |
---|---|
Asking for unreasonable roots, for example, |
Full code of this example is at root_finding_n_example.cpp.