- Series Evaluation diff --git a/doc/html/math_toolkit/internals1/cf.html b/doc/html/math_toolkit/internals1/cf.html index 2fc185849..d2f575c51 100644 --- a/doc/html/math_toolkit/internals1/cf.html +++ b/doc/html/math_toolkit/internals1/cf.html @@ -5,7 +5,7 @@ - + diff --git a/doc/html/math_toolkit/internals1/minima.html b/doc/html/math_toolkit/internals1/minima.html index 178d43f08..2eedcea2c 100644 --- a/doc/html/math_toolkit/internals1/minima.html +++ b/doc/html/math_toolkit/internals1/minima.html @@ -5,7 +5,7 @@ - + diff --git a/doc/html/math_toolkit/internals1/rational.html b/doc/html/math_toolkit/internals1/rational.html index 2740dff87..940056928 100644 --- a/doc/html/math_toolkit/internals1/rational.html +++ b/doc/html/math_toolkit/internals1/rational.html @@ -5,7 +5,7 @@ - + diff --git a/doc/html/math_toolkit/internals1/roots.html b/doc/html/math_toolkit/internals1/roots.html index 292b51f35..b5e73f0dd 100644 --- a/doc/html/math_toolkit/internals1/roots.html +++ b/doc/html/math_toolkit/internals1/roots.html @@ -5,7 +5,7 @@ - + @@ -61,7 +61,13 @@ Description
-
@@ -117,11 +123,13 @@
- int digits
- The desired number of binary digits. + The desired number of binary digits of precision.
- uintmax_t& max_iter
- An optional maximum number of iterations to perform. On exit this is + An optional maximum number of iterations to perform. Default
max_iteris(std::numeric_limits<boost::uintmax_t>::max)(),+ perhaps18446744073709551615+ which may be more than you want for wait for!
On exit this is set to the actual number of iterations performed.-
- Default max_iter =
(std::numeric_limits<boost::uintmax_t>::max)()is effectively 'iterate for ever'!. + Defaultmax_iter = + (std::numeric_limits<boost::uintmax_t>::max)()is effectively 'iterate for ever'!. - - They may be very sensitive to the initial guess, typically they converge - very rapidly if the initial guess has two or three decimal digits correct. - However convergence can be no better than bisection, or in some rare - cases, even worse than bisection if the initial guess is a long way from - the correct value and the derivatives are close to zero. + The functions may be very sensitive to the initial guess. Typically they + converge very rapidly if the initial guess has two or three decimal digits + correct. However, convergence can be no better than bisection, + or in some rare cases, even worse than bisection + if the initial guess is a long way from the correct value and the derivatives + are close to zero.
- These functions include special cases to handle zero first (and second - where appropriate) derivatives, and fall back to bisection in this case. - However, it is helpful if functor F is defined to return an arbitrarily - small value of the correct sign rather than zero. + where appropriate) derivatives, and fall back to bisection + in this case. However, it is helpful if functor F + is defined to return an arbitrarily small value of + the correct sign rather than zero.
- If the derivative at the current best guess for the result is infinite (or very close to being infinite) then these functions may terminate prematurely. A large first derivative leads to a very small next step, - triggering the termination condition. Derivative based iteration may + triggering the termination condition. Derivative-based iteration may not be appropriate in such cases.
-
If the function is 'Really Well Behaved' (monotonic and has only one
- root) the bracket bounds min and max may as well be set to the widest
- limits like zero and
numeric_limits<T>::max(). + root) the bracket bounds min and max + may as well be set to the widest limits like zero andnumeric_limits<T>::max(). - But if the function more complex and may have more than one root or a - pole, the choice of bounds is protection against jumping out to seek - the 'wrong' root. + pole, then the choice of bounds provides protection against jumping out + to seek the 'wrong' root.
- - These functions fall back to bisection if the next computed step would - take the next value out of bounds. The bounds are updated after each - step to ensure this leads to convergence. However, a good initial guess - backed up by asymptotically-tight bounds will improve performance no - end - rather than relying on bisection. + These functions fall back to bisection + if the next computed step would take the next value out of bounds. The + bounds are updated after each step to ensure this leads to convergence. + However, a good initial guess backed up by asymptotically-tight bounds + will improve performance no end - rather than relying on bisection.
- The value of digits is crucial to good performance - of these functions, if it is set too high then at best you will get one - extra (unnecessary) iteration, and at worst the last few steps will proceed - by bisection. Remember that the returned value can never be more accurate - than f(x) can be evaluated, and that if f(x) suffers from cancellation - errors as it tends to zero then the computed steps will be effectively - random. The value of digits should be set so that - iteration terminates before this point: remember that for second and - third order methods the number of correct digits in the result is increasing - quite substantially with each iteration, digits + of these functions, if it is set too high then, at best, you will get + one extra (unnecessary) iteration, and, at worst, the last few steps + will proceed by bisection. + Remember that the returned value can never be more accurate than f(x) + can be evaluated, and that if f(x) suffers from + cancellation errors as it tends to zero, then the computed steps will + be effectively random. The value of digits should + be set so that iteration terminates before this point: remember that + for second and third order methods the number of correct digits in the + result is increasing quite substantially with each iteration, digits should be set by experiment so that the final iteration just takes the - next value into the zone where f(x) becomes inaccurate. + next value into the zone where f(x) becomes inaccurate.
-
- To get the binary digits of accuracy, use policies::get_max_root_iterations<Policy>()).
+ To get the binary digits of accuracy, use
policies::get_max_root_iterations<Policy>()). -
If you need some diagnostic output to see what is going on, you can
@@ -208,13 +220,15 @@
Raphson Method
- Given an initial guess x0 the subsequent values are computed using: + Given an initial guess x0 the subsequent values are + computed using:
- Out of bounds steps revert to bisection of the current bounds. + Out of bounds steps revert to bisection + of the current bounds.
Under ideal conditions, the number of correct digits doubles with each iteration. @@ -235,7 +249,8 @@ wrong direction) causes the method to revert to a Newton-Raphson step.
- Out of bounds steps revert to bisection of the current bounds. + Out-of-bounds steps revert to bisection + of the current bounds.
Under ideal conditions, the number of correct digits trebles with each iteration. @@ -246,7 +261,8 @@ Method
- Given an initial guess x0 the subsequent values are computed using: + Given an initial guess x0 the subsequent values are + computed using:
@@ -258,7 +274,8 @@
by more than 10%.
- Out of bounds steps revert to bisection of the current bounds. + Out of bounds steps revert to bisection + of the current bounds.
Under ideal conditions, the number of correct digits trebles with each iteration. @@ -275,90 +292,140 @@
- To begin with lets solve the problem using Newton-Raphson iterations, we'll - begin by defining a function object (functor) that returns the evaluation - of the function to solve, along with its first derivative f'(x): + To begin with let us solve the problem using Newton-Raphson iterations. We + will begin by defining a function object (functor) that returns the evaluation + of the function to solve, along with its first derivative f'(x):
template <class T> -struct cbrt_functor -{ - cbrt_functor(T const& target) : a(target) - { // Constructor stores value to be 'cube-rooted'. - } - boost::math::tuple<T, T> operator()(T const& z) - { // z is estimate so far. - return boost::math::make_tuple( - z*z*z - a, // return both f(x) - 3 * z*z); // and f'(x) - } +struct cbrt_functor_1stderiv +{ // Functor returning function and 1st derivative. + + cbrt_functor_1stderiv(T const& target) : value(target) + { // Constructor stores the value to be 'cube_rooted'. + } + + std::pair<T, T> operator()(T const& z) // z is best estimate so far. + { // Return both f(x) and derivative f'(x). + T fx = z*z*z - value; // Difference estimate fx = x^3 - value. + T d1x = 3 * z*z; // 1st derivative d1x = 3x^2. + return std::make_pair(fx, d1x); // 'return' both fx and d1x. + } private: - T a; // to be 'cube-rooted'. -}; + T value; // to be 'cube_rooted'. +}; // cbrt_functor_1stderiv
Implementing the cube root 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: + 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(T z) -{ - using namespace std; // for frexp, ldexp, numeric_limits. - using namespace boost::math::tools; +T cbrt_1deriv(T x) +{ //! \return cube root of x using 1st derivative and Newton_Raphson. + using namespace std; // for frexp, ldexp, numeric_limits. + using namespace boost::math::tools; - int exp; - frexp(z, &exp); // Get exponent of z (ignore mantissa). - T min = ldexp(0.5, exp/3); - T max = ldexp(2.0, exp/3); - T guess = ldexp(1.0, exp/3); // Rough guess is to divide the exponent by three. - int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. - return newton_raphson_iterate(detail::cbrt_functor<T>(z), guess, min, max, digits); -} + int exponent; + frexp(x, &exponent); // Get exponent of x (ignore mantissa). + T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three. + // Set an initial bracket interval. + 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. + // digits used to control how accurate to try to make the result. + int digits = std::numeric_limits<T>::digits; // Maximum possible binary digits accuracy for type T. + + boost::uintmax_t maxit = 20; // Optionally limit the number of iterations. + //cout << "Max Iterations " << maxit << endl; // + T result = newton_raphson_iterate(cbrt_functor_1stderiv<T>(x), guess, min, max, digits, maxit); + // Can check and show how many iterations (updated by newton_raphson_iterate). + cout << "Iterations " << maxit << endl; + return result; +} // cbrt_1deriv
- Using the test data in
+libs/math/test/cbrt_test.cppthis + Using the test data in ../../test/cbrt_data.ipp, + the test ../../test/test_cbrt.cpp 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 tostd::numeric_limits<T>::digits * 2 / 3+ than 6 iterations at double precision. ++ However, you will note that a high (maximum) 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 / 3gave 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.+ So we might decide to trade the accuracy and limit iterations for speed. +
+int get_digits = (digits * 2) /3; // Two thirds of maximum possible accuracy. +boost::uintmax_t maxit = 7; // Limit the number of iterations. +
Note also that the above code omits a probably optimization by computing - z², and reusing it, omits error handling, and does not handle negative values - of z correctly. (These are left as an exercise for the reader!) + z², and reusing it, omits error handling, and does not + handle negative values of z correctly. (As is traditional, these are left + as an exercise for the reader!)
- The
boost::math::cbrtfunction also includes these and other - improvements. + Theboost::math::cbrtfunction cbrt + also includes these and other improvements.- Now let's adapt the functor slightly to return the second derivative as well: + Now let's adapt the functor slightly to return the second derivative f''(x) + as well:
template <class T> -struct cbrt_functor -{ - cbrt_functor(T const& target) : a(target){} - boost::math::tuple<T, T, T> operator()(T const& z) - { - return boost::math::make_tuple( - z*z*z - a, - 3 * z*z, - 6 * z); - } +struct cbrt_functor_2deriv +{ // Functor returning both 1st and 2nd derivatives. + cbrt_functor_2deriv(T const& to_find_root_of) : value(to_find_root_of) + { // Constructor stores value to find root of, for example: + } + + // using boost::math::tuple; // to return three values. + boost::math::tuple<T, T, T> operator()(T const& x) + { // Return both f(x) and f'(x) and f''(x). + using boost::math::make_tuple; + T fx = x*x*x - value; // Difference (estimate x^3 - value). + T dx = 3 * x*x; // 1st derivative = 3x^2. + T d2x = 6 * x; // 2nd derivative = 6x. + return make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x. + } private: - T a; -}; + T value; // to be 'cube_rooted'. +}; // struct cbrt_functor_2deriv
And then adapt the
cbrtfunction to use Halley iterations: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; + 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. + + int digits = std::numeric_limits<T>::digits /2 ; // Half maximum possible binary digits accuracy for type T. + + //boost::uintmax_t maxit = (std::numeric_limits<boost::uintmax_t>::max)(); + // the default is (std::numeric_limits<boost::uintmax_t>::max)() = 18446744073709551615 + // which is more than we might wish to wait for!!! so we can reduce it. + boost::uintmax_t maxit = 10; + T result = halley_iterate(cbrt_functor_2deriv<T>(x), guess, min, max, digits, maxit); + // Can show how many iterations (updated by newton_raphson_iterate). + cout << "Iterations " << maxit << endl; + return result; +} // cbrt_2deriv(x) +
+template <class T> T cbrt(T z) { using namespace std; @@ -382,21 +449,27 @@ Just to complete the picture, we could have called
schroeder_iteratein 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". + 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, 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 + Finally, had we called cbrt + with Boost.Multiprecision + or NTL::RR set to 1000 + bit precision, then full precision still 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! + Or to put it another way: nothing beats a really good + initial guess! +
++ See root_finding_example.cpp + for full source code.
diff --git a/doc/html/math_toolkit/internals1/roots2.html b/doc/html/math_toolkit/internals1/roots2.html index aa2817c8b..4ab8b112c 100644 --- a/doc/html/math_toolkit/internals1/roots2.html +++ b/doc/html/math_toolkit/internals1/roots2.html @@ -5,7 +5,7 @@ - + @@ -143,19 +143,47 @@ Description - These functions solve the root of some function f(x) - without the need for any derivatives of f(x). + These functions
+bisect, + bracket_and_solve_root, andtoms748_solvesolve the root of some function + f(x) without the need for any + derivatives of f(x).+ (Algorithms that can use the derivative of f(x) are + nearly always faster and root-finding + with derivatives should be used whenever possible). +
++ Variants of these functions allow control via a tolerance for stopping iteration, + setting a maximum number of iterations, or provision of a guess. +
++ Variants with and without Policy are provided + to allow: +
+-
+
- + Control of the precision used via precision + policy. + +
- + Handling of errors via error + handling policies. + +
The
bracket_and_solve_root- functions use TOMS Algorithm 748 that is asymptotically the most efficient - known, and have been shown to be optimal for a certain classes of smooth - functions. Variants with and without __policies are provided. + functions use TOMS + Algorithm 748: enclosing zeros of continuous functions that is asymptotically + the most efficient known, and have been shown to be optimal for a certain + classes of smooth functions.- Alternatively, there is a simple bisection routine which can be useful in - its own right in some situations, or alternatively for narrowing down the - range containing the root, prior to calling a more advanced algorithm. + Alternatively, there is a simple bisection + routine
bisectwhich can + be useful in its own right in some situations, or alternatively for narrowing + down the range containing the root, to get a good guess prior to calling + a more advanced algorithm.All the algorithms in this section reduce the diameter of the enclosing interval @@ -165,15 +193,28 @@ the root. This is also in contrast to some other derivative-free methods (for example the methods of Brent or Dekker) which only reduce the enclosing interval on the final - step. Therefore these methods return a std::pair containing the enclosing - interval found, and accept a function object specifying the termination condition. - Three function objects are provided for ready-made termination conditions: - eps_tolerance causes termination when the relative error - in the enclosing interval is below a certain threshold, while equal_floor - and equal_ceil are useful for certain statistical applications - where the result is known to be an integer. Other user-defined termination - conditions are likely to be used only rarely, but may be useful in some specific - circumstances. + step. +
++ Therefore, these methods return a
+std::pair+ containing the enclosing interval found, and accept a function object specifying + the termination condition. Three function objects are provided for ready-made + termination conditions: ++-
+
- + eps_tolerance causes termination when the relative + error in the enclosing interval is below a certain threshold, + +
- + equal_floor, equal_ceil and + equal_nearest_integer are useful for certain statistical + applications where the result is known to be an integer. + +
+ Other user-defined termination conditions are likely to be used only rarely, + but might be useful in some specific circumstances.
@@ -207,7 +248,8 @@ const Policy&);
- These functions locate the root using bisection: function arguments are: + These functions locate the root using bisection: + their function arguments are:
@@ -302,24 +346,36 @@ const Policy&);@@ -218,27 +260,29 @@
- min
- The left bracket of the interval known to contain the root. + The left bracket of an interval known to contain the root.
- max
- The right bracket of the interval known to contain the root. It is - a precondition that min < max and f(min)*f(max) - <= 0, the function signals evaluation error if these - preconditions are violated. The action taken is controlled by the evaluation - error policy. A best guess may be returned, perhaps significantly wrong. + The right bracket of an interval known to contain the root. It is a + precondition that min < max and f(min) + * f(max) <= 0. The function signals evaluation_error + if these preconditions are violated. The action taken is controlled + by the error + handling policies. A best guess may be returned, perhaps significantly + wrong.
- tol
A binary functor that specifies the termination condition: the function will return the current brackets enclosing the root when tol(min,max) - becomes true. + becomes
true.- max_iter
The maximum number of invocations of f(x) to make - while searching for the root. On exist this is set to actual number - of invocations performed. + while searching for the root.
On exit this is set to actual + number of invocations performed.- This is a convenience function that calls toms748_solve - internally to find the root of f(x). It's usable only - when f(x) is a monotonic function, and the location - of the root is known approximately, and in particular it is known whether - the root is occurs for positive or negative x. The parameters - are: +
+bracket_and_solve_rootis + a convenience function that calls TOMS + Algorithm 748: enclosing zeros of continuous functions internally + to find the root of f(x). It's usable only when: ++-
+
- + Derivative of f(x) is a monotonic function. + +
- + The location of the root is known approximately (a guess). + +
- + It is known whether the root occurs for positive or negative x. + +
+ The parameters are:
@@ -431,10 +488,16 @@ const Policy&);- f
- A unary functor that is the function whose root is to be solved. f(x) + A unary functor that is the function whose root is to be solved. f(x) must be uniformly increasing or decreasing on x.
- guess
- An initial approximation to the root + An initial approximation to the root.
- factor
@@ -340,13 +396,14 @@
A binary functor that determines the termination condition for the search for the root. tol is passed the current - brackets at each step, when it returns true then the current brackets - are returned as the result. + brackets at each step, when it returns
true+ then the current brackets are returned as the result.- max_iter
The maximum number of function invocations to perform in the search - for the root. On exit is set to the actual number of invocations performed. + for the root.
On exit max_iter is set to + the actual number of invocations performed.
- These two functions implement TOMS Algorithm 748: it uses a mixture of cubic, - quadratic and linear (secant) interpolation to locate the root of f(x). - The two functions differ only by whether values for f(a) - and f(b) are already available. The toms748_solve parameters + The functions
+toms748_solve+ implement TOMS + Algorithm 748: enclosing zeros of continuous functions : using a + mixture of cubic, quadratic and linear (secant) interpolation to locate the + root of f(x). Function variants differ only by whether + values for f(a) and f(b) are already + available, and if special policies apply. ++ The
toms748_solveparameters are:@@ -454,8 +517,8 @@@@ -487,10 +550,24 @@ documentation for more details.The upper bound for the initial bracket of the root. It is a precondition that a < b and that a - and b bracket the root to find so that f(a)*f(b) - < 0. + and b bracket the root to find so that f(a) + * f(b) < 0.
- fa
@@ -469,14 +532,14 @@
A binary functor that determines the termination condition for the search for the root. tol is passed the current - brackets at each step, when it returns true, then the current brackets - are returned as the result. + brackets at each step, when it returns
true, + then the current brackets are returned as the result.- max_iter
The maximum number of function invocations to perform in the search - for the root. On exit max_iter is set to actual - number of function invocations used. + for the root.
On exit max_iter is set to + actual number of function invocations used.- toms748_solve returns: a pair of values r that bracket - the root so that: f(r.first) * f(r.second) <= 0 and either tol(r.first, - r.second) == true or max_iter >= m where m is the - initial value of max_iter passed to the function. +
+toms748_solvereturns a pair + of values r that bracket the root so that: +f(r.first) * f(r.second) <= 0 +
++ and either +
+tol(r.first, r.second) == true +
++ or +
+max_iter >= m +
++ where m is the initial value of max_iter + passed to the function.
In other words, it's up to the caller to verify whether termination occurred @@ -508,13 +585,18 @@
eps_toleranceis the usual - termination condition used with these root finding functions. Its operator() + termination condition used with these root-finding functions. Itsoperator()will return true when the relative distance between a - and b is less than twice the machine epsilon for T, - or 21-bits, whichever is the larger. In other words, you set bits - to the number of bits of precision you want in the result. The minimal tolerance - of twice the machine epsilon of T is required to ensure that we get back - a bracketing interval: since this must clearly be at least 1 epsilon in size. + and b is less than twice the machine + epsilon for T, or 21-bits, whichever is the larger. In other words, + you set bits to the number of bits of precision you + want in the result. The minimal tolerance of twice the machine + epsilon ofT, typically +2 * + std::numeric_limits<T>::epsilon, is required to ensure that we + get back a bracketing interval > zero, since this must clearly be at least + 1 machine epsilon + in size.struct equal_floor { diff --git a/doc/html/math_toolkit/internals1/series_evaluation.html b/doc/html/math_toolkit/internals1/series_evaluation.html index 784afcf91..1e32b62e7 100644 --- a/doc/html/math_toolkit/internals1/series_evaluation.html +++ b/doc/html/math_toolkit/internals1/series_evaluation.html @@ -5,8 +5,8 @@ - - + + diff --git a/doc/html/math_toolkit/internals1/tuples.html b/doc/html/math_toolkit/internals1/tuples.html index e1f697033..32c22edd2 100644 --- a/doc/html/math_toolkit/internals1/tuples.html +++ b/doc/html/math_toolkit/internals1/tuples.html @@ -5,7 +5,7 @@ - + diff --git a/doc/html/math_toolkit/internals2.html b/doc/html/math_toolkit/internals2.html index 429c7d996..929641028 100644 --- a/doc/html/math_toolkit/internals2.html +++ b/doc/html/math_toolkit/internals2.html @@ -5,7 +5,7 @@ - + diff --git a/doc/html/math_toolkit/internals_overview.html b/doc/html/math_toolkit/internals_overview.html index 69ef023d4..66e9f0570 100644 --- a/doc/html/math_toolkit/internals_overview.html +++ b/doc/html/math_toolkit/internals_overview.html @@ -5,9 +5,9 @@ - - - + + +
@@ -28,8 +28,9 @@ This section contains internal utilities used by the library's implementation - along with tools used in development and testing. These tools have only minimal - documentation, and crucially do not have stable interfaces. + along with tools used in development and testing. These tools have limitied + documentation, but now have quite stable interfaces and may also be useful + outside Boost.Math.
There is no doubt that these components can be improved, but they are also diff --git a/doc/html/math_toolkit/main_intro.html b/doc/html/math_toolkit/main_intro.html index 8d4ef04ed..a8a701faa 100644 --- a/doc/html/math_toolkit/main_intro.html +++ b/doc/html/math_toolkit/main_intro.html @@ -87,7 +87,7 @@ Toolkit
- Provides many of the tools required to implement + Provides many of the tools required to implement mathematical special functions: hopefully the presence of these will encourage other authors to contribute more special function implementations in the future. These tools are currently considered experimental: they are "exposed implementation diff --git a/doc/html/math_toolkit/navigation.html b/doc/html/math_toolkit/navigation.html index ff15a1729..3675488b4 100644 --- a/doc/html/math_toolkit/navigation.html +++ b/doc/html/math_toolkit/navigation.html @@ -27,7 +27,7 @@ Navigation
Boost.Math documentation is provided in both HTML and PDF formats. diff --git a/doc/html/math_toolkit/special_tut/special_tut_impl.html b/doc/html/math_toolkit/special_tut/special_tut_impl.html index 136ca0814..6c3c459fd 100644 --- a/doc/html/math_toolkit/special_tut/special_tut_impl.html +++ b/doc/html/math_toolkit/special_tut/special_tut_impl.html @@ -182,7 +182,7 @@
- Where possible, evaluation of series, continued fractions, polynomials, - or root finding should use one of the boiler + or root finding should use one of the boiler plate functions. In any case, after any iterative method, you should verify that the number of iterations did not exceed the maximum specified in the Policy type, and if it diff --git a/doc/html/toolkit.html b/doc/html/toolkit.html index 3d920f616..50a16988a 100644 --- a/doc/html/toolkit.html +++ b/doc/html/toolkit.html @@ -1,7 +1,7 @@ -
Chapter 12. Internals (Series, Rationals and Continued Fractions, Root Finding, Function Minimization, Testing and Development Tools) +Chapter 12. Series, Rationals and Continued Fractions, Root Finding, Function Minimization, Testing, and Development Tools @@ -24,13 +24,13 @@+Chapter 12. Series, Rationals and Continued Fractions, Root Finding, Function Minimization, + Testing, and Development ToolsTable of Contents
- Overview -
- Utilities & internals +
- Utilities
- Series Evaluation
- Continued Fraction Evaluation diff --git a/doc/internals/internals_overview.qbk b/doc/internals/internals_overview.qbk index b49bebf96..b0327f62f 100644 --- a/doc/internals/internals_overview.qbk +++ b/doc/internals/internals_overview.qbk @@ -2,7 +2,8 @@ This section contains internal utilities used by the library's implementation along with tools used in development and testing. These tools have -only minimal documentation, and crucially ['do not have stable interfaces]. +limitied documentation, but now have quite stable interfaces +and may also be useful outside Boost.Math. There is no doubt that these components can be improved, but they are also largely incidental to the main purpose of this library. @@ -14,7 +15,7 @@ submissions to this library. [endsect] [/section:internals_overview Overview] [/ - Copyright 2006, 2010, 2013 John Maddock and Paul A. Bristow. + Copyright 2006, 2010, 2013, 2014 John Maddock and Paul A. Bristow. 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). diff --git a/doc/internals/roots.qbk b/doc/internals/roots.qbk index fb0950926..f6a178607 100644 --- a/doc/internals/roots.qbk +++ b/doc/internals/roots.qbk @@ -1,5 +1,7 @@ [section:roots Root Finding With Derivatives: Newton-Raphson, Halley & Schroeder] +[import ../../example/root_finding_example.cpp] + [h4 Synopsis] `` @@ -31,7 +33,11 @@ [h4 Description] -These functions all perform iterative root finding using derivatives: +These functions `newton_raphson_iterate, halley_iterate`, and `schroeder_iterate` +all perform iterative root finding using derivatives. + +(These functions should be used whenever derivative(s) are available, +otherwise use __root_finding_without_derivatives). * `newton_raphson_iterate` performs second order [link math_toolkit.internals1.roots.newton Newton-Raphson iteration], @@ -46,7 +52,7 @@ The functions all take the same parameters: returns a __tuple: For the second order iterative methods ([@http://en.wikipedia.org/wiki/Newton_Raphson Newton Raphson]) - the __tuple should have *two* elements containing the evaluation + the __tuple should have [*two] elements containing the evaluation of the function and its first derivative. For the third order methods @@ -57,46 +63,49 @@ Schroeder) [[T guess] [The initial starting value. A good guess is crucial to quick convergence!]] [[T min] [The minimum possible value for the result, this is used as an initial lower bracket.]] [[T max] [The maximum possible value for the result, this is used as an initial upper bracket.]] -[[int digits] [The desired number of binary digits.]] -[[uintmax_t& max_iter] [An optional maximum number of iterations to perform. On exit this is set to the actual number of iterations performed.]] +[[int digits] [The desired number of binary digits of precision.]] +[[uintmax_t& max_iter] [An optional maximum number of iterations to perform. +Default `max_iter` is `(std::numeric_limits
::max()`. *But if the function more complex and may have more than one root or a pole, -the choice of bounds is protection against jumping out to seek the 'wrong' root. -* These functions fall back to bisection if the next computed step would take the +then the choice of bounds provides protection against jumping out to seek the 'wrong' root. +* These functions fall back to __bisection if the next computed step would take the next value out of bounds. The bounds are updated after each step to ensure this leads to convergence. However, a good initial guess backed up by asymptotically-tight -bounds will improve performance no end - rather than relying on bisection. -* The value of /digits/ is crucial to good performance of these functions, -if it is set too high then at best you will get one extra (unnecessary) -iteration, and at worst the last few steps will proceed by bisection. -Remember that the returned value can never be more accurate than f(x) can be -evaluated, and that if f(x) suffers from cancellation errors as it -tends to zero then the computed steps will be effectively random. The -value of /digits/ should be set so that iteration terminates before this point: +bounds will improve performance no end - rather than relying on __bisection. +* The value of ['digits] is crucial to good performance of these functions, +if it is set too high then, at best, you will get one extra (unnecessary) +iteration, and, at worst, the last few steps will proceed by __bisection. +Remember that the returned value can never be more accurate than ['f(x)] can be +evaluated, and that if ['f(x)] suffers from cancellation errors as it +tends to zero, then the computed steps will be effectively random. The +value of ['digits] should be set so that iteration terminates before this point: remember that for second and third order methods the number of correct -digits in the result is increasing quite -substantially with each iteration, /digits/ should be set by experiment so that the final -iteration just takes the next value into the zone where f(x) becomes inaccurate. -* To get the binary digits of accuracy, use policies::get_max_root_iterations ()). +digits in the result is increasing quite substantially with each iteration, +['digits] should be set by experiment so that the final +iteration just takes the next value into the zone where ['f(x)] becomes inaccurate. +* To get the binary digits of accuracy, use `policies::get_max_root_iterations ()`). * If you need some diagnostic output to see what is going on, you can `#define BOOST_MATH_INSTRUMENT` before the `#include `, and also ensure that display of all the possibly significant digits with @@ -109,11 +118,11 @@ than computing the derivatives themselves. As ever, algebraic simplification ca be a big win. [h4:newton Newton Raphson Method] -Given an initial guess x0 the subsequent values are computed using: +Given an initial guess ['x0] the subsequent values are computed using: [equation roots1] -Out of bounds steps revert to bisection of the current bounds. +Out of bounds steps revert to __bisection of the current bounds. Under ideal conditions, the number of correct digits doubles with each iteration. @@ -127,13 +136,13 @@ Over-compensation by the second derivative (one which would proceed in the wrong direction) causes the method to revert to a Newton-Raphson step. -Out of bounds steps revert to bisection of the current bounds. +Out-of-bounds steps revert to __bisection of the current bounds. Under ideal conditions, the number of correct digits trebles with each iteration. [h4:schroeder Schroeder's Method] -Given an initial guess x0 the subsequent values are computed using: +Given an initial guess ['x0] the subsequent values are computed using: [equation roots3] @@ -142,7 +151,7 @@ in the wrong direction) causes the method to revert to a Newton-Raphson step. Likewise a Newton step is used whenever that Newton step would change the next value by more than 10%. -Out of bounds steps revert to bisection of the current bounds. +Out of bounds steps revert to __bisection of the current bounds. Under ideal conditions, the number of correct digits trebles with each iteration. @@ -155,78 +164,47 @@ solve along with its derivatives are: To begin with let us solve the problem using Newton-Raphson iterations. We will begin by defining a function object (functor) that returns the evaluation -of the function to solve, along with its first derivative f'(x): +of the function to solve, along with its first derivative ['f'(x)]: - template - struct cbrt_functor - { - cbrt_functor(T const& target) : a(target) - { // Constructor stores value to be 'cube-rooted'. - } - ``__tuple`` - T cbrt(T z) - { - using namespace std; // for frexp, ldexp, numeric_limits. - using namespace boost::math::tools; +[root_finding_cbrt_1deriv] - int exp; - frexp(z, &exp); // Get exponent of z (ignore mantissa). - T min = ldexp(0.5, exp/3); - T max = ldexp(2.0, exp/3); - T guess = ldexp(1.0, exp/3); // Rough guess is to divide the exponent by three. - int digits = std::numeric_limits ::digits; // Maximum possible binary digits accuracy for type T. - return newton_raphson_iterate(detail::cbrt_functor (z), guess, min, max, digits); - } +Using the test data in [@../../test/cbrt_data.ipp], the test [@../../test/test_cbrt.cpp] +found the cube root exact to the last digit in every case, +and in no more than 6 iterations at double precision. -Using the test data in `libs/math/test/cbrt_test.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 +However, you will note that a high (maximum) 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 +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 ::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 probably optimization by computing z[sup2], -and reusing it, omits error handling, and does not handle -negative values of z correctly. (These are left as an exercise for the reader!) +So we might decide to trade the accuracy and limit iterations for speed. -The `boost::math::cbrt` function also includes these and other improvements. + int get_digits = (digits * 2) /3; // Two thirds of maximum possible accuracy. + boost::uintmax_t maxit = 7; // Limit the number of iterations. -Now let's adapt the functor slightly to return the second derivative as well: +Note also that the above code omits a probably optimization by computing ['z[sup2]], +and reusing it, omits error handling, and does not handle negative values of z correctly. +(As is traditional, these are left as an exercise for the reader!) - template - struct cbrt_functor - { - cbrt_functor(T const& target) : a(target){} - ``__tuple`` T cbrt(T z) { @@ -249,19 +227,23 @@ What's more, the maximum number of iterations was now just 4. Just to complete the picture, we could have called `schroeder_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 +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 [@http://shoup.net/ntl/doc/RR.txt NTL::RR] -set to 1000 bit precision, then full precision can be obtained with just 7 iterations. +Finally, had we called __cbrt with +__multiprecision or [@http://shoup.net/ntl/doc/RR.txt NTL::RR] +set to 1000 bit precision, then full precision still 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 +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!] +Or to put it another way: [*nothing beats a really good initial guess!] + +See [@../../example/root_finding_example.cpp root_finding_example.cpp] for full source code. + [endsect] [/section:roots Root Finding With Derivatives] diff --git a/doc/internals/roots_without_derivatives.qbk b/doc/internals/roots_without_derivatives.qbk index d1ad6f399..ee1f3fea3 100644 --- a/doc/internals/roots_without_derivatives.qbk +++ b/doc/internals/roots_without_derivatives.qbk @@ -5,99 +5,99 @@ `` #include `` - + namespace boost{ namespace math{ namespace tools{ // Bisection template - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair - std::pair struct eps_tolerance @@ -364,25 +391,26 @@ rather than because the termination condition /tol/ was satisfied. eps_tolerance(int bits); bool operator()(const T& a, const T& b)const; }; - -`eps_tolerance` is the usual termination condition used with these root finding functions. -Its operator() will return true when the relative distance between /a/ and /b/ -is less than twice the machine epsilon for T, or 2[super 1-bits], whichever is -the larger. In other words, you set /bits/ to the number of bits of precision you -want in the result. The minimal tolerance of twice the machine epsilon of T is -required to ensure that we get back a bracketing interval: since this must clearly -be at least 1 epsilon in size. - + +`eps_tolerance` is the usual termination condition used with these root-finding functions. +Its `operator()` will return true when the relative distance between ['a] and ['b] +is less than twice the __epsilon for T, or 2[super 1-bits], whichever is +the larger. In other words, you set ['bits] to the number of bits of precision you +want in the result. The minimal tolerance of twice the __epsilon of `T`, +typically `2 * std::numeric_limits ::epsilon`, is +required to ensure that we get back a bracketing interval > zero, since this must clearly +be at least 1 __epsilon in size. + struct equal_floor { equal_floor(); template bool operator()(const T& a, const T& b)const; }; - + This termination condition is used when you want to find an integer result -that is the /floor/ of the true root. It will terminate as soon as both ends -of the interval have the same /floor/. - +that is the ['floor] of the true root. It will terminate as soon as both ends +of the interval have the same ['floor]. + struct equal_ceil { equal_ceil(); @@ -390,8 +418,8 @@ of the interval have the same /floor/. }; This termination condition is used when you want to find an integer result -that is the /ceil/ of the true root. It will terminate as soon as both ends -of the interval have the same /ceil/. +that is the ['ceil] of the true root. It will terminate as soon as both ends +of the interval have the same ['ceil]. struct equal_nearest_integer { @@ -400,7 +428,7 @@ of the interval have the same /ceil/. }; This termination condition is used when you want to find an integer result -that is the /closest/ to the true root. It will terminate as soon as both ends +that is the ['closest] to the true root. It will terminate as soon as both ends of the interval round to the same nearest integer. [h4 Implementation] @@ -408,17 +436,17 @@ of the interval round to the same nearest integer. The implementation of the bisection algorithm is extremely straightforward and not detailed here. __TOMS748 is described in detail in: -['Algorithm 748: Enclosing Zeros of Continuous Functions, -G. E. Alefeld, F. A. Potra and Yixun Shi, -ACM Transactions on Mathematica1 Software, Vol. 21. No. 3. September 1995. +['Algorithm 748: Enclosing Zeros of Continuous Functions, +G. E. Alefeld, F. A. Potra and Yixun Shi, +ACM Transactions on Mathematica1 Software, Vol. 21. No. 3. September 1995. Pages 327-344.] The implementation here is a faithful translation of this paper into C++. [endsect] [/section:roots2 Root Finding Without Derivatives] -[/ - Copyright 2006, 2010 John Maddock and Paul A. Bristow. +[/ + Copyright 2006, 2010, 2014 John Maddock and Paul A. Bristow. 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). diff --git a/doc/math.qbk b/doc/math.qbk index 3459816ef..0d275b6ad 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -93,9 +93,10 @@ and use the function's name as the link text.] [def __lanczos [link math_toolkit.lanczos Lanczos approximation]] [def __zero_error [link math_toolkit.relative_error.zero_error effectively zero error]] [def __relative_error [link math_toolkit.relative_error relative zero error]] +[/ links to roots] [def __root_finding_with_derivatives [link math_toolkit.internals1.roots root-finding with derivatives]] [def __root_finding_without_derivatives [link math_toolkit.internals1.roots2 root-finding without derivatives]] - +[def __bisection [@http://en.wikipedia.org/wiki/Bisection_method bisection]] [/gammas] [def __lgamma [link math_toolkit.sf_gamma.lgamma lgamma]] [def __digamma [link math_toolkit.sf_gamma.digamma digamma]] @@ -546,11 +547,11 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include octonion/math-octonion.qbk] [include gcd/math-gcd.qbk] -[mathpart toolkit Internals (Series, Rationals and Continued Fractions, Root Finding, Function Minimization, Testing and Development Tools)] +[mathpart toolkit Series, Rationals and Continued Fractions, Root Finding, Function Minimization, Testing, and Development Tools] [include internals/internals_overview.qbk] -[section:internals1 Utilities & internals] +[section:internals1 Utilities] [include internals/series.qbk] [include internals/fraction.qbk] [include internals/rational.qbk] @@ -558,7 +559,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include internals/roots_without_derivatives.qbk] [include internals/minima.qbk] [include internals/tuple.qbk] -[endsect] [/section:internals1 Utilities & internals] +[endsect] [/section:internals1 Utilities] [section:internals2 Testing and Development] [include internals/polynomial.qbk] diff --git a/example/root_finding_example.cpp b/example/root_finding_example.cpp index 5829d981c..e10d1d420 100644 --- a/example/root_finding_example.cpp +++ b/example/root_finding_example.cpp @@ -1,25 +1,22 @@ // root_finding_example.cpp -// Copyright Paul A. Bristow 2010. +// Copyright Paul A. Bristow 2010, 2014. // Use, modification and distribution are subject to 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) -// Example of finding roots using Newton-Raphson, Halley +// Example of finding roots using Newton-Raphson, Halley, Schroeder, TOMS748 . // Note that this file contains Quickbook mark-up as well as code // and comments, don't change any of the special comment mark-ups! -//#ifdef _MSC_VER -//# pragma warning(disable: 4180) // qualifier has no effect (in Fusion). -//#endif - +// To get (copious) diagnostic output, add make this define here or elsewhere. //#define BOOST_MATH_INSTRUMENT -//[root_finding_example1 -/*` +//[root_finding_headers +/* This example demonstrates how to use the various tools for root finding taking the simple cube root function (cbrt) as an example. It shows how use of derivatives can improve the speed. @@ -46,8 +43,11 @@ using boost::math::tools::toms748_solve; // using boost::math::tie; // which provide convenient aliases for various implementations, // including std::tr1, depending on what is available. +#include // pair, make_pair -//] [/root_finding_example1] +#include + +//] [/root_finding_headers] #include using std::cout; using std::endl; @@ -56,8 +56,7 @@ using std::setw; using std::setprecision; #include using std::numeric_limits; -//[root_finding_example2 -/*` +/* Let's suppose we want to find the cube root of a number. @@ -79,15 +78,17 @@ __spaces ['f]\'(x) = 2x[sup2] __spaces ['f]\'\'(x) = 6x */ +//[root_finding_cbrt_functor_noderiv + template -struct cbrt_functor_1 -{ // cube root of x using only function - no derivatives. - cbrt_functor_1(T const& to_find_root_of) : value(to_find_root_of) +struct cbrt_functor_noderiv +{ // Cube root of x using only function - no derivatives. + cbrt_functor_noderiv(T const& to_find_root_of) : value(to_find_root_of) { // Constructor stores value to find root of. - // For example: calling cbrt_functor_ (x) to get cube root of x. + // For example: calling cbrt_functor (x) to get cube root of x. } T operator()(T const& x) - { // Return both f(x)only. + { //! \returns f(x) - value. T fx = x*x*x - value; // Difference (estimate x^3 - value). return fx; } @@ -95,6 +96,12 @@ private: T value; // to be 'cube_rooted'. }; +//] [/root_finding_cbrt_functor_noderiv] + +//cout << ", std::numeric_limits<" << typeid(T).name() << ">::digits = " << digits +// << ", accuracy " << get_digits << " bits."<< endl; + + /*`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. @@ -104,46 +111,61 @@ 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'). */ +//[root_finding_cbrt_noderiv + template -T cbrt_1(T x) -{ // return cube root of x using bracket_and_solve (no derivatives). +T cbrt_noderiv(T x) +{ //! \returns cube root of x using bracket_and_solve (no derivatives). using namespace std; // Help ADL of std functions. - using namespace boost::math; + using namespace boost::math; // 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; // To multiply - int digits = std::numeric_limits ::digits; // Maximum possible binary digits accuracy for type T. + T factor = 2; // To multiply and divide guess to bracket. // digits used to control how accurate to try to make the result. - int get_digits = (digits * 3) /4; // Near maximum (3/4) possible accuracy. - //cout << ", std::numeric_limits<" << typeid(T).name() << ">::digits = " << digits - // << ", accuracy " << get_digits << " bits."<< endl; + // int digits = 3 * std::numeric_limits ::digits / 4; // 3/4 maximum possible binary digits accuracy for type T. + int digits = std::numeric_limits ::digits; // Maximum possible binary digits accuracy for type T. //boost::uintmax_t maxit = (std::numeric_limits (); bool is_rising = true; // So if result if guess^3 is too low, try increasing guess. - eps_tolerance tol(get_digits); + eps_tolerance tol(digits); std::pair (x), guess, factor, is_rising, tol, it); + bracket_and_solve_root(cbrt_functor_noderiv (x), guess, factor, is_rising, tol, it); - // Can show how many iterations (this information is lost outside cbrt_1). + // Can show how many iterations (this information is lost outside cbrt_noderiv). cout << "Iterations " << maxit << endl; if(it >= maxit) { // cout << "Unable to locate solution in chosen iterations:" " Current best guess is between " << r.first << " and " << r.second << endl; } - return r.first + (r.second - r.first)/2; // Midway between brackets. -} // T cbrt_1(T x) + T distance = float_distance(r.first, r.second); + std::cout << distance << " bits separate brackets." << std::endl; + for (int i = 0; i < distance; i++) + { + std::cout << float_advance(r.first, i) << std::endl; + } + + return r.first + (r.second - r.first) / 2; // Midway between bracketed interval. +} // T cbrt_noderiv(T x) + +//] [/root_finding_cbrt_noderiv] + + + +// maxit = 10 +// Unable to locate solution in chosen iterations: Current best guess is between 3.0365889718756613 and 3.0365889718756627 -//[root_finding_example2 /*` 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. @@ -164,47 +186,64 @@ both the evaluation of the function to solve, along with its first derivative: To \'return\' two values, we use a pair of floating-point values: */ + +//[root_finding_cbrt_functor_1stderiv + template -struct cbrt_functor_2 -{ // Functor also returning 1st derviative. - cbrt_functor_2(T const& to_find_root_of) : value(to_find_root_of) - { // Constructor stores value to find root of, - // for example: calling cbrt_functor_2 (x) to use to get cube root of x. +struct cbrt_functor_1stderiv +{ // Functor returning function and 1st derivative. + + cbrt_functor_1stderiv(T const& target) : value(target) + { // Constructor stores the value to be 'cube_rooted'. } - std::pair -T cbrt_2(T x) -{ // return cube root of x using 1st derivative and Newton_Raphson. +T cbrt_1deriv(T x) +{ //! \return cube root of x using 1st derivative and Newton_Raphson. + using namespace std; // For frexp, ldexp, numeric_limits. + using namespace boost::math::tools; // For newton_raphson_iterate. + int exponent; - frexp(x, &exponent); // Get exponent of z (ignore mantissa). + frexp(x, &exponent); // Get exponent of x (ignore mantissa). T guess = ldexp(1., exponent/3); // Rough guess is to divide the exponent by three. + // Set an initial bracket interval. 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. - int digits = std::numeric_limits ::digits; // Maximum possible binary digits accuracy for type T. // digits used to control how accurate to try to make the result. - int get_digits = (digits * 3) /4; // Near maximum (3/4) possible accuracy. + int digits = std::numeric_limits ::digits; // Maximum possible binary digits accuracy for type T. + + boost::uintmax_t maxit = 20; // Optionally limit the number of iterations. + //cout << "Max Iterations " << maxit << endl; // + T result = newton_raphson_iterate(cbrt_functor_1stderiv (x), guess, min, max, digits, maxit); + // Can check and show how many iterations (updated by newton_raphson_iterate). + // cout << "Iterations " << maxit << endl; + return result; +} // cbrt_1deriv + +//] [/root_finding_cbrt_1deriv] + +// int get_digits = (digits * 2) /3; // Two thirds of maximum possible accuracy. //boost::uintmax_t maxit = (std::numeric_limits (x), guess, min, max, get_digits, maxit); - // Can show how many iterations (updated by newton_raphson_iterate) but lost on exit. - // cout << "Iterations " << maxit << endl; - return result; -} /*` Finally need to define yet another functor that returns @@ -216,12 +255,13 @@ f''(x) = 3 * 3x To \'return\' three values, we use a tuple of three floating-point values: */ +//[root_finding_cbrt_functor_2deriv + template -struct cbrt_functor_3 +struct cbrt_functor_2deriv { // Functor returning both 1st and 2nd derivatives. - cbrt_functor_3(T const& to_find_root_of) : value(to_find_root_of) + cbrt_functor_2deriv(T const& to_find_root_of) : value(to_find_root_of) { // Constructor stores value to find root of, for example: - // calling cbrt_functor_3 (x) to get cube root of x, } // using boost::math::tuple; // to return three values. @@ -235,72 +275,73 @@ struct cbrt_functor_3 } private: T value; // to be 'cube_rooted'. -}; // struct cbrt_functor_3 +}; // struct cbrt_functor_2deriv + +//] [/root_finding_cbrt_functor_2deriv] + /*`Our cube function is now:*/ +//[root_finding_cbrt_2deriv + template -T cbrt_3(T x) +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; + using namespace std; // Help ADL of std functions. + using namespace boost::math; // halley_iterate + 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. - int digits = std::numeric_limits ::digits; // Maximum possible binary digits accuracy for type T. - // digits used to control how accurate to try to make the result. - int get_digits = (digits * 3) /4; // Near maximum (3/4) possible accuracy. - //cout << "Value " << x << ", guess " << guess - // << ", min " << min << ", max " << max - // << ", std::numeric_limits<" << typeid(T).name() << ">::digits = " << digits - // << ", accuracy " << get_digits << " bits."<< endl; + T max = ldexp(2., exponent/3); // Maximum possible value is twice our guess. - //boost::uintmax_t maxit = (std::numeric_limits ::digits /2 ; // Half maximum possible binary digits accuracy for type T. boost::uintmax_t maxit = 10; - //cout << "Max Iterations " << maxit << endl; // - T result = halley_iterate(cbrt_functor_3 (x), guess, min, max, digits, maxit); - // Can show how many iterations (updated by newton_raphson_iterate). - cout << "Iterations " << maxit << endl; + T result = halley_iterate(cbrt_functor_2deriv (x), guess, min, max, digits, maxit); + // Can show how many iterations (updated by halley_iterate). + // cout << "Iterations " << maxit << endl; return result; -} // cbrt_3(x) +} // cbrt_2deriv(x) +//] [/root_finding_cbrt_2deriv] + + // Default is + //boost::uintmax_t maxit = (std::numeric_limits ::max_digits10); - // Show all possibly significant decimal digits. - try - { - double v27 = 27; // that has an exact integer cube root. - double v28 = 28; // whose cube root is not exactly representable. +//[root_finding_example_1 + cout << "Cube Root finding (cbrt) Example." << endl; + // Show all possibly significant decimal digits. + cout.precision(std::numeric_limits ::max_digits10); + // or use cout.precision(max_digits10 = 2 + std::numeric_limits ::digits * 3010/10000); + try + { // Always use try'n'catch blocks with Boost.Math to get any error messages. + + double v27 = 27; // Example of a value that has an exact integer cube root. + double v28 = 28; // Example of a value whose cube root is *not* exactly representable. // Using bracketing: - double r = cbrt_1(v27); - cout << "cbrt_1(" << v27 << ") = " << r << endl; - r = cbrt_1(v28); - cout << "cbrt_1(" << v28 << ") = " << r << endl; + double r = cbrt_noderiv(v27); + cout << "cbrt_noderiv(" << v27 << ") = " << r << endl; + r = cbrt_noderiv(v28); + cout << "cbrt_noderiv(" << v28 << ") = " << r << endl; // Using 1st differential Newton-Raphson: - r = cbrt_2(v27); - cout << "cbrt_1(" << v27 << ") = " << r << endl; - r = cbrt_2(v28); - cout << "cbrt_2(" << v28 << ") = " << r << endl; + r = cbrt_1deriv(v27); + cout << "cbrt_1deriv(" << v27 << ") = " << r << endl; + r = cbrt_1deriv(v28); + cout << "cbrt_1deriv(" << v28 << ") = " << r << endl; // Using Halley with 1st and 2nd differentials. - r = cbrt_3(v27); - cout << "cbrt_3(" << v27 << ") = " << r << endl; - r = cbrt_3(v28); - cout << "cbrt_3(" << v28 << ") = " << r << endl; - - - - //] [/root_finding_example2] + r = cbrt_2deriv(v27); + cout << "cbrt_2deriv(" << v27 << ") = " << r << endl; + r = cbrt_2deriv(v28); + cout << "cbrt_2deriv(" << v28 << ") = " << r << endl; } catch(const std::exception& e) { // Always useful to include try & catch blocks because default policies @@ -310,6 +351,7 @@ int main() std::cout << "\n""Message from thrown exception was:\n " << e.what() << std::endl; } +//] [/root_finding_example_1 return 0; } // int main() @@ -318,24 +360,18 @@ int main() Normal output is: [pre - root_finding_example.cpp - Generating code - Finished generating code - root_finding_example.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Release\root_finding_example.exe - Cube Root finding (cbrt) Example. + Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Release\root_finding_example.exe"1> Cube Root finding (cbrt) Example. Iterations 10 - cbrt_1(27) = 3 + cbrt_noderiv(27) = 3 Iterations 10 Unable to locate solution in chosen iterations: Current best guess is between 3.0365889718756613 and 3.0365889718756627 - cbrt_1(28) = 3.0365889718756618 - cbrt_1(27) = 3 - cbrt_2(28) = 3.0365889718756627 - Iterations 4 - cbrt_3(27) = 3 - Iterations 5 - cbrt_3(28) = 3.0365889718756627 - -] [/pre] + cbrt_noderiv(28) = 3.0365889718756618 + cbrt_1deriv(27) = 3 + cbrt_1deriv(28) = 3.0365889718756627 + cbrt_2deriv(27) = 3 + cbrt_2deriv(28) = 3.0365889718756627 +] +[/pre] to get some (much!) diagnostic output we can add
- These functions all perform iterative root finding using derivatives:
+ These functions newton_raphson_iterate, halley_iterate,
+ and schroeder_iterate all
+ perform iterative root finding using derivatives.
+
+ (These functions should be used whenever derivative(s) are available, otherwise + use root-finding without derivatives).