1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Finding the Cubed Root With and Without Derivatives</title> 5<link rel="stylesheet" href="../../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../root_finding_examples.html" title="Examples of Root-Finding (with and without derivatives)"> 9<link rel="prev" href="../root_finding_examples.html" title="Examples of Root-Finding (with and without derivatives)"> 10<link rel="next" href="lambda.html" title="Using C++11 Lambda's"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="../root_finding_examples.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../root_finding_examples.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="lambda.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h3 class="title"> 27<a name="math_toolkit.root_finding_examples.cbrt_eg"></a><a class="link" href="cbrt_eg.html" title="Finding the Cubed Root With and Without Derivatives">Finding the 28 Cubed Root With and Without Derivatives</a> 29</h3></div></div></div> 30<p> 31 First some <code class="computeroutput"><span class="preprocessor">#includes</span></code> that 32 will be needed. 33 </p> 34<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">tools</span><span class="special">/</span><span class="identifier">roots</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 35<span class="comment">//using boost::math::policies::policy;</span> 36<span class="comment">//using boost::math::tools::newton_raphson_iterate;</span> 37<span class="comment">//using boost::math::tools::halley_iterate; //</span> 38<span class="comment">//using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits.</span> 39<span class="comment">//using boost::math::tools::bracket_and_solve_root;</span> 40<span class="comment">//using boost::math::tools::toms748_solve;</span> 41 42<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">next</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> <span class="comment">// For float_distance.</span> 43<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">tuple</span><span class="special">></span> <span class="comment">// for std::tuple and std::make_tuple.</span> 44<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">cbrt</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> <span class="comment">// For boost::math::cbrt.</span> 45</pre> 46<div class="tip"><table border="0" summary="Tip"> 47<tr> 48<td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../../doc/src/images/tip.png"></td> 49<th align="left">Tip</th> 50</tr> 51<tr><td align="left" valign="top"><p> 52 For clarity, <code class="computeroutput"><span class="keyword">using</span></code> statements 53 are provided to list what functions are being used in this example: you 54 can, of course, partly or fully qualify the names in other ways. (For your 55 application, you may wish to extract some parts into header files, but 56 you should never use <code class="computeroutput"><span class="keyword">using</span></code> 57 statements globally in header files). 58 </p></td></tr> 59</table></div> 60<p> 61 Let's suppose we want to find the root of a number <span class="emphasis"><em>a</em></span>, 62 and to start, compute the cube root. 63 </p> 64<p> 65 So the equation we want to solve is: 66 </p> 67<div class="blockquote"><blockquote class="blockquote"><p> 68 <span class="serif_italic"><span class="emphasis"><em>f(x) = x³ -a</em></span></span> 69 </p></blockquote></div> 70<p> 71 We will first solve this without using any information about the slope or 72 curvature of the cube root function. 73 </p> 74<p> 75 Fortunately, the cube-root function is 'Really Well Behaved' in that it is 76 monotonic and has only one root (we leave negative values 'as an exercise 77 for the student'). 78 </p> 79<p> 80 We then show how adding what we can know about this function, first just 81 the slope or 1st derivative <span class="emphasis"><em>f'(x)</em></span>, will speed homing 82 in on the solution. 83 </p> 84<p> 85 Lastly, we show how adding the curvature <span class="emphasis"><em>f''(x)</em></span> too 86 will speed convergence even more. 87 </p> 88<h4> 89<a name="math_toolkit.root_finding_examples.cbrt_eg.h0"></a> 90 <span class="phrase"><a name="math_toolkit.root_finding_examples.cbrt_eg.cbrt_no_derivatives"></a></span><a class="link" href="cbrt_eg.html#math_toolkit.root_finding_examples.cbrt_eg.cbrt_no_derivatives">Cube 91 root function without derivatives</a> 92 </h4> 93<p> 94 First we define a function object (functor): 95 </p> 96<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 97<span class="keyword">struct</span> <span class="identifier">cbrt_functor_noderiv</span> 98<span class="special">{</span> 99 <span class="comment">// cube root of x using only function - no derivatives.</span> 100 <span class="identifier">cbrt_functor_noderiv</span><span class="special">(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">to_find_root_of</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">a</span><span class="special">(</span><span class="identifier">to_find_root_of</span><span class="special">)</span> 101 <span class="special">{</span> <span class="comment">/* Constructor just stores value a to find root of. */</span> <span class="special">}</span> 102 <span class="identifier">T</span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">x</span><span class="special">)</span> 103 <span class="special">{</span> 104 <span class="identifier">T</span> <span class="identifier">fx</span> <span class="special">=</span> <span class="identifier">x</span><span class="special">*</span><span class="identifier">x</span><span class="special">*</span><span class="identifier">x</span> <span class="special">-</span> <span class="identifier">a</span><span class="special">;</span> <span class="comment">// Difference (estimate x^3 - a).</span> 105 <span class="keyword">return</span> <span class="identifier">fx</span><span class="special">;</span> 106 <span class="special">}</span> 107<span class="keyword">private</span><span class="special">:</span> 108 <span class="identifier">T</span> <span class="identifier">a</span><span class="special">;</span> <span class="comment">// to be 'cube_rooted'.</span> 109<span class="special">};</span> 110</pre> 111<p> 112 Implementing the cube-root function itself is fairly trivial now: the hardest 113 part is finding a good approximation to begin with. In this case we'll just 114 divide the exponent by three. (There are better but more complex guess algorithms 115 used in 'real life'.) 116 </p> 117<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 118<span class="identifier">T</span> <span class="identifier">cbrt_noderiv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">x</span><span class="special">)</span> 119<span class="special">{</span> 120 <span class="comment">// return cube root of x using bracket_and_solve (no derivatives).</span> 121 <span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">std</span><span class="special">;</span> <span class="comment">// Help ADL of std functions.</span> 122 <span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">;</span> <span class="comment">// For bracket_and_solve_root.</span> 123 124 <span class="keyword">int</span> <span class="identifier">exponent</span><span class="special">;</span> 125 <span class="identifier">frexp</span><span class="special">(</span><span class="identifier">x</span><span class="special">,</span> <span class="special">&</span><span class="identifier">exponent</span><span class="special">);</span> <span class="comment">// Get exponent of z (ignore mantissa).</span> 126 <span class="identifier">T</span> <span class="identifier">guess</span> <span class="special">=</span> <span class="identifier">ldexp</span><span class="special">(</span><span class="number">1.</span><span class="special">,</span> <span class="identifier">exponent</span><span class="special">/</span><span class="number">3</span><span class="special">);</span> <span class="comment">// Rough guess is to divide the exponent by three.</span> 127 <span class="identifier">T</span> <span class="identifier">factor</span> <span class="special">=</span> <span class="number">2</span><span class="special">;</span> <span class="comment">// How big steps to take when searching.</span> 128 129 <span class="keyword">const</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">maxit</span> <span class="special">=</span> <span class="number">20</span><span class="special">;</span> <span class="comment">// Limit to maximum iterations.</span> 130 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span> <span class="comment">// Initially our chosen max iterations, but updated with actual.</span> 131 <span class="keyword">bool</span> <span class="identifier">is_rising</span> <span class="special">=</span> <span class="keyword">true</span><span class="special">;</span> <span class="comment">// So if result if guess^3 is too low, then try increasing guess.</span> 132 <span class="keyword">int</span> <span class="identifier">digits</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">digits</span><span class="special">;</span> <span class="comment">// Maximum possible binary digits accuracy for type T.</span> 133 <span class="comment">// Some fraction of digits is used to control how accurate to try to make the result.</span> 134 <span class="keyword">int</span> <span class="identifier">get_digits</span> <span class="special">=</span> <span class="identifier">digits</span> <span class="special">-</span> <span class="number">3</span><span class="special">;</span> <span class="comment">// We have to have a non-zero interval at each step, so</span> 135 <span class="comment">// maximum accuracy is digits - 1. But we also have to</span> 136 <span class="comment">// allow for inaccuracy in f(x), otherwise the last few</span> 137 <span class="comment">// iterations just thrash around.</span> 138 <span class="identifier">eps_tolerance</span><span class="special"><</span><span class="identifier">T</span><span class="special">></span> <span class="identifier">tol</span><span class="special">(</span><span class="identifier">get_digits</span><span class="special">);</span> <span class="comment">// Set the tolerance.</span> 139 <span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="identifier">T</span><span class="special">,</span> <span class="identifier">T</span><span class="special">></span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">bracket_and_solve_root</span><span class="special">(</span><span class="identifier">cbrt_functor_noderiv</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="identifier">x</span><span class="special">),</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">factor</span><span class="special">,</span> <span class="identifier">is_rising</span><span class="special">,</span> <span class="identifier">tol</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span> 140 <span class="keyword">return</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special">+</span> <span class="special">(</span><span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span> <span class="special">-</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span><span class="special">)/</span><span class="number">2</span><span class="special">;</span> <span class="comment">// Midway between brackets is our result, if necessary we could</span> 141 <span class="comment">// return the result as an interval here.</span> 142<span class="special">}</span> 143</pre> 144<div class="note"><table border="0" summary="Note"> 145<tr> 146<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../../doc/src/images/note.png"></td> 147<th align="left">Note</th> 148</tr> 149<tr><td align="left" valign="top"> 150<p> 151 The final parameter specifying a maximum number of iterations is optional. 152 However, it defaults to <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> 153 <span class="identifier">maxit</span> <span class="special">=</span> 154 <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">>::</span><span class="identifier">max</span><span class="special">)();</span></code> which is <code class="computeroutput"><span class="number">18446744073709551615</span></code> 155 and is more than anyone would wish to wait for! 156 </p> 157<p> 158 So it may be wise to chose some reasonable estimate of how many iterations 159 may be needed, In this case the function is so well behaved that we can 160 chose a low value of 20. 161 </p> 162<p> 163 Internally when Boost.Math uses these functions, it sets the maximum iterations 164 to <code class="computeroutput"><span class="identifier">policies</span><span class="special">::</span><span class="identifier">get_max_root_iterations</span><span class="special"><</span><span class="identifier">Policy</span><span class="special">>();</span></code>. 165 </p> 166</td></tr> 167</table></div> 168<p> 169 Should we have wished we can show how many iterations were used in <code class="computeroutput"><span class="identifier">bracket_and_solve_root</span></code> (this information 170 is lost outside <code class="computeroutput"><span class="identifier">cbrt_noderiv</span></code>), 171 for example with: 172 </p> 173<pre class="programlisting"><span class="keyword">if</span> <span class="special">(</span><span class="identifier">it</span> <span class="special">>=</span> <span class="identifier">maxit</span><span class="special">)</span> 174<span class="special">{</span> 175 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Unable to locate solution in "</span> <span class="special"><<</span> <span class="identifier">maxit</span> <span class="special"><<</span> <span class="string">" iterations:"</span> 176 <span class="string">" Current best guess is between "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">" and "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 177<span class="special">}</span> 178<span class="keyword">else</span> 179<span class="special">{</span> 180 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Converged after "</span> <span class="special"><<</span> <span class="identifier">it</span> <span class="special"><<</span> <span class="string">" (from maximum of "</span> <span class="special"><<</span> <span class="identifier">maxit</span> <span class="special"><<</span> <span class="string">" iterations)."</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 181<span class="special">}</span> 182</pre> 183<p> 184 for output like 185 </p> 186<pre class="programlisting"><span class="identifier">Converged</span> <span class="identifier">after</span> <span class="number">11</span> <span class="special">(</span><span class="identifier">from</span> <span class="identifier">maximum</span> <span class="identifier">of</span> <span class="number">20</span> <span class="identifier">iterations</span><span class="special">).</span> 187</pre> 188<p> 189 This snippet from <code class="computeroutput"><span class="identifier">main</span><span class="special">()</span></code> 190 in <a href="../../../../example/root_finding_example.cpp" target="_top">root_finding_example.cpp</a> 191 shows how it can be used. 192 </p> 193<pre class="programlisting"><span class="keyword">try</span> 194<span class="special">{</span> 195 <span class="keyword">double</span> <span class="identifier">threecubed</span> <span class="special">=</span> <span class="number">27.</span><span class="special">;</span> <span class="comment">// Value that has an *exactly representable* integer cube root.</span> 196 <span class="keyword">double</span> <span class="identifier">threecubedp1</span> <span class="special">=</span> <span class="number">28.</span><span class="special">;</span> <span class="comment">// Value whose cube root is *not* exactly representable.</span> 197 198 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"cbrt(28) "</span> <span class="special"><<</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cbrt</span><span class="special">(</span><span class="number">28.</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// boost::math:: version of cbrt.</span> 199 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"std::cbrt(28) "</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cbrt</span><span class="special">(</span><span class="number">28.</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// std:: version of cbrt.</span> 200 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span><span class="string">" cast double "</span> <span class="special"><<</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">double</span><span class="special">>(</span><span class="number">3.0365889718756625194208095785056696355814539772481111</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 201 202 <span class="comment">// Cube root using bracketing:</span> 203 <span class="keyword">double</span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">cbrt_noderiv</span><span class="special">(</span><span class="identifier">threecubed</span><span class="special">);</span> 204 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"cbrt_noderiv("</span> <span class="special"><<</span> <span class="identifier">threecubed</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">r</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 205 <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">cbrt_noderiv</span><span class="special">(</span><span class="identifier">threecubedp1</span><span class="special">);</span> 206 <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"cbrt_noderiv("</span> <span class="special"><<</span> <span class="identifier">threecubedp1</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">r</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> 207</pre> 208<pre class="programlisting"> cbrt_noderiv(27) = 3 209 cbrt_noderiv(28) = 3.0365889718756618 210</pre> 211<p> 212 The result of <code class="computeroutput"><span class="identifier">bracket_and_solve_root</span></code> 213 is a <a href="http://www.cplusplus.com/reference/utility/pair/" target="_top">pair</a> 214 of values that could be displayed. 215 </p> 216<p> 217 The number of bits separating them can be found using <code class="computeroutput"><span class="identifier">float_distance</span><span class="special">(</span><span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span><span class="special">,</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span><span class="special">)</span></code>. The distance is zero (closest representable) 218 for 3<sup>3</sup> = 27 but <code class="computeroutput"><span class="identifier">float_distance</span><span class="special">(</span><span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span><span class="special">,</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span><span class="special">)</span> <span class="special">=</span> <span class="number">3</span></code> 219 for cube root of 28 with this function. The result (avoiding overflow) is 220 midway between these two values. 221 </p> 222<h4> 223<a name="math_toolkit.root_finding_examples.cbrt_eg.h1"></a> 224 <span class="phrase"><a name="math_toolkit.root_finding_examples.cbrt_eg.cbrt_1st_derivative"></a></span><a class="link" href="cbrt_eg.html#math_toolkit.root_finding_examples.cbrt_eg.cbrt_1st_derivative">Cube 225 root function with 1st derivative (slope)</a> 226 </h4> 227<p> 228 We now solve the same problem, but using more information about the function, 229 to show how this can speed up finding the best estimate of the root. 230 </p> 231<p> 232 For the root function, the 1st differential (the slope of the tangent to 233 a curve at any point) is known. 234 </p> 235<p> 236 This algorithm is similar to this <a href="http://en.wikipedia.org/wiki/Nth_root_algorithm" target="_top">nth 237 root algorithm</a>. 238 </p> 239<p> 240 If you need some reminders, then <a href="http://en.wikipedia.org/wiki/Derivative#Derivatives_of_elementary_functions" target="_top">derivatives 241 of elementary functions</a> may help. 242 </p> 243<p> 244 Using the rule that the derivative of <span class="emphasis"><em>x<sup>n</sup></em></span> for positive 245 n (actually all nonzero n) is <span class="emphasis"><em>n x<sup>n-1</sup></em></span>, allows us to get 246 the 1st differential as <span class="emphasis"><em>3x<sup>2</sup></em></span>. 247 </p> 248<p> 249 To see how this extra information is used to find a root, view <a href="http://en.wikipedia.org/wiki/Newton%27s_method" target="_top">Newton-Raphson 250 iterations</a> and the <a href="http://en.wikipedia.org/wiki/Newton%27s_method#mediaviewer/File:NewtonIteration_Ani.gif" target="_top">animation</a>. 251 </p> 252<p> 253 We define a better functor <code class="computeroutput"><span class="identifier">cbrt_functor_deriv</span></code> 254 that returns both the evaluation of the function to solve, along with its 255 first derivative: 256 </p> 257<p> 258 To '<span class="emphasis"><em>return</em></span>' two values, we use a <a href="http://en.cppreference.com/w/cpp/utility/pair" target="_top">std::pair</a> 259 of floating-point values. 260 </p> 261<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 262<span class="keyword">struct</span> <span class="identifier">cbrt_functor_deriv</span> 263<span class="special">{</span> <span class="comment">// Functor also returning 1st derivative.</span> 264 <span class="identifier">cbrt_functor_deriv</span><span class="special">(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">to_find_root_of</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">a</span><span class="special">(</span><span class="identifier">to_find_root_of</span><span class="special">)</span> 265 <span class="special">{</span> <span class="comment">// Constructor stores value a to find root of,</span> 266 <span class="comment">// for example: calling cbrt_functor_deriv<T>(a) to use to get cube root of a.</span> 267 <span class="special">}</span> 268 <span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="identifier">T</span><span class="special">,</span> <span class="identifier">T</span><span class="special">></span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">x</span><span class="special">)</span> 269 <span class="special">{</span> 270 <span class="comment">// Return both f(x) and f'(x).</span> 271 <span class="identifier">T</span> <span class="identifier">fx</span> <span class="special">=</span> <span class="identifier">x</span><span class="special">*</span><span class="identifier">x</span><span class="special">*</span><span class="identifier">x</span> <span class="special">-</span> <span class="identifier">a</span><span class="special">;</span> <span class="comment">// Difference (estimate x^3 - value).</span> 272 <span class="identifier">T</span> <span class="identifier">dx</span> <span class="special">=</span> <span class="number">3</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">*</span><span class="identifier">x</span><span class="special">;</span> <span class="comment">// 1st derivative = 3x^2.</span> 273 <span class="keyword">return</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">make_pair</span><span class="special">(</span><span class="identifier">fx</span><span class="special">,</span> <span class="identifier">dx</span><span class="special">);</span> <span class="comment">// 'return' both fx and dx.</span> 274 <span class="special">}</span> 275<span class="keyword">private</span><span class="special">:</span> 276 <span class="identifier">T</span> <span class="identifier">a</span><span class="special">;</span> <span class="comment">// Store value to be 'cube_rooted'.</span> 277<span class="special">};</span> 278</pre> 279<p> 280 Our cube root function is now: 281 </p> 282<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 283<span class="identifier">T</span> <span class="identifier">cbrt_deriv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">x</span><span class="special">)</span> 284<span class="special">{</span> 285 <span class="comment">// return cube root of x using 1st derivative and Newton_Raphson.</span> 286 <span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">;</span> 287 <span class="keyword">int</span> <span class="identifier">exponent</span><span class="special">;</span> 288 <span class="identifier">frexp</span><span class="special">(</span><span class="identifier">x</span><span class="special">,</span> <span class="special">&</span><span class="identifier">exponent</span><span class="special">);</span> <span class="comment">// Get exponent of z (ignore mantissa).</span> 289 <span class="identifier">T</span> <span class="identifier">guess</span> <span class="special">=</span> <span class="identifier">ldexp</span><span class="special">(</span><span class="number">1.</span><span class="special">,</span> <span class="identifier">exponent</span><span class="special">/</span><span class="number">3</span><span class="special">);</span> <span class="comment">// Rough guess is to divide the exponent by three.</span> 290 <span class="identifier">T</span> <span class="identifier">min</span> <span class="special">=</span> <span class="identifier">ldexp</span><span class="special">(</span><span class="number">0.5</span><span class="special">,</span> <span class="identifier">exponent</span><span class="special">/</span><span class="number">3</span><span class="special">);</span> <span class="comment">// Minimum possible value is half our guess.</span> 291 <span class="identifier">T</span> <span class="identifier">max</span> <span class="special">=</span> <span class="identifier">ldexp</span><span class="special">(</span><span class="number">2.</span><span class="special">,</span> <span class="identifier">exponent</span><span class="special">/</span><span class="number">3</span><span class="special">);</span> <span class="comment">// Maximum possible value is twice our guess.</span> 292 <span class="keyword">const</span> <span class="keyword">int</span> <span class="identifier">digits</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">digits</span><span class="special">;</span> <span class="comment">// Maximum possible binary digits accuracy for type T.</span> 293 <span class="keyword">int</span> <span class="identifier">get_digits</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">int</span><span class="special">>(</span><span class="identifier">digits</span> <span class="special">*</span> <span class="number">0.6</span><span class="special">);</span> <span class="comment">// Accuracy doubles with each step, so stop when we have</span> 294 <span class="comment">// just over half the digits correct.</span> 295 <span class="keyword">const</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">maxit</span> <span class="special">=</span> <span class="number">20</span><span class="special">;</span> 296 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span> 297 <span class="identifier">T</span> <span class="identifier">result</span> <span class="special">=</span> <span class="identifier">newton_raphson_iterate</span><span class="special">(</span><span class="identifier">cbrt_functor_deriv</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="identifier">x</span><span class="special">),</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">max</span><span class="special">,</span> <span class="identifier">get_digits</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span> 298 <span class="keyword">return</span> <span class="identifier">result</span><span class="special">;</span> 299<span class="special">}</span> 300</pre> 301<p> 302 The result of <a href="../../../../../../libs/math/include/boost/math/tools/roots.hpp" target="_top"><code class="computeroutput"><span class="identifier">newton_raphson_iterate</span></code></a> function 303 is a single value. 304 </p> 305<div class="tip"><table border="0" summary="Tip"> 306<tr> 307<td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../../doc/src/images/tip.png"></td> 308<th align="left">Tip</th> 309</tr> 310<tr><td align="left" valign="top"><p> 311 There is a compromise between accuracy and speed when choosing the value 312 of <code class="computeroutput"><span class="identifier">digits</span></code>. It is tempting 313 to simply chose <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">digits</span></code>, but this may mean some inefficient 314 and unnecessary iterations as the function thrashes around trying to locate 315 the last bit. In theory, since the precision doubles with each step it 316 is sufficient to stop when half the bits are correct: as the last step 317 will have doubled that to full precision. Of course the function has no 318 way to tell if that is actually the case unless it does one more step to 319 be sure. In practice setting the precision to slightly more than <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">digits</span> <span class="special">/</span> <span class="number">2</span></code> is a good choice. 320 </p></td></tr> 321</table></div> 322<p> 323 Note that it is up to the caller of the function to check the iteration count 324 after the call to see if iteration stoped as a result of running out of iterations 325 rather than meeting the required precision. 326 </p> 327<p> 328 Using the test data in <a href="../../../../test/test_cbrt.cpp" target="_top">/test/test_cbrt.cpp</a> 329 this found the cube root exact to the last digit in every case, and in no 330 more than 6 iterations at double precision. However, you will note that a 331 high precision was used in this example, exactly what was warned against 332 earlier on in these docs! In this particular case it is possible to compute 333 <span class="emphasis"><em>f(x)</em></span> exactly and without undue cancellation error, so 334 a high limit is not too much of an issue. 335 </p> 336<p> 337 However, reducing the limit to <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">digits</span> <span class="special">*</span> <span class="number">2</span> <span class="special">/</span> <span class="number">3</span></code> 338 gave full precision in all but one of the test cases (and that one was out 339 by just one bit). The maximum number of iterations remained 6, but in most 340 cases was reduced by one. 341 </p> 342<p> 343 Note also that the above code omits a probable optimization by computing 344 z² 345and reusing it, omits error handling, and does not handle negative values 346 of z correctly. (These are left as the customary exercise for the reader!) 347 </p> 348<p> 349 The <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cbrt</span></code> function also includes these and other 350 improvements: most importantly it uses a much better initial guess which 351 reduces the iteration count to just 1 in almost all cases. 352 </p> 353<h4> 354<a name="math_toolkit.root_finding_examples.cbrt_eg.h2"></a> 355 <span class="phrase"><a name="math_toolkit.root_finding_examples.cbrt_eg.cbrt_2_derivatives"></a></span><a class="link" href="cbrt_eg.html#math_toolkit.root_finding_examples.cbrt_eg.cbrt_2_derivatives">Cube 356 root with 1st & 2nd derivative (slope & curvature)</a> 357 </h4> 358<p> 359 Next we define yet another even better functor <code class="computeroutput"><span class="identifier">cbrt_functor_2deriv</span></code> 360 that returns both the evaluation of the function to solve, along with its 361 first <span class="bold"><strong>and second</strong></span> derivative: 362 </p> 363<div class="blockquote"><blockquote class="blockquote"><p> 364 <span class="serif_italic"><span class="emphasis"><em>f''(x) = 6x</em></span></span> 365 </p></blockquote></div> 366<p> 367 using information about both slope and curvature to speed convergence. 368 </p> 369<p> 370 To <span class="emphasis"><em>'return'</em></span> three values, we use a <code class="computeroutput"><span class="identifier">tuple</span></code> 371 of three floating-point values: 372 </p> 373<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 374<span class="keyword">struct</span> <span class="identifier">cbrt_functor_2deriv</span> 375<span class="special">{</span> 376 <span class="comment">// Functor returning both 1st and 2nd derivatives.</span> 377 <span class="identifier">cbrt_functor_2deriv</span><span class="special">(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">to_find_root_of</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">a</span><span class="special">(</span><span class="identifier">to_find_root_of</span><span class="special">)</span> 378 <span class="special">{</span> <span class="comment">// Constructor stores value a to find root of, for example:</span> 379 <span class="comment">// calling cbrt_functor_2deriv<T>(x) to get cube root of x,</span> 380 <span class="special">}</span> 381 <span class="identifier">std</span><span class="special">::</span><span class="identifier">tuple</span><span class="special"><</span><span class="identifier">T</span><span class="special">,</span> <span class="identifier">T</span><span class="special">,</span> <span class="identifier">T</span><span class="special">></span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">x</span><span class="special">)</span> 382 <span class="special">{</span> 383 <span class="comment">// Return both f(x) and f'(x) and f''(x).</span> 384 <span class="identifier">T</span> <span class="identifier">fx</span> <span class="special">=</span> <span class="identifier">x</span><span class="special">*</span><span class="identifier">x</span><span class="special">*</span><span class="identifier">x</span> <span class="special">-</span> <span class="identifier">a</span><span class="special">;</span> <span class="comment">// Difference (estimate x^3 - value).</span> 385 <span class="identifier">T</span> <span class="identifier">dx</span> <span class="special">=</span> <span class="number">3</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">*</span><span class="identifier">x</span><span class="special">;</span> <span class="comment">// 1st derivative = 3x^2.</span> 386 <span class="identifier">T</span> <span class="identifier">d2x</span> <span class="special">=</span> <span class="number">6</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">;</span> <span class="comment">// 2nd derivative = 6x.</span> 387 <span class="keyword">return</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">make_tuple</span><span class="special">(</span><span class="identifier">fx</span><span class="special">,</span> <span class="identifier">dx</span><span class="special">,</span> <span class="identifier">d2x</span><span class="special">);</span> <span class="comment">// 'return' fx, dx and d2x.</span> 388 <span class="special">}</span> 389<span class="keyword">private</span><span class="special">:</span> 390 <span class="identifier">T</span> <span class="identifier">a</span><span class="special">;</span> <span class="comment">// to be 'cube_rooted'.</span> 391<span class="special">};</span> 392</pre> 393<p> 394 Our cube root function is now: 395 </p> 396<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 397<span class="identifier">T</span> <span class="identifier">cbrt_2deriv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">x</span><span class="special">)</span> 398<span class="special">{</span> 399 <span class="comment">// return cube root of x using 1st and 2nd derivatives and Halley.</span> 400 <span class="comment">//using namespace std; // Help ADL of std functions.</span> 401 <span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">;</span> 402 <span class="keyword">int</span> <span class="identifier">exponent</span><span class="special">;</span> 403 <span class="identifier">frexp</span><span class="special">(</span><span class="identifier">x</span><span class="special">,</span> <span class="special">&</span><span class="identifier">exponent</span><span class="special">);</span> <span class="comment">// Get exponent of z (ignore mantissa).</span> 404 <span class="identifier">T</span> <span class="identifier">guess</span> <span class="special">=</span> <span class="identifier">ldexp</span><span class="special">(</span><span class="number">1.</span><span class="special">,</span> <span class="identifier">exponent</span><span class="special">/</span><span class="number">3</span><span class="special">);</span> <span class="comment">// Rough guess is to divide the exponent by three.</span> 405 <span class="identifier">T</span> <span class="identifier">min</span> <span class="special">=</span> <span class="identifier">ldexp</span><span class="special">(</span><span class="number">0.5</span><span class="special">,</span> <span class="identifier">exponent</span><span class="special">/</span><span class="number">3</span><span class="special">);</span> <span class="comment">// Minimum possible value is half our guess.</span> 406 <span class="identifier">T</span> <span class="identifier">max</span> <span class="special">=</span> <span class="identifier">ldexp</span><span class="special">(</span><span class="number">2.</span><span class="special">,</span> <span class="identifier">exponent</span><span class="special">/</span><span class="number">3</span><span class="special">);</span> <span class="comment">// Maximum possible value is twice our guess.</span> 407 <span class="keyword">const</span> <span class="keyword">int</span> <span class="identifier">digits</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">digits</span><span class="special">;</span> <span class="comment">// Maximum possible binary digits accuracy for type T.</span> 408 <span class="comment">// digits used to control how accurate to try to make the result.</span> 409 <span class="keyword">int</span> <span class="identifier">get_digits</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">int</span><span class="special">>(</span><span class="identifier">digits</span> <span class="special">*</span> <span class="number">0.4</span><span class="special">);</span> <span class="comment">// Accuracy triples with each step, so stop when just</span> 410 <span class="comment">// over one third of the digits are correct.</span> 411 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">maxit</span> <span class="special">=</span> <span class="number">20</span><span class="special">;</span> 412 <span class="identifier">T</span> <span class="identifier">result</span> <span class="special">=</span> <span class="identifier">halley_iterate</span><span class="special">(</span><span class="identifier">cbrt_functor_2deriv</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="identifier">x</span><span class="special">),</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">max</span><span class="special">,</span> <span class="identifier">get_digits</span><span class="special">,</span> <span class="identifier">maxit</span><span class="special">);</span> 413 <span class="keyword">return</span> <span class="identifier">result</span><span class="special">;</span> 414<span class="special">}</span> 415</pre> 416<p> 417 The function <code class="computeroutput"><span class="identifier">halley_iterate</span></code> 418 also returns a single value, and the number of iterations will reveal if 419 it met the convergence criterion set by <code class="computeroutput"><span class="identifier">get_digits</span></code>. 420 </p> 421<p> 422 The no-derivative method gives a result of 423 </p> 424<pre class="programlisting"><span class="identifier">cbrt_noderiv</span><span class="special">(</span><span class="number">28</span><span class="special">)</span> <span class="special">=</span> <span class="number">3.0365889718756618</span> 425</pre> 426<p> 427 with a 3 bits distance between the bracketed values, whereas the derivative 428 methods both converge to a single value 429 </p> 430<pre class="programlisting"><span class="identifier">cbrt_2deriv</span><span class="special">(</span><span class="number">28</span><span class="special">)</span> <span class="special">=</span> <span class="number">3.0365889718756627</span> 431</pre> 432<p> 433 which we can compare with the <a href="../../../../../../libs/math/doc/html/math_toolkit/powers/cbrt.html" target="_top">boost::math::cbrt</a> 434 </p> 435<pre class="programlisting"><span class="identifier">cbrt</span><span class="special">(</span><span class="number">28</span><span class="special">)</span> <span class="special">=</span> <span class="number">3.0365889718756627</span> 436</pre> 437<p> 438 Note that the iterations are set to stop at just one-half of full precision, 439 and yet, even so, not one of the test cases had a single bit wrong. What's 440 more, the maximum number of iterations was now just 4. 441 </p> 442<p> 443 Just to complete the picture, we could have called <a class="link" href="../roots_deriv.html#math_toolkit.roots_deriv.schroder"><code class="computeroutput"><span class="identifier">schroder_iterate</span></code></a> in the last example: 444 and in fact it makes no difference to the accuracy or number of iterations 445 in this particular case. However, the relative performance of these two methods 446 may vary depending upon the nature of <span class="emphasis"><em>f(x)</em></span>, and the 447 accuracy to which the initial guess can be computed. There appear to be no 448 generalisations that can be made except "try them and see". 449 </p> 450<p> 451 Finally, had we called <code class="computeroutput"><span class="identifier">cbrt</span></code> 452 with <a href="http://shoup.net/ntl/doc/RR.txt" target="_top">NTL::RR</a> set to 453 1000 bit precision (about 300 decimal digits), then full precision can be 454 obtained with just 7 iterations. To put that in perspective, an increase 455 in precision by a factor of 20, has less than doubled the number of iterations. 456 That just goes to emphasise that most of the iterations are used up getting 457 the first few digits correct: after that these methods can churn out further 458 digits with remarkable efficiency. 459 </p> 460<p> 461 Or to put it another way: <span class="emphasis"><em>nothing beats a really good initial guess!</em></span> 462 </p> 463<p> 464 Full code of this example is at <a href="../../../../example/root_finding_example.cpp" target="_top">root_finding_example.cpp</a>, 465 </p> 466</div> 467<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 468<td align="left"></td> 469<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 470 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 471 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 472 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 473 Daryle Walker and Xiaogang Zhang<p> 474 Distributed under the Boost Software License, Version 1.0. (See accompanying 475 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) 476 </p> 477</div></td> 478</tr></table> 479<hr> 480<div class="spirit-nav"> 481<a accesskey="p" href="../root_finding_examples.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../root_finding_examples.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="lambda.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a> 482</div> 483</body> 484</html> 485
