![]() |
Home | Libraries | People | FAQ | More |
First some #includes
that
will be needed.
#include <boost/math/tools/roots.hpp> //using boost::math::policies::policy; //using boost::math::tools::newton_raphson_iterate; //using boost::math::tools::halley_iterate; // //using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits. //using boost::math::tools::bracket_and_solve_root; //using boost::math::tools::toms748_solve; #include <boost/math/special_functions/next.hpp> // For float_distance. #include <tuple> // for std::tuple and std::make_tuple. #include <boost/math/special_functions/cbrt.hpp> // For boost::math::cbrt.
![]() |
Tip |
---|---|
For clarity, |
Let's suppose we want to find the root of a number a, and to start, compute the cube root.
So the equation we want to solve is:
f(x) = x³ -a
We will first solve this without using any information about the slope or curvature of the cube root function.
Fortunately, the cube-root function is 'Really Well Behaved' in that it is monotonic and has only one root (we leave negative values 'as an exercise for the student').
We then show how adding what we can know about this function, first just the slope or 1st derivative f'(x), will speed homing in on the solution.
Lastly, we show how adding the curvature f''(x) too will speed convergence even more.
First we define a function object (functor):
template <class T> struct cbrt_functor_noderiv { // cube root of x using only function - no derivatives. cbrt_functor_noderiv(T const& to_find_root_of) : a(to_find_root_of) { /* Constructor just stores value a to find root of. */ } T operator()(T const& x) { T fx = x*x*x - a; // Difference (estimate x^3 - a). return fx; } private: T a; // to be 'cube_rooted'. };
Implementing the cube-root function itself is fairly trivial now: the hardest part is finding a good approximation to begin with. In this case we'll just divide the exponent by three. (There are better but more complex guess algorithms used in 'real life'.)
template <class T> T cbrt_noderiv(T x) { // return cube root of x using bracket_and_solve (no derivatives). using namespace std; // Help ADL of std functions. using namespace boost::math::tools; // For bracket_and_solve_root. int exponent; frexp(x, &exponent); // Get exponent of z (ignore mantissa). T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three. T factor = 2; // How big steps to take when searching. const boost::uintmax_t maxit = 20; // Limit to maximum iterations. boost::uintmax_t it = maxit; // Initally our chosen max iterations, but updated with actual. bool is_rising = true; // So if result if guess^3 is too low, then try increasing guess. int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. // Some fraction of digits is used to control how accurate to try to make the result. int get_digits = digits - 3; // We have to have a non-zero interval at each step, so // maximum accuracy is digits - 1. But we also have to // allow for inaccuracy in f(x), otherwise the last few // iterations just thrash around. eps_tolerance<T> tol(get_digits); // Set the tolerance. std::pair<T, T> r = bracket_and_solve_root(cbrt_functor_noderiv<T>(x), guess, factor, is_rising, tol, it); return r.first + (r.second - r.first)/2; // Midway between brackets is our result, if necessary we could // return the result as an interval here. }
![]() |
Note |
---|---|
The final parameter specifying a maximum number of iterations is optional.
However, it defaults to So it may be wise to chose some reasonable estimate of how many iterations may be needed, In this case the function is so well behaved that we can chose a low value of 20.
Internally when Boost.Math uses these functions, it sets the maximum iterations
to |
Should we have wished we can show how many iterations were used in bracket_and_solve_root
(this information
is lost outside cbrt_noderiv
),
for example with:
if (it >= maxit) { std::cout << "Unable to locate solution in " << maxit << " iterations:" " Current best guess is between " << r.first << " and " << r.second << std::endl; } else { std::cout << "Converged after " << it << " (from maximum of " << maxit << " iterations)." << std::endl; }
for output like
Converged after 11 (from maximum of 20 iterations).
This snippet from main()
in root_finding_example.cpp
shows how it can be used.
try { double threecubed = 27.; // Value that has an *exactly representable* integer cube root. double threecubedp1 = 28.; // Value whose cube root is *not* exactly representable. std::cout << "cbrt(28) " << boost::math::cbrt(28.) << std::endl; // boost::math:: version of cbrt. std::cout << "std::cbrt(28) " << std::cbrt(28.) << std::endl; // std:: version of cbrt. std::cout <<" cast double " << static_cast<double>(3.0365889718756625194208095785056696355814539772481111) << std::endl; // Cube root using bracketing: double r = cbrt_noderiv(threecubed); std::cout << "cbrt_noderiv(" << threecubed << ") = " << r << std::endl; r = cbrt_noderiv(threecubedp1); std::cout << "cbrt_noderiv(" << threecubedp1 << ") = " << r << std::endl;
cbrt_noderiv(27) = 3 cbrt_noderiv(28) = 3.0365889718756618
The result of bracket_and_solve_root
is a pair
of values that could be displayed.
The number of bits separating them can be found using float_distance(r.first, r.second)
. The distance is zero (closest representable)
for 33 = 27 but float_distance(r.first, r.second) = 3
for cube root of 28 with this function. The result (avoiding overflow) is
midway between these two values.
We now solve the same problem, but using more information about the function, to show how this can speed up finding the best estimate of the root.
For the root function, the 1st differential (the slope of the tangent to a curve at any point) is known.
This algorithm is similar to this nth root algorithm.
If you need some reminders, then derivatives of elementary functions may help.
Using the rule that the derivative of xn for positive n (actually all nonzero n) is n xn-1, allows us to get the 1st differential as 3x2.
To see how this extra information is used to find a root, view Newton-Raphson iterations and the animation.
We define a better functor cbrt_functor_deriv
that returns both the evaluation of the function to solve, along with its
first derivative:
To 'return' two values, we use a std::pair of floating-point values.
template <class T> struct cbrt_functor_deriv { // Functor also returning 1st derivative. cbrt_functor_deriv(T const& to_find_root_of) : a(to_find_root_of) { // Constructor stores value a to find root of, // for example: calling cbrt_functor_deriv<T>(a) to use to get cube root of a. } std::pair<T, T> operator()(T const& x) { // Return both f(x) and f'(x). T fx = x*x*x - a; // Difference (estimate x^3 - value). T dx = 3 * x*x; // 1st derivative = 3x^2. return std::make_pair(fx, dx); // 'return' both fx and dx. } private: T a; // Store value to be 'cube_rooted'. };
Our cube root function is now:
template <class T> T cbrt_deriv(T x) { // return cube root of x using 1st derivative and Newton_Raphson. using namespace boost::math::tools; int exponent; frexp(x, &exponent); // Get exponent of z (ignore mantissa). T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three. T min = ldexp(0.5, exponent/3); // Minimum possible value is half our guess. T max = ldexp(2., exponent/3); // Maximum possible value is twice our guess. const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. int get_digits = static_cast<int>(digits * 0.6); // Accuracy doubles with each step, so stop when we have // just over half the digits correct. const boost::uintmax_t maxit = 20; boost::uintmax_t it = maxit; T result = newton_raphson_iterate(cbrt_functor_deriv<T>(x), guess, min, max, get_digits, it); return result; }
The result of newton_raphson_iterate
function
is a single value.
![]() |
Tip |
---|---|
There is a compromise between accuracy and speed when chosing the value
of |
Note that it is up to the caller of the function to check the iteration count after the call to see if iteration stoped as a result of running out of iterations rather than meeting the required precision.
Using the test data in /test/test_cbrt.cpp this found the cube root exact to the last digit in every case, and in no more than 6 iterations at double precision. However, you will note that a high precision was used in this example, exactly what was warned against earlier on in these docs! In this particular case it is possible to compute f(x) exactly and without undue cancellation error, so a high limit is not too much of an issue.
However, reducing the limit to std::numeric_limits<T>::digits * 2 / 3
gave full precision in all but one of the test cases (and that one was out
by just one bit). The maximum number of iterations remained 6, but in most
cases was reduced by one.
Note also that the above code omits a probable optimization by computing z² and reusing it, omits error handling, and does not handle negative values of z correctly. (These are left as the customary exercise for the reader!)
The boost::math::cbrt
function also includes these and other
improvements: most importantly it uses a much better initial guess which
reduces the iteration count to just 1 in almost all cases.
Next we define yet another even better functor cbrt_functor_2deriv
that returns both the evaluation of the function to solve, along with its
first and second derivative:
f''(x) = 6x
using information about both slope and curvature to speed convergence.
To 'return' three values, we use a tuple
of three floating-point values:
template <class T> struct cbrt_functor_2deriv { // Functor returning both 1st and 2nd derivatives. cbrt_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of) { // Constructor stores value a to find root of, for example: // calling cbrt_functor_2deriv<T>(x) to get cube root of x, } std::tuple<T, T, T> operator()(T const& x) { // Return both f(x) and f'(x) and f''(x). T fx = x*x*x - a; // Difference (estimate x^3 - value). T dx = 3 * x*x; // 1st derivative = 3x^2. T d2x = 6 * x; // 2nd derivative = 6x. return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x. } private: T a; // to be 'cube_rooted'. };
Our cube root function is now:
template <class T> T cbrt_2deriv(T x) { // return cube root of x using 1st and 2nd derivatives and Halley. //using namespace std; // Help ADL of std functions. using namespace boost::math::tools; int exponent; frexp(x, &exponent); // Get exponent of z (ignore mantissa). T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three. T min = ldexp(0.5, exponent/3); // Minimum possible value is half our guess. T max = ldexp(2., exponent/3); // Maximum possible value is twice our guess. const int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. // digits used to control how accurate to try to make the result. int get_digits = static_cast<int>(digits * 0.4); // Accuracy triples with each step, so stop when just // over one third of the digits are correct. boost::uintmax_t maxit = 20; T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, get_digits, maxit); return result; }
The function halley_iterate
also returns a single value, and the number of iterations will reveal if
it met the convergence criterion set by get_digits
.
The no-derivative method gives a result of
cbrt_noderiv(28) = 3.0365889718756618
with a 3 bits distance between the bracketed values, whereas the derivative methods both converge to a single value
cbrt_2deriv(28) = 3.0365889718756627
which we can compare with the boost::math::cbrt
cbrt(28) = 3.0365889718756627
Note that the iterations are set to stop at just one-half of full precision, and yet, even so, not one of the test cases had a single bit wrong. What's more, the maximum number of iterations was now just 4.
Just to complete the picture, we could have called schroder_iterate
in the last example:
and in fact it makes no difference to the accuracy or number of iterations
in this particular case. However, the relative performance of these two methods
may vary depending upon the nature of f(x), and the
accuracy to which the initial guess can be computed. There appear to be no
generalisations that can be made except "try them and see".
Finally, had we called cbrt
with NTL::RR set to
1000 bit precision (about 300 decimal digits), then full precision can be
obtained with just 7 iterations. To put that in perspective, an increase
in precision by a factor of 20, has less than doubled the number of iterations.
That just goes to emphasise that most of the iterations are used up getting
the first few digits correct: after that these methods can churn out further
digits with remarkable efficiency.
Or to put it another way: nothing beats a really good initial guess!
Full code of this example is at root_finding_example.cpp,