diff --git a/doc/background/implementation.qbk b/doc/background/implementation.qbk index e043d5e8d..e67c1775e 100644 --- a/doc/background/implementation.qbk +++ b/doc/background/implementation.qbk @@ -23,10 +23,10 @@ has received considerable effort. (It also required much CPU effort - there was some danger of molten plastic dripping from the bottom of JM's laptop, -so instead, PAB's Dual-core desktop was kept 50% busy for *days* +so instead, PAB's Dual-core desktop was kept 50% busy for [*days] calculating some tables of test values!) -For a specific RealType, say float or double, +For a specific RealType, say `float` or `double`, it may be possible to find approximations for some functions that are simpler and thus faster, but less accurate (perhaps because there are no refining iterations, @@ -55,7 +55,7 @@ constant values are given to 50 decimal digits if available Values are specified as long double types by appending L, unless they are exactly representable, for example integers, or binary fractions like 0.125. This avoids the risk of loss of accuracy converting from double, the default type. -Values are used after static_cast(1.2345L) +Values are used after `static_cast(1.2345L)` to provide the appropriate RealType for spot tests. Functions that return constants values, like kurtosis for example, are written as diff --git a/doc/background/special_tut.qbk b/doc/background/special_tut.qbk index 71036ebe0..4a1b5fdef 100644 --- a/doc/background/special_tut.qbk +++ b/doc/background/special_tut.qbk @@ -2,15 +2,15 @@ [section:special_tut_impl Implementation] -In this section we'll provide a "recipe" for adding a new special function to this library to make life easier for -future authors wishing to contribute. We'll assume the function returns a single floating point result, and takes -two floating point arguments. For the sake of exposistion we'll give the function the name "my_special". +In this section, we'll provide a "recipe" for adding a new special function to this library to make life easier for +future authors wishing to contribute. We'll assume the function returns a single floating-point result, and takes +two floating-point arguments. For the sake of exposition we'll give the function the name [~my_special]. -Normally the implementation of such a function is split into two layers - a public user layer, and an internal -implementation layer that does the actual work. The implementation layer is declared inside a "detail" namespace -and has a simple signature: +Normally, the implementation of such a function is split into two layers - a public user layer, and an internal +implementation layer that does the actual work. +The implementation layer is declared inside a `detail` namespace and has a simple signature: - namespace boost{ namespace math{ namespace detail{ + namespace boost { namespace math { namespace detail { template T my_special_imp(const T& a, const T&b, const Policy& pol) @@ -109,7 +109,7 @@ what each line does: We're now almost there, we just need to flesh out the details of the implementation layer: - namespace boost{ namespace math{ namespace detail{ + namespace boost { namespace math { namespace detail { template T my_special_imp(const T& a, const T&b, const Policy& pol) @@ -125,13 +125,13 @@ The following guidelines indicate what (other than basic arithmetic) can go in t [link math_toolkit.error_handling.finding_more_information policy based error handlers]. * Calls to standard library functions should be made unqualified (this allows argument dependent lookup to find standard library functions for user-defined floating point -types such as those from Boost.Multiprecision). In addition the macro `BOOST_MATH_STD_USING` +types such as those from __multiprecision). In addition, the macro `BOOST_MATH_STD_USING` should appear at the start of the function (note no semi-colon afterwards!) so that all the math functions in `namespace std` are visible in the current scope. * Calls to other special functions should be made as fully qualified calls, and include the policy parameter as the last argument, for example `boost::math::tgamma(a, pol)`. * Where possible, evaluation of series, continued fractions, polynomials, or root -finding should use one of the [link toolkit boiler plate functions]. In any case, after +finding should use one of the [link math_toolkit.internals_overview 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 did terminate as a result of exceeding the maximum, then the appropriate error handler should be called (see existing code for examples). @@ -139,7 +139,7 @@ maximum, then the appropriate error handler should be called (see existing code for example: `constants::pi()`. * Where tables of coefficients are used (for example for rational approximations), care should be taken to ensure these are initialized at program startup to ensure thread safety when using user-defined number types. -See for example the use of `erf_initializer` in boost/math/special_functions/erf.hpp. +See for example the use of `erf_initializer` in [@../include/boost/math/special_functions/erf.hpp erf.hpp]. Here are some other useful internal functions: diff --git a/doc/distributions/skew_normal.qbk b/doc/distributions/skew_normal.qbk index 9dde567c5..ca1267556 100644 --- a/doc/distributions/skew_normal.qbk +++ b/doc/distributions/skew_normal.qbk @@ -2,14 +2,14 @@ ``#include `` - namespace boost{ namespace math{ - - template class skew_normal_distribution; - + typedef skew_normal_distribution<> normal; - + template class skew_normal_distribution { @@ -24,17 +24,17 @@ RealType shape()const; // The distribution is right skewed if shape > 0 and is left skewed if shape < 0. // The distribution is normal if shape is zero. }; - + }} // namespaces - + The skew normal distribution is a variant of the most well known -Gaussian statistical distribution. +Gaussian statistical distribution. The skew normal distribution with shape zero resembles the [@http://en.wikipedia.org/wiki/Normal_distribution Normal Distribution], hence the latter can be regarded as a special case of the more generic skew normal distribution. -If the standard (mean = 0, scale = 1) normal distribution probability density function is +If the standard (mean = 0, scale = 1) normal distribution probability density function is [space][space][equation normal01_pdf] @@ -51,7 +51,7 @@ with shape parameter [alpha], defined by O'Hagan and Leonhard (1976) is Given [@http://en.wikipedia.org/wiki/Location_parameter location] [xi], [@http://en.wikipedia.org/wiki/Scale_parameter scale] [omega], and [@http://en.wikipedia.org/wiki/Shape_parameter shape] [alpha], -it can be +it can be [@http://en.wikipedia.org/wiki/Skew_normal_distribution transformed], to the form: @@ -72,24 +72,24 @@ in the following graphs: [h4 Member Functions] skew_normal_distribution(RealType location = 0, RealType scale = 1, RealType shape = 0); - + Constructs a skew_normal distribution with location [xi], scale [omega] and shape [alpha]. Requires scale > 0, otherwise __domain_error is called. - RealType location()const; - + RealType location()const; + returns the location [xi] of this distribution, - + RealType scale()const; returns the scale [omega] of this distribution, - + RealType shape()const; returns the shape [alpha] of this distribution. - + (Location and scale function match other similar distributions, allowing the functions `find_location` and `find_scale` to be used generically). @@ -137,7 +137,7 @@ and at [@http://cran.r-project.org/web/packages/sn/sn.pd R skew-normal(sn) packa Package sn provides functions related to the skew-normal (SN) and the skew-t (ST) probability distributions, both for the univariate and for the the multivariate case, -including regression models. +including regression models. __Mathematica was also used to generate some more accurate spot test data. @@ -145,7 +145,7 @@ __Mathematica was also used to generate some more accurate spot test data. The skew_normal distribution with shape = zero is implemented as a special case, equivalent to the normal distribution in terms of the -[link math_toolkit.sf_erf.error_function error function], +[link math_toolkit.sf_erf.error_function error function], and therefore should have excellent accuracy. The PDF and mean, variance, skewness and kurtosis are also accurately evaluated using @@ -165,7 +165,7 @@ and [omega] is its scale, and [alpha] is its shape. [table [[Function][Implementation Notes]] [[pdf][Using:[equation skew_normal_pdf] ]] -[[cdf][Using: [equation skew_normal_cdf]\n +[[cdf][Using: [equation skew_normal_cdf][br] where ['T(h,a)] is Owen's T function, and ['[Phi](x)] is the normal distribution. ]] [[cdf complement][Using: complement of normal distribution + 2 * Owens_t]] [[quantile][Maximum of the pdf is sought through searching the root of f'(x)=0]] diff --git a/doc/math.qbk b/doc/math.qbk index 56c84d097..b82948308 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -102,11 +102,13 @@ and use the function's name as the link text.] [def __bracket_solve [link math_toolkit.roots.roots_noderiv.bracket_solve bracket and solve]] [def __root_finding_examples [link math_toolkit.roots.root_finding_examples root-finding examples]] [def __brent_minima [link math_toolkit.roots.brent_minima Brent's method]] +[def __brent_minima_example [link math_toolkit.roots.brent_minima.example Brent's method example]] [def __newton[link math_toolkit.roots.roots_deriv.newton Newton-Raphson iteration]], [def __halley [link math_toolkit.roots.roots_deriv.halley Halley]] [def __schroeder [link math_toolkit.roots.roots_deriv.schroeder Schroeder]] [def __brent_minima_example [link math_toolkit.roots.brent_minima.example Brent minima finding example]] -[def __bisection [link math_toolkit.roots.roots_noderiv.bisection Bisection]] +[def __bisection [link math_toolkit.roots.roots_noderiv.bisect bisection]] +[def __bisect [link math_toolkit.roots.roots_noderiv.bisect bisect]] [def __bisection_wikipedia [@https://en.wikipedia.org/wiki/Bisection bisection]] [def __root_finding_example_without_derivatives [link math_toolkit.roots.root_finding_examples.no_derivatives root-finding without derivatives]] [def __root_finding_example_cbrt_without_derivatives [link math_toolkit.roots.root_finding_examples.cbrt_no_derivatives cube root-finding without derivatives]] @@ -362,7 +364,8 @@ and use the function's name as the link text.] [def __MPFR [@http://www.mpfr.org/ GNU MPFR library]] [def __GMP [@http://gmplib.org/ GNU Multiple Precision Arithmetic Library]] [def __multiprecision [@boost:/libs/multiprecision/doc/html/index.html Boost.Multiprecision]] -[def __cpp_dec_float [@http://www.boost.org/doc/libs/1_53_0_beta1/libs/multiprecision/doc/html/boost_multiprecision/tut/floats/cpp_dec_float.html cpp_dec_float]] +[def __cpp_dec_float [@boost:/libs/multiprecision/doc/html/boost_multiprecision/tut/floats/cpp_dec_float.html cpp_dec_float]] +[def __cpp_bin_float [@boost:/libs/multiprecision/doc/html/boost_multiprecision/tut/floats/cpp_bin_float.html cpp_bin_float]] [def __R [@http://www.r-project.org/ The R Project for Statistical Computing]] [def __godfrey [link godfrey Godfrey]] [def __pugh [link pugh Pugh]] diff --git a/doc/overview/overview.qbk b/doc/overview/overview.qbk index 16e46a776..47f437be5 100644 --- a/doc/overview/overview.qbk +++ b/doc/overview/overview.qbk @@ -49,24 +49,28 @@ sizes: typically `float`, `double` or `long double`. [h4 Implementation Toolkit] -Provides [link toolkit many of the tools] required to implement +The section [link math_toolkit.internals_overview Internal tools] +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 details" -whose interfaces and\/or implementations may change. +function implementations in the future. -There are helpers for the -[link math_toolkit.internals.series_evaluation evaluation of infinite series], -[link math_toolkit.internals.cf continued -fractions] and [link math_toolkit.internals.rational -rational approximations]. +Some tools are now considered well-tried and their signatures stable and unlikely to change. There is a fairly comprehensive set of root finding both __root_finding_without_derivatives and __root_finding_with_derivatives with derivative support, and function minimization using __brent_minima. +Other [link math_toolkit.internals_overview Internal tools] +are currently still considered experimental: they are "exposed implementation details" +whose interfaces and\/or implementations may change without notice. + +There are helpers for the +[link math_toolkit.internals.series_evaluation evaluation of infinite series], +[link math_toolkit.internals.cf continued +fractions] and [link math_toolkit.internals.rational +rational approximations]. A [link math_toolkit.internals.minimax Remez algorithm implementation] allows for the locating of minimax rational approximations. @@ -82,7 +86,7 @@ external graphing application. [endsect] [/section:intro Introduction] [/ - Copyright 2006, 2012 John Maddock and Paul A. Bristow. + Copyright 2006, 2012, 2015 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/roots/minima.qbk b/doc/roots/minima.qbk index 4ccdc6b66..0d6e47647 100644 --- a/doc/roots/minima.qbk +++ b/doc/roots/minima.qbk @@ -58,7 +58,7 @@ of ['f(x)] at the minima. [h4:example Brent Minimisation Example] -As a demonstration, we replicate this [@http://en.wikipedia.org/wiki/Brent%27s_method#Example example] +As a demonstration, we replicate this [@http://en.wikipedia.org/wiki/Brent%27s_method#Example Wikipedia example] minimising the function ['y= (x+3)(x-1)[super 2]]. It is obvious from the equation and the plot that there is a diff --git a/doc/roots/root_finding_examples.qbk b/doc/roots/root_finding_examples.qbk index f4d2963fd..11d07004a 100644 --- a/doc/roots/root_finding_examples.qbk +++ b/doc/roots/root_finding_examples.qbk @@ -7,13 +7,14 @@ The examples demonstrate how to use the various tools for [@http://en.wikipedia.org/wiki/Root-finding_algorithm root finding]. -We start with the simple cube root function [cbrt] ( C++ standard function name +We start with the simple cube root function `cbrt` ( C++ standard function name [@http://en.cppreference.com/w/cpp/numeric/math/cbrt cbrt]) showing __cbrt_no_derivatives. We then show how use of derivatives can improve the speed of convergence. -(But these examples are only a demonstration and does not try to make the ultimate improvements of a 'real-life' +(But these examples are only a demonstration and do not try to make +the ultimate improvements of an 'industrial-strength' implementation, for example, of `boost::math::cbrt`, mainly by using a better computed initial 'guess' at [@boost:/libs/math/include/boost/math/special_functions/cbrt.hpp cbrt.hpp]). @@ -47,10 +48,10 @@ 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, -the 1st derivation ['f'(x)], will speed homing in on the solution. +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. +Lastly, we show how adding the curvature ['f''(x)] too will speed convergence even more. [h3:cbrt_no_derivatives Cube root function without derivatives] @@ -193,11 +194,11 @@ Or to put it another way: ['nothing beats a really good initial guess!] [h3:fifth_root Fifth-root function] -Let's now suppose we want to find the [*fifth root] of a number a. +Let's now suppose we want to find the [*fifth root] of a number ['a]. The equation we want to solve is : -__spaces['f](x) = x[super 5] -a +__spaces['f](x) = ['x[super 5] -a] If your differentiation is a little rusty (or you are faced with an function whose complexity makes differentiation daunting), @@ -240,7 +241,7 @@ __spaces['f]\'\'(x) = 20x[super 3] The apocryphally astute reader might, by now, be asking "How do we know if this computes the 'right' answer?". For most values, there is, sadly, no 'right' answer. -This is because values can only rarely be exactly represented by C++ floating-point types. +This is because values can only rarely be ['exactly represented] by C++ floating-point types. What we do want is the 'rightest' representation - one that is the nearest representable value. (For more about how numbers are represented see __floating_point). @@ -256,20 +257,20 @@ to be the [*nearest representable] `double` value. For example, the cube root functions in our example for cbrt(28.) return - std::cbrt(28) 3.0365889718756627 +`std::cbrt(28.) = 3.0365889718756627` - WolframAlpha says 3.036588971875662519420809578505669635581453977248111123242141... +WolframAlpha says `3.036588971875662519420809578505669635581453977248111123242141...` - static_cast(3.03658897187566251942080957850) - 3.0365889718756627 +`static_cast(3.03658897187566251942080957850) = 3.0365889718756627` - example cbrt(28) 3.0365889718756627 +This example `cbrt(28.) = 3.0365889718756627` [tip To ensure that all potentially significant decimal digits are displayed use `std::numeric_limits::max_digits10` -or if not available on older platforms or compilers use `2 + std::numeric_limits::digits * 3010/10000`. +(or if not available on older platforms or compilers use `2+std::numeric_limits::digits*3010/10000`).[br] + Ideally, values should agree to `std::numeric-limits::digits10` decimal digits. -This also means that a 'reference' value to be [*input ]or `static_cast` should have +This also means that a 'reference' value to be [*input] or `static_cast` should have at least `max_digits10` decimal digits (17 for 64-bit `double`). ] @@ -281,17 +282,17 @@ Almost all platforms can easily use __multiprecision, for example, __cpp_dec_float or a binary type __cpp_bin_float types, to compute values at very much higher precision. -[note With multiprecision types, it is debatable whether to use the type ['T] for computing the initial guesses. +[note With multiprecision types, it is debatable whether to use the type `T` for computing the initial guesses. Type `double` is like to be accurate enough for the method used in these examples. This would limit the range of possible values to that of `double`. -There is also the cost of conversion to and from type T to consider. +There is also the cost of conversion to and from type `T` to consider. In these examples, `double` is used via `typedef double guess_type`.] Since the functors and functions used above are templated on the value type, we can very simply use them with any of the __multiprecision types. Some examples below are 50 decimal digit decimal and binary types -(and on some platforms a much faster `float128 type` or `quad_float`) +(and on some platforms a much faster `float128` or `quad_float` type ) that we can use with these includes: [root_finding_multiprecision_include_1] @@ -326,7 +327,7 @@ value = 2, cube root =1.25992104989487 value = 2, cube root =1.2599210498948731647672106072782283505702514647015 ] -[tip Be [*very careful] about the floating-point type ['T] that is passed to the root-finding function. +[tip Be [*very careful] about the floating-point type `T` that is passed to the root-finding function. Carelessly passing a integer by writing `cpp_dec_float_50 r = cbrt_2deriv(2);` or `show_cube_root(2);` will provoke many warnings and compile errors. @@ -337,13 +338,13 @@ used to compute the guess and bracket values as `double`. Even more treacherous is passing a `double` as in `cpp_dec_float_50 r = cbrt_2deriv(2.);` which silently gives the 'wrong' result, computing a `double` result and [*then] converting to `cpp_dec_float_50`! All digits beyond `max_digits10` will be incorrect. -Making the `cbrt` type explicit with `cbrt_2deriv(2.);` will give you the desired result. +Making the `cbrt` type explicit with `cbrt_2deriv(2.);` will give you the desired 5o decimal digit precision result. ] [/tip] [h3:nth_root Generalizing to Compute the nth root] If desired, we can now further generalize to compute the ['n]th root by computing the derivatives [*at compile-time] -using the rules for differentiation and the compile-time static function `pow` +using the rules for differentiation and the compile-time `static` function `pow` where template parameter `N` is an integer. Our functor and function now have an additional template parameter `N`, for the root required. @@ -374,7 +375,7 @@ as noted above.] [warning Asking for unreasonable roots, for example, `show_nth_root<1000000>(2.);` may lead to [@http://en.wikipedia.org/wiki/Loss_of_significance Loss of significance] like `Type double value = 2, 1000000th root = 1.00000069314783`. -Use of the the `pow` function is more sensible for this need. +Use of the the `pow` function is more sensible for this unusual need. ] Full code of these examples is at diff --git a/doc/roots/roots.qbk b/doc/roots/roots.qbk index c6ed0d795..4cc14bfcb 100644 --- a/doc/roots/roots.qbk +++ b/doc/roots/roots.qbk @@ -6,36 +6,38 @@ #include `` - namespace boost{ namespace math{ - namespace tools{ - + namespace boost { namespace math { + namespace tools { // Note namespace boost::math::tools. + // Newton-Raphson template T newton_raphson_iterate(F f, T guess, T min, T max, int digits); template T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter); + // Halley template T halley_iterate(F f, T guess, T min, T max, int digits); template T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter); + // Schroeder template T schroeder_iterate(F f, T guess, T min, T max, int digits); template T schroeder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter); - }}} // namespaces + }}} // namespaces boost::math::tools. [h4 Description] These functions all perform iterative root-finding [*using derivatives]: -* `newton_raphson_iterate` performs second order __newton. +* `newton_raphson_iterate` performs second-order __newton. -* `halley_iterate` and`schroeder_iterate` perform third order +* `halley_iterate` and`schroeder_iterate` perform third-order __halley and __schroeder iteration. The functions all take the same parameters: @@ -44,11 +46,11 @@ The functions all take the same parameters: [[F f] [Type F must be a callable function object that accepts one parameter and returns a __tuple: -For the second order iterative methods ([@http://en.wikipedia.org/wiki/Newton_Raphson Newton Raphson]) +For second-order iterative method ([@http://en.wikipedia.org/wiki/Newton_Raphson Newton Raphson]) the __tuple should have [*two] elements containing the evaluation of the function and its first derivative. -For the third order methods +For the third-order methods ([@http://en.wikipedia.org/wiki/Halley%27s_method Halley] and Schroeder) the __tuple should have *three* elements containing the evaluation of @@ -58,7 +60,7 @@ Schroeder) [[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 precision.]] [[uintmax_t& max_iter] [An optional maximum number of iterations to perform. On exit, this is updated to the actual number of iterations performed.]] -[[T] [Floating-point type, typically `double` (often deduced by __ADL)]] +[[T] [Floating-point type, typically `double`, often deduced by __ADL.]] [[F] [Floating-point types of the functor, usually deduced from T by __ADL, so it is rarely necessary to define F.]] ] @@ -67,10 +69,10 @@ When using these functions you should note that: * Default `max_iter = (std::numeric_limits::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 +can be no better than __bisect, or in some rare cases, even worse than __bisect 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 +derivatives, and fall back to __bisect 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 @@ -82,7 +84,7 @@ the bracket bounds ['min] and ['max] may as well be set to the widest limits like zero and `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 +* These functions fall back to __bisect 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. @@ -100,7 +102,9 @@ iteration just takes the next value into the zone where ['f(x)] becomes inaccura * 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 +and also ensure that display of all the significant digits with +` cout.precision(std::numeric_limits::digits10)`: +or even possibly significant digits with ` cout.precision(std::numeric_limits::max_digits10)`: but be warned, this may produce copious output! * Finally: you may well be able to do better than these functions by hand-coding @@ -114,7 +118,7 @@ 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. diff --git a/doc/roots/roots_overview.qbk b/doc/roots/roots_overview.qbk index 563e7a967..803388335 100644 --- a/doc/roots/roots_overview.qbk +++ b/doc/roots/roots_overview.qbk @@ -2,14 +2,14 @@ Several tools are provided to aid finding minima and roots of functions. Some __root_finding_without_derivatives methods are __bisection, -__bracket_solve including use of __root_finding_TOMS748. +__bracket_solve, including use of __root_finding_TOMS748. For __root_finding_with_derivatives the methods of -__newton __halley and __schroeder are implemented. +__newton, __halley, and __schroeder are implemented. -For locating minima of a function, __brent_minima is provided. +For locating minima of a function, a __brent_minima_example is provided. -Also provided are __root_finding_examples: +There are several fully-worked __root_finding_examples, including: * __root_finding_example_cbrt_without_derivatives * __root_finding_example_cbrt_with_1_derivative diff --git a/doc/roots/roots_without_derivatives.qbk b/doc/roots/roots_without_derivatives.qbk index 2d7ad7529..6b9ca9894 100644 --- a/doc/roots/roots_without_derivatives.qbk +++ b/doc/roots/roots_without_derivatives.qbk @@ -6,8 +6,8 @@ #include `` - namespace boost{ namespace math{ - namespace tools{ + namespace boost { namespace math { + namespace tools { // Note namespace boost::math::tools. // Bisection template std::pair @@ -113,15 +113,15 @@ [h4 Description] -These functions solve the root of some function ['f(x)] ['without the -need for any derivatives of ['f(x)]]. +These functions solve the root of some function ['f(x)] +['without the need for any derivatives of ['f(x)]]. The `bracket_and_solve_root` functions use __root_finding_TOMS748 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 __policy_section are provided. -Alternatively, __bisection is a simple __bisection_wikipedia routine which can be useful +Alternatively, __bisect is a simple __bisection_wikipedia 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. @@ -146,28 +146,28 @@ 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. -[h6:bisection Bisection] +[h6:bisect Bisection] template std::pair - bisect( + bisect( // Unlimited iterations. + F f, + T min, + T max, + Tol tol); + + template + std::pair + bisect( // Limited iterations. F f, T min, T max, Tol tol, boost::uintmax_t& max_iter); - template - std::pair - bisect( - F f, - T min, - T max, - Tol tol); - template std::pair - bisect( + bisect( // Specified policy. F f, T min, T max, @@ -177,18 +177,18 @@ only rarely, but may be useful in some specific circumstances. These functions locate the root using __bisection_wikipedia. -Function arguments are: +`bisect` function arguments are: [variablelist [[f] [A unary functor which is the function ['f(x)] whose root is to be found.]] [[min] [The left bracket of the interval known to contain the root.]] -[[max] [The right bracket of the interval known to contain the root. +[[max] [The right bracket of the interval known to contain the root.[br] 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. + The action taken is controlled by the evaluation error policy.[br] 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.]] + will return the current brackets enclosing the root when ['tol(min, max)] becomes true.]] [[max_iter][The maximum number of invocations of ['f(x)] to make while searching for the root. On exit, this is updated to the actual number of invocations performed.]] [[T] [Floating-point type, typically `double` (often deduced by __ADL)]] [[F] [Floating-point types of the functor, usually deduced from T by __ADL, so it is rarely necessary to define F.]] @@ -196,7 +196,7 @@ Function arguments are: [optional_policy] -Returns: a pair of values ['r] that bracket the root so that: +[*Returns]: a pair of values ['r] that bracket the root so that: f(r.first) * f(r.second) <= 0 @@ -211,9 +211,9 @@ or 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 -as a result of exceeding /max_iter/ function invocations (easily done by -checking the updated value of /max_iter/ when the function returns), rather than -because the termination condition /tol/ was satisfied. +as a result of exceeding ['max_iter] function invocations (easily done by +checking the updated value of ['max_iter] when the function returns), rather than +because the termination condition ['tol] was satisfied. [h6:bracket_solve Bracket and solve] @@ -261,13 +261,13 @@ particular it is known whether the root is occurs for positive or negative 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.]] -[[T] [Floating-point type, typically `double` (often deduced by __ADL)]] +[[T] [Floating-point type, typically `double`, often deduced by __ADL]] [[F] [Floating-point types of the functor, usually deduced from T by __ADL, so it is rarely necessary to define F.]] ] [optional_policy] -Returns: a pair of values /r/ that bracket the root so that: +[*Returns]: a pair of values ['r] that bracket the root so that: f(r.first) * f(r.second) <= 0 @@ -279,12 +279,12 @@ or max_iter >= m -where /m/ is the initial value of /max_iter/ passed to the function. +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 -as a result of exceeding /max_iter/ function invocations (easily done by -checking the value of /max_iter/ when the function returns), rather than -because the termination condition /tol/ was satisfied. +as a result of exceeding ['max_iter] function invocations (easily done by +checking the value of ['max_iter] when the function returns), rather than +because the termination condition ['tol] was satisfied. [h6:TOMS748 Algorithm TOMS 748: Alefeld, Potra and Shi: Enclosing zeros of continuous functions] @@ -332,32 +332,32 @@ because the termination condition /tol/ was satisfied. 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 __root_finding_TOMS748 parameters are: +['f(x)]. The two functions differ only by whether values for ['f(a)] and +['f(b)] are already available. The __root_finding_TOMS748 parameters are: [variablelist [[f] [A unary functor that is the function whose root is to be solved. - f(x) need not be uniformly increasing or decreasing on /x/ and + f(x) need not be uniformly increasing or decreasing on ['x] and may have multiple roots.]] [[a] [ The lower bound for the initial bracket of the root.]] [[b] [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/.]] -[[fa] [Optional: the value of /f(a)/.]] -[[fb] [Optional: the value of /f(b)/.]] + It is a precondition that ['a < b] and that ['a] and ['b] + bracket the root to find so that ['f(a) * f(b) < 0].]] +[[fa] [Optional: the value of ['f(a)].]] +[[fb] [Optional: the value of ['f(b)].]] [[tol] [A binary functor that determines the termination condition for the search - for the root. /tol/ is passed the current brackets at each step, + 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.]] [[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 + for the root. On exit, ['max_iter] is set to actual number of function invocations used.]] -[[T] [Floating-point type, typically `double` (often deduced by __ADL)]] +[[T] [Floating-point type, typically `double`, often deduced by __ADL]] [[F] [Floating-point types of the functor, usually deduced from T by __ADL, so it is rarely necessary to define F.]] ] [optional_policy] -`toms748_solve` returns: a pair of values /r/ that bracket the root so that: +`toms748_solve` returns: a pair of values ['r] that bracket the root so that: f(r.first) * f(r.second) <= 0 @@ -369,13 +369,13 @@ or max_iter >= m -where /m/ is the initial value of /max_iter/ passed to the function. +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 -as a result of exceeding /max_iter/ function invocations (easily done by -checking the updated value of /max_iter/ +as a result of exceeding ['max_iter] function invocations (easily done by +checking the updated value of ['max_iter] against its previous value passed as parameter), -rather than because the termination condition /tol/ was satisfied. +rather than because the termination condition ['tol] was satisfied. template struct eps_tolerance @@ -385,11 +385,11 @@ rather than because the termination condition /tol/ was satisfied. }; `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/ +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 +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 type T] is +required to ensure that we get back a bracketing interval, since this must clearly be at least 1 epsilon in size. struct equal_floor @@ -399,8 +399,8 @@ be at least 1 epsilon in size. }; 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 { @@ -409,8 +409,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 { @@ -429,7 +429,9 @@ See __root_finding_examples. [h4 Implementation] The implementation of the bisection algorithm is extremely straightforward -and not detailed here. __TOMS748 is described in detail in: +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, diff --git a/example/root_finding_example.cpp b/example/root_finding_example.cpp index 67655e20f..a8ddbc31d 100644 --- a/example/root_finding_example.cpp +++ b/example/root_finding_example.cpp @@ -1,369 +1,490 @@ // root_finding_example.cpp -// Copyright Paul A. Bristow 2010, 2014. +// Copyright Paul A. Bristow 2010, 2015 // 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, Schroeder, TOMS748 . +// Example of finding roots using Newton-Raphson, Halley. // Note that this file contains Quickbook mark-up as well as code // and comments, don't change any of the special comment mark-ups! -// To get (copious) diagnostic output, add make this define here or elsewhere. +//#ifdef _MSC_VER +//# pragma warning(disable: 4180) // qualifier has no effect (in Fusion). +//#endif + //#define BOOST_MATH_INSTRUMENT -//[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. +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. (But is only a demonstration and does not try to make the ultimate improvements of 'real-life' -implementation of boost::math::cbrt; mainly by using a better computed initial 'guess' -at `` ). +implementation of `boost::math::cbrt`, mainly by using a better computed initial 'guess' +at ``). -First some includes that will be needed. -Using statements are provided to list what functions are being used in this example: -you can of course qualify the names in other ways. +Then we show how a high roots (fifth) can be computed, +and in root_finding_n_example.cpp a generic method +for the ['n]th root that constructs the derivatives at compile-time, + +These methods should be applicable to other functions that can be differentiated easily. + +First some `#includes` that will be needed. + +[tip For clarity, `using` statements are provided to list what functions are being used in this example: +you can of course partly or fully qualify the names in other ways. +(For your application, you may wish to extract some parts into header files, +but you should never use `using` statements globally in header files).] */ +//[root_finding_include_1 + #include -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; +//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 -#include // pair, make_pair +#include // For float_distance. +#include // for tuple and make_tuple. +#include // For boost::math::cbrt. -//] [/root_finding_headers] +//] [/root_finding_include_1] + +// using boost::math::tuple; +// using boost::math::make_tuple; +// using boost::math::tie; +// which provide convenient aliases for various implementations, +// including std::tr1, depending on what is available. #include -using std::cout; using std::endl; +//using std::cout; using std::endl; #include -using std::setw; using std::setprecision; +//using std::setw; using std::setprecision; #include -using std::numeric_limits; +//using std::numeric_limits; /* -Let's suppose we want to find the cube root of a number. +Let's suppose we want to find the root of a number ['a], and to start, compute the cube root. -The equation we want to solve is: +So the equation we want to solve is: -__spaces ['f](x) = x[cubed] +__spaces ['f](x) = x[cubed] - a We will first solve this without using any information -about the slope or curvature of the cbrt function. +about the slope or curvature of the cube root function. -We then show how adding what we can know, for this function, about the slope, -the 1st derivation /f'(x)/, will speed homing in on the solution, -and then finally how adding the curvature /f''(x)/ as well will improve even more. +We then show how adding what we can know about this function, first just the slope, +the 1st derivation /f'(x)/, will speed homing in on the solution. -The 1st and 2nd derivatives of x[cubed] are: +Lastly we show how adding the curvature /f''(x)/ too will speed convergence even more. -__spaces ['f]\'(x) = 2x[sup2] - -__spaces ['f]\'\'(x) = 6x */ -//[root_finding_cbrt_functor_noderiv +//[root_finding_noderiv_1 template 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. +{ // 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) - { //! \returns f(x) - value. - T fx = x*x*x - value; // Difference (estimate x^3 - value). + { + T fx = x*x*x - a; // Difference (estimate x^3 - a). return fx; } private: - T value; // to be 'cube_rooted'. + T a; // to be 'cube_rooted'. }; +//] [/root_finding_noderiv_1 -//] [/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: +/* +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'.) -Cube root function is 'Really Well Behaved' in that it is monotonic +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 +//[root_finding_noderiv_2 template T cbrt_noderiv(T x) -{ //! \returns cube root of x using bracket_and_solve (no derivatives). +{ // return cube root of x using bracket_and_solve (no derivatives). using namespace std; // Help ADL of std functions. - using namespace boost::math; // For bracket_and_solve_root. + 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; // To multiply and divide guess to bracket. - // digits used to control how accurate to try to make the result. - // 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. + T factor = 2; // How big steps to take when searching. - //boost::uintmax_t maxit = (std::numeric_limits::max)(); - // (std::numeric_limits::max)() = 18446744073709551615 - // which is more than we might wish to wait for!!! - // so we can choose some reasonable estimate of how many iterations may be needed. - - const boost::uintmax_t maxit = 10; + const boost::uintmax_t maxit = 20; // Limit to maximum iterations. boost::uintmax_t it = maxit; // Initally our chosen max iterations, but updated with actual. - // We could also have used a maximum iterations provided by any policy: - // boost::uintmax_t max_it = policies::get_max_root_iterations(); - bool is_rising = true; // So if result if guess^3 is too low, try increasing guess. - eps_tolerance tol(digits); - std::pair r = + bool is_rising = true; // So if result if guess^3 is too low, then try increasing guess. + int digits = std::numeric_limits::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) /4; // Near maximum (3/4) possible accuracy. + eps_tolerance tol(get_digits); // Set the tolerance. + std::pair r = 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_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; - } - 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 - + return r.first + (r.second - r.first)/2; // Midway between brackets. +} /*` + +[note Default `boost::uintmax_t maxit = (std::numeric_limits::max)();` +and `(std::numeric_limits::max)() = 18446744073709551615` +which is more than anyone would wish to wait for! + +So it may be wise to chose some reasonable estimate of how many iterations may be needed, here only 20.] + +[note We could also have used a maximum iterations provided by any policy: +`boost::uintmax_t max_it = policies::get_max_root_iterations();`] + +[tip One can show how many iterations 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). +*/ +//] [/root_finding_noderiv_2] + + +// Cube root with 1st derivative (slope) + +/* 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 this function, the 1st differential (the slope of the tangent to a curve at any point) is known. +For the root function, the 1st differential (the slope of the tangent to a curve at any point) is known. -[@http://en.wikipedia.org/wiki/Derivative#Derivatives_of_elementary_functions derivatives] -gives some reminders. +If you need some reminders then +[@http://en.wikipedia.org/wiki/Derivative#Derivatives_of_elementary_functions Derivatives of elementary functions] +may help. -Using the rule that the derivative of x^n for positive n (actually all nonzero n) is nx^n-1, -allows use to get the 1st differential as 3x^2. +Using the rule that the derivative of ['x[super n]] for positive n (actually all nonzero n) is ['n x[super n-1]], +allows us to get the 1st differential as ['3x[super 2]]. -To see how this extra information is used to find the root, view this demo: -[@http://en.wikipedia.org/wiki/Newton%27s_methodNewton Newton-Raphson iterations]. +To see how this extra information is used to find a root, view +[@http://en.wikipedia.org/wiki/Newton%27s_method Newton-Raphson iterations] +and the [@http://en.wikipedia.org/wiki/Newton%27s_method#mediaviewer/File:NewtonIteration_Ani.gif animation]. -We need to define a different functor that returns +We need to define a different 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 pair of floating-point values: +To \'return\' two values, we use a `std::pair` of floating-point values: */ -//[root_finding_cbrt_functor_1stderiv +//[root_finding_1_deriv_1 template -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'. +struct cbrt_functor_deriv +{ // Functor also returning 1st derviative. + 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(x) to use to get cube root of x. } - - std::pair 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. + std::pair 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 value; // to be 'cube_rooted'. -}; // cbrt_functor_1stderiv + T a; // to be 'cube_rooted'. +}; -//] [/root_finding_cbrt_functor_1stderiv] - - -/*`Our cube root function using cbrt_functor_1stderiv is now:*/ - -//[root_finding_cbrt_1deriv +/*`Our cube root function is now:*/ template -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 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::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::max)(); - // the default (std::numeric_limits::max)() = 18446744073709551615 - // which is more than we might wish to wait for!!! so we can reduce it - -/*` -Finally need to define yet another functor that returns -both the evaluation of the function to solve, -along with its first and second derivatives: - -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_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. - std::tuple operator()(T const& x) - { // Return both f(x) and f'(x) and f''(x). - 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 std:: make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x. - } -private: - T value; // to be 'cube_rooted'. -}; // struct cbrt_functor_2deriv - -//] [/root_finding_cbrt_functor_2deriv] - - -/*`Our cube function is now:*/ - -//[root_finding_cbrt_2deriv - -template -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; // halley_iterate - +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. - - int digits = std::numeric_limits::digits /2 ; // Half maximum possible binary digits accuracy for type T. - boost::uintmax_t maxit = 10; - 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; + const int digits = std::numeric_limits::digits; // Maximum possible binary digits accuracy for type T. + int get_digits = (digits * 3) /4; // Near maximum (3/4) possible accuracy. + const boost::uintmax_t maxit = 20; + boost::uintmax_t it = maxit; + T result = newton_raphson_iterate(cbrt_functor_deriv(x), guess, min, max, get_digits, it); return result; -} // cbrt_2deriv(x) +} -//] [/root_finding_cbrt_2deriv] +//] [/root_finding_1_deriv_1] + + +/* +[h3:cbrt_2_derivatives Cube root with 1st & 2nd derivative (slope & curvature)] + +Finally we define yet another functor `cbrt_functor_2deriv` that returns +both the evaluation of the function to solve, +along with its first *and second* derivatives: + +__spaces[''f](x) = 6x + +To \'return\' three values, we use a `tuple` of three floating-point values: +*/ + +//[root_finding_2deriv_1 + +template +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(x) to get cube root of x, + } + boost::math::tuple operator()(T const& x) + { // Return both f(x) and f'(x) and f''(x). + using boost::math::make_tuple; + 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 make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x. + } +private: + T a; // to be 'cube_rooted'. +}; + +/*`Our cube root function is now:*/ + +template +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::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/2; // Near maximum (1/2) possible accuracy. + boost::uintmax_t maxit = 20; + T result = halley_iterate(cbrt_functor_2deriv(x), guess, min, max, get_digits, maxit); + return result; +} + +//] [/root_finding_2deriv_1] + +/* + +[h3 Fifth-root function] +Let's now suppose we want to find the [*fifth root] of a number ['a]. + +The equation we want to solve is : + +__spaces['f](x) = x[super 5] - a + +If your differentiation is a little rusty +(or you are faced with an equation whose complexity is daunting), +then you can get help, for example from the invaluable +[@http://www.wolframalpha.com/ WolframAlpha site.] + +For example, entering the commmand: `differentiate x ^ 5` + +or the Wolfram Language command: ` D[x ^ 5, x]` + +gives the output: `d/dx(x ^ 5) = 5 x ^ 4` + +and to get the second differential, enter: `second differentiate x ^ 5` + +or the Wolfram Language command: `D[x ^ 5, { x, 2 }]` + +to get the output: `d ^ 2 / dx ^ 2(x ^ 5) = 20 x ^ 3` + +To get a reference value, we can enter: [^fifth root 3126] + +or: `N[3126 ^ (1 / 5), 50]` + +to get a result with a precision of 50 decimal digits: + +5.0003199590478625588206333405631053401128722314376 + +(We could also get a reference value using Boost.Multiprecision - see below). + +The 1st and 2nd derivatives of x[super 5] are: + +__spaces['f]\'(x) = 5x[super 4] + +__spaces['f]\'\'(x) = 20x[super 3] + +*/ + +//[root_finding_fifth_1 +//] [/root_finding_fifth_1] + + +//[root_finding_fifth_functor_2deriv + +/*`Using these expressions for the derivatives, the functor is: +*/ + +template +struct fifth_functor_2deriv +{ // Functor returning both 1st and 2nd derivatives. + fifth_functor_2deriv(T const& to_find_root_of) : a(to_find_root_of) + { // Constructor stores value a to find root of, for example: + } + + std::tuple operator()(T const& x) + { // Return both f(x) and f'(x) and f''(x). + T fx = x*x*x*x*x - a; // Difference (estimate x^3 - value). + T dx = 5 * x*x*x*x; // 1st derivative = 5x^4. + T d2x = 20 * x*x*x; // 2nd derivative = 20 x^3 + return std::make_tuple(fx, dx, d2x); // 'return' fx, dx and d2x. + } +private: + T a; // to be 'fifth_rooted'. +}; // struct fifth_functor_2deriv + +//] [/root_finding_fifth_functor_2deriv] + +//[root_finding_fifth_2deriv + +/*`Our fifth-root function is now: +*/ + +template +T fifth_2deriv(T x) +{ // return fifth 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. + + int exponent; + frexp(x, &exponent); // Get exponent of z (ignore mantissa). + T guess = ldexp(1., exponent / 5); // Rough guess is to divide the exponent by five. + T min = ldexp(0.5, exponent / 5); // Minimum possible value is half our guess. + T max = ldexp(2., exponent / 5); // Maximum possible value is twice our guess. + + const int digits = std::numeric_limits::digits / 2; // Half maximum possible binary digits accuracy for type T. + const boost::uintmax_t maxit = 50; + boost::uintmax_t it = maxit; + T result = halley_iterate(fifth_functor_2deriv(x), guess, min, max, digits, it); + return result; +} + +//] [/root_finding_fifth_2deriv] - // Default is - //boost::uintmax_t maxit = (std::numeric_limits::max)(); - // the default is (std::numeric_limits::max)() = 18446744073709551615 - // which is more than we might wish to wait for!!! so we can reduce it. int main() { + std::cout << "Root finding Examples." << std::endl; + std::cout.precision(std::numeric_limits::max_digits10); + // Show all possibly significant decimal digits for double. + //std::cout.precision(std::numeric_limits::digits10); + // Show all guaranteed significant decimal digits for double. -//[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); + +//[root_finding_main_1 try - { // Always use try'n'catch blocks with Boost.Math to get any error messages. + { + double threecubed = 27.; // value that has an exactly representable integer cube root. + double threecubedp1 = 28.; // whose cube root is *not* exactly representable. - 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. + 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(3.0365889718756625194208095785056696355814539772481111) << std::endl; - // Using bracketing: - double r = cbrt_noderiv(v27); - cout << "cbrt_noderiv(" << v27 << ") = " << r << endl; - r = cbrt_noderiv(v28); - cout << "cbrt_noderiv(" << v28 << ") = " << r << 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; +//] [/root_finding_main_1] + //[root_finding_main_2 - // Using 1st differential Newton-Raphson: - r = cbrt_1deriv(v27); - cout << "cbrt_1deriv(" << v27 << ") = " << r << endl; - r = cbrt_1deriv(v28); - cout << "cbrt_1deriv(" << v28 << ") = " << r << endl; + // Cube root using 1st differential Newton-Raphson: + r = cbrt_deriv(threecubed); + std::cout << "cbrt_deriv(" << threecubed << ") = " << r << std::endl; + r = cbrt_deriv(threecubedp1); + std::cout << "cbrt_deriv(" << threecubedp1 << ") = " << r << std::endl; + + // Cube root using Halley with 1st and 2nd differentials. + r = cbrt_2deriv(threecubed); + std::cout << "cbrt_2deriv(" << threecubed << ") = " << r << std::endl; + r = cbrt_2deriv(threecubedp1); + std::cout << "cbrt_2deriv(" << threecubedp1 << ") = " << r << std::endl; + + // Fifth root. + + double fivepowfive = 3125; // Example of a value that has an exact integer fifth root. + // Exact value of fifth root is exactly 5. + std::cout << "Fifth root of " << fivepowfive << " is " << 5 << std::endl; + + double fivepowfivep1 = fivepowfive + 1; // Example of a value whose fifth root is *not* exactly representable. + // Value of fifth root is 5.0003199590478625588206333405631053401128722314376 (50 decimal digits precision) + // and to std::numeric_limits::max_digits10 double precision (usually 17) is + + double root5v2 = static_cast(5.0003199590478625588206333405631053401128722314376); + std::cout << "Fifth root of " << fivepowfivep1 << " is " << root5v2 << std::endl; // Using Halley with 1st and 2nd differentials. - r = cbrt_2deriv(v27); - cout << "cbrt_2deriv(" << v27 << ") = " << r << endl; - r = cbrt_2deriv(v28); - cout << "cbrt_2deriv(" << v28 << ") = " << r << endl; + r = fifth_2deriv(fivepowfive); + std::cout << "fifth_2deriv(" << fivepowfive << ") = " << r << std::endl; + r = fifth_2deriv(fivepowfivep1); + std::cout << "fifth_2deriv(" << fivepowfivep1 << ") = " << r << std::endl; +//] [/root_finding_main_?] } catch(const std::exception& e) - { // Always useful to include try & catch blocks because default policies - // are to throw exceptions on arguments that cause errors like underflow, overflow. + { // Always useful to include try & catch blocks because default policies + // are to throw exceptions on arguments that cause errors like underflow, overflow. // Lacking try & catch blocks, the program will abort without a message below, // which may give some helpful clues as to the cause of the exception. std::cout << "\n""Message from thrown exception was:\n " << e.what() << std::endl; } -//] [/root_finding_example_1 return 0; } // int main() //[root_finding_example_output /*` -Normal output is: +Normal output is: [pre - Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Release\root_finding_example.exe"1> Cube Root finding (cbrt) Example. + 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. Iterations 10 - cbrt_noderiv(27) = 3 + cbrt_1(27) = 3 Iterations 10 Unable to locate solution in chosen iterations: Current best guess is between 3.0365889718756613 and 3.0365889718756627 - cbrt_noderiv(28) = 3.0365889718756618 - cbrt_1deriv(27) = 3 - cbrt_1deriv(28) = 3.0365889718756627 - cbrt_2deriv(27) = 3 - cbrt_2deriv(28) = 3.0365889718756627 -] -[/pre] + 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] to get some (much!) diagnostic output we can add @@ -374,3 +495,10 @@ to get some (much!) diagnostic output we can add ] */ //] [/root_finding_example_output] + +/* + +cbrt(28) 3.0365889718756622 +std::cbrt(28) 3.0365889718756627 + +*/ diff --git a/example/root_finding_n_example.cpp b/example/root_finding_n_example.cpp index f704a27f0..c01198536 100644 --- a/example/root_finding_n_example.cpp +++ b/example/root_finding_n_example.cpp @@ -62,7 +62,7 @@ struct nth_functor_2deriv } private: T a; // to be 'nth_rooted'. -}; +}; //] [/root_finding_nth_functor_2deriv] @@ -119,8 +119,8 @@ T show_nth_root(T value) //std::cout.precision(std::numeric_limits::max_digits10); // or use cout.precision(max_digits10 = 2 + std::numeric_limits::digits * 3010/10000); // Or guaranteed significant digits: - std::cout.precision(std::numeric_limits::digits10); - + std::cout.precision(std::numeric_limits::digits10); + T r = nth_2deriv(value); std::cout << "Type " << typeid(T).name() << " value = " << value << ", " << N << "th root = " << r << std::endl; return r; @@ -165,8 +165,8 @@ int main() // Type double value = 2, 1000000th root = 1.00000069314783 } catch (const std::exception& e) - { // Always useful to include try & catch blocks because default policies - // are to throw exceptions on arguments that cause errors like underflow, overflow. + { // Always useful to include try & catch blocks because default policies + // are to throw exceptions on arguments that cause errors like underflow, overflow. // Lacking try & catch blocks, the program will abort without a message below, // which may give some helpful clues as to the cause of the exception. std::cout << @@ -177,8 +177,15 @@ int main() /* +//[root_finding_example_output_1 +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,1> value = 2, + 5th root = 1.1486983549970350067986269467779275894438508890978 +Type class boost::multiprecision::number,0> value = 2, + 5th root = 1.1486983549970350067986269467779275894438508890978 - - +//] [/root_finding_example_output_1] */ \ No newline at end of file