1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Generalizing to Compute the nth root</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="multiprecision_root.html" title="Root-finding using Boost.Multiprecision"> 10<link rel="next" href="elliptic_eg.html" title="A More complex example - Inverting the Elliptic Integrals"> 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="multiprecision_root.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="elliptic_eg.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.nth_root"></a><a class="link" href="nth_root.html" title="Generalizing to Compute the nth root">Generalizing 28 to Compute the nth root</a> 29</h3></div></div></div> 30<p> 31 If desired, we can now further generalize to compute the <span class="emphasis"><em>n</em></span>th 32 root by computing the derivatives <span class="bold"><strong>at compile-time</strong></span> 33 using the rules for differentiation and <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">pow</span><span class="special"><</span><span class="identifier">N</span><span class="special">></span></code> where 34 template parameter <code class="computeroutput"><span class="identifier">N</span></code> is an 35 integer and a compile time constant. Our functor and function now have an 36 additional template parameter <code class="computeroutput"><span class="identifier">N</span></code>, 37 for the root required. 38 </p> 39<div class="note"><table border="0" summary="Note"> 40<tr> 41<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../../doc/src/images/note.png"></td> 42<th align="left">Note</th> 43</tr> 44<tr><td align="left" valign="top"><p> 45 Since the powers and derivatives are fixed at compile time, the resulting 46 code is as efficient as as if hand-coded as the cube and fifth-root examples 47 above. A good compiler should also optimise any repeated multiplications. 48 </p></td></tr> 49</table></div> 50<p> 51 Our <span class="emphasis"><em>n</em></span>th root functor is 52 </p> 53<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">int</span> <span class="identifier">N</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span> <span class="special">=</span> <span class="keyword">double</span><span class="special">></span> 54<span class="keyword">struct</span> <span class="identifier">nth_functor_2deriv</span> 55<span class="special">{</span> <span class="comment">// Functor returning both 1st and 2nd derivatives.</span> 56 <span class="identifier">BOOST_STATIC_ASSERT_MSG</span><span class="special">(</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">is_integral</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">value</span> <span class="special">==</span> <span class="keyword">false</span><span class="special">,</span> <span class="string">"Only floating-point type types can be used!"</span><span class="special">);</span> 57 <span class="identifier">BOOST_STATIC_ASSERT_MSG</span><span class="special">((</span><span class="identifier">N</span> <span class="special">></span> <span class="number">0</span><span class="special">)</span> <span class="special">==</span> <span class="keyword">true</span><span class="special">,</span> <span class="string">"root N must be > 0!"</span><span class="special">);</span> 58 59 <span class="identifier">nth_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> 60 <span class="special">{</span> <span class="comment">/* Constructor stores value a to find root of, for example: */</span> <span class="special">}</span> 61 62 <span class="comment">// using boost::math::tuple; // to return three values.</span> 63 <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> 64 <span class="special">{</span> 65 <span class="comment">// Return f(x), f'(x) and f''(x).</span> 66 <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">pow</span><span class="special">;</span> 67 <span class="identifier">T</span> <span class="identifier">fx</span> <span class="special">=</span> <span class="identifier">pow</span><span class="special"><</span><span class="identifier">N</span><span class="special">>(</span><span class="identifier">x</span><span class="special">)</span> <span class="special">-</span> <span class="identifier">a</span><span class="special">;</span> <span class="comment">// Difference (estimate x^n - a).</span> 68 <span class="identifier">T</span> <span class="identifier">dx</span> <span class="special">=</span> <span class="identifier">N</span> <span class="special">*</span> <span class="identifier">pow</span><span class="special"><</span><span class="identifier">N</span> <span class="special">-</span> <span class="number">1</span><span class="special">>(</span><span class="identifier">x</span><span class="special">);</span> <span class="comment">// 1st derivative f'(x).</span> 69 <span class="identifier">T</span> <span class="identifier">d2x</span> <span class="special">=</span> <span class="identifier">N</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">N</span> <span class="special">-</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="identifier">pow</span><span class="special"><</span><span class="identifier">N</span> <span class="special">-</span> <span class="number">2</span> <span class="special">>(</span><span class="identifier">x</span><span class="special">);</span> <span class="comment">// 2nd derivative f''(x).</span> 70 71 <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> 72 <span class="special">}</span> 73<span class="keyword">private</span><span class="special">:</span> 74 <span class="identifier">T</span> <span class="identifier">a</span><span class="special">;</span> <span class="comment">// to be 'nth_rooted'.</span> 75<span class="special">};</span> 76</pre> 77<p> 78 and our <span class="emphasis"><em>n</em></span>th root function is 79 </p> 80<pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">int</span> <span class="identifier">N</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span> <span class="special">=</span> <span class="keyword">double</span><span class="special">></span> 81<span class="identifier">T</span> <span class="identifier">nth_2deriv</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">x</span><span class="special">)</span> 82<span class="special">{</span> <span class="comment">// return nth root of x using 1st and 2nd derivatives and Halley.</span> 83 84 <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> 85 <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 halley_iterate.</span> 86 87 <span class="identifier">BOOST_STATIC_ASSERT_MSG</span><span class="special">(</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">is_integral</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">value</span> <span class="special">==</span> <span class="keyword">false</span><span class="special">,</span> <span class="string">"Only floating-point type types can be used!"</span><span class="special">);</span> 88 <span class="identifier">BOOST_STATIC_ASSERT_MSG</span><span class="special">((</span><span class="identifier">N</span> <span class="special">></span> <span class="number">0</span><span class="special">)</span> <span class="special">==</span> <span class="keyword">true</span><span class="special">,</span> <span class="string">"root N must be > 0!"</span><span class="special">);</span> 89 <span class="identifier">BOOST_STATIC_ASSERT_MSG</span><span class="special">((</span><span class="identifier">N</span> <span class="special">></span> <span class="number">1000</span><span class="special">)</span> <span class="special">==</span> <span class="keyword">false</span><span class="special">,</span> <span class="string">"root N is too big!"</span><span class="special">);</span> 90 91 <span class="keyword">typedef</span> <span class="keyword">double</span> <span class="identifier">guess_type</span><span class="special">;</span> <span class="comment">// double may restrict (exponent) range for a multiprecision T?</span> 92 93 <span class="keyword">int</span> <span class="identifier">exponent</span><span class="special">;</span> 94 <span class="identifier">frexp</span><span class="special">(</span><span class="keyword">static_cast</span><span class="special"><</span><span class="identifier">guess_type</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> 95 <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="keyword">static_cast</span><span class="special"><</span><span class="identifier">guess_type</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="identifier">N</span><span class="special">);</span> <span class="comment">// Rough guess is to divide the exponent by n.</span> 96 <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="keyword">static_cast</span><span class="special"><</span><span class="identifier">guess_type</span><span class="special">>(</span><span class="number">1.</span><span class="special">)</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="identifier">N</span><span class="special">);</span> <span class="comment">// Minimum possible value is half our guess.</span> 97 <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="keyword">static_cast</span><span class="special"><</span><span class="identifier">guess_type</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="identifier">N</span><span class="special">);</span> <span class="comment">// Maximum possible value is twice our guess.</span> 98 99 <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="number">0.4</span><span class="special">;</span> <span class="comment">// Accuracy triples with each step, so stop when</span> 100 <span class="comment">// slightly more than one third of the digits are correct.</span> 101 <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> 102 <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> 103 <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">nth_functor_2deriv</span><span class="special"><</span><span class="identifier">N</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">digits</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span> 104 <span class="keyword">return</span> <span class="identifier">result</span><span class="special">;</span> 105<span class="special">}</span> 106</pre> 107<pre class="programlisting"> <span class="identifier">show_nth_root</span><span class="special"><</span><span class="number">5</span><span class="special">,</span> <span class="keyword">double</span><span class="special">>(</span><span class="number">2.</span><span class="special">);</span> 108 <span class="identifier">show_nth_root</span><span class="special"><</span><span class="number">5</span><span class="special">,</span> <span class="keyword">long</span> <span class="keyword">double</span><span class="special">>(</span><span class="number">2.</span><span class="special">);</span> 109<span class="preprocessor">#ifndef</span> <span class="identifier">_MSC_VER</span> <span class="comment">// float128 is not supported by Microsoft compiler 2013.</span> 110 <span class="identifier">show_nth_root</span><span class="special"><</span><span class="number">5</span><span class="special">,</span> <span class="identifier">float128</span><span class="special">>(</span><span class="number">2</span><span class="special">);</span> 111<span class="preprocessor">#endif</span> 112 <span class="identifier">show_nth_root</span><span class="special"><</span><span class="number">5</span><span class="special">,</span> <span class="identifier">cpp_dec_float_50</span><span class="special">>(</span><span class="number">2</span><span class="special">);</span> <span class="comment">// dec</span> 113 <span class="identifier">show_nth_root</span><span class="special"><</span><span class="number">5</span><span class="special">,</span> <span class="identifier">cpp_bin_float_50</span><span class="special">>(</span><span class="number">2</span><span class="special">);</span> <span class="comment">// bin</span> 114</pre> 115<p> 116 produces an output similar to this 117 </p> 118<pre class="programlisting"> <span class="identifier">Using</span> <span class="identifier">MSVC</span> <span class="number">2013</span> 119 120<span class="identifier">nth</span> <span class="identifier">Root</span> <span class="identifier">finding</span> <span class="identifier">Example</span><span class="special">.</span> 121<span class="identifier">Type</span> <span class="keyword">double</span> <span class="identifier">value</span> <span class="special">=</span> <span class="number">2</span><span class="special">,</span> <span class="number">5</span><span class="identifier">th</span> <span class="identifier">root</span> <span class="special">=</span> <span class="number">1.14869835499704</span> 122<span class="identifier">Type</span> <span class="keyword">long</span> <span class="keyword">double</span> <span class="identifier">value</span> <span class="special">=</span> <span class="number">2</span><span class="special">,</span> <span class="number">5</span><span class="identifier">th</span> <span class="identifier">root</span> <span class="special">=</span> <span class="number">1.14869835499704</span> 123<span class="identifier">Type</span> <span class="keyword">class</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">number</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">backends</span><span class="special">::</span><span class="identifier">cpp_dec_float</span><span class="special"><</span><span class="number">50</span><span class="special">,</span><span class="keyword">int</span><span class="special">,</span><span class="keyword">void</span><span class="special">>,</span><span class="number">1</span><span class="special">></span> <span class="identifier">value</span> <span class="special">=</span> <span class="number">2</span><span class="special">,</span> 124 <span class="number">5</span><span class="identifier">th</span> <span class="identifier">root</span> <span class="special">=</span> <span class="number">1.1486983549970350067986269467779275894438508890978</span> 125<span class="identifier">Type</span> <span class="keyword">class</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">number</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">backends</span><span class="special">::</span><span class="identifier">cpp_bin_float</span><span class="special"><</span><span class="number">50</span><span class="special">,</span><span class="number">10</span><span class="special">,</span><span class="keyword">void</span><span class="special">,</span><span class="keyword">int</span><span class="special">,</span><span class="number">0</span><span class="special">,</span><span class="number">0</span><span class="special">>,</span><span class="number">0</span><span class="special">></span> <span class="identifier">value</span> <span class="special">=</span> <span class="number">2</span><span class="special">,</span> 126 <span class="number">5</span><span class="identifier">th</span> <span class="identifier">root</span> <span class="special">=</span> <span class="number">1.1486983549970350067986269467779275894438508890978</span> 127</pre> 128<div class="tip"><table border="0" summary="Tip"> 129<tr> 130<td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../../doc/src/images/tip.png"></td> 131<th align="left">Tip</th> 132</tr> 133<tr><td align="left" valign="top"> 134<p> 135 Take care with the type passed to the function. It is best to pass a <code class="computeroutput"><span class="keyword">double</span></code> or greater-precision floating-point 136 type. 137 </p> 138<p> 139 Passing an integer value, for example, <code class="computeroutput"><span class="identifier">nth_2deriv</span><span class="special"><</span><span class="number">5</span><span class="special">>(</span><span class="number">2</span><span class="special">)</span></code> will be 140 rejected, while <code class="computeroutput"><span class="identifier">nth_2deriv</span><span class="special"><</span><span class="number">5</span><span class="special">,</span> 141 <span class="keyword">double</span><span class="special">>(</span><span class="number">2</span><span class="special">)</span></code> converts 142 the integer to <code class="computeroutput"><span class="keyword">double</span></code>. 143 </p> 144<p> 145 Avoid passing a <code class="computeroutput"><span class="keyword">float</span></code> value 146 that will provoke warnings (actually spurious) from the compiler about 147 potential loss of data, as noted above. 148 </p> 149</td></tr> 150</table></div> 151<div class="warning"><table border="0" summary="Warning"> 152<tr> 153<td rowspan="2" align="center" valign="top" width="25"><img alt="[Warning]" src="../../../../../../doc/src/images/warning.png"></td> 154<th align="left">Warning</th> 155</tr> 156<tr><td align="left" valign="top"><p> 157 Asking for unreasonable roots, for example, <code class="computeroutput"><span class="identifier">show_nth_root</span><span class="special"><</span><span class="number">1000000</span><span class="special">>(</span><span class="number">2.</span><span class="special">);</span></code> 158 may lead to <a href="http://en.wikipedia.org/wiki/Loss_of_significance" target="_top">Loss 159 of significance</a> like <code class="computeroutput"><span class="identifier">Type</span> 160 <span class="keyword">double</span> <span class="identifier">value</span> 161 <span class="special">=</span> <span class="number">2</span><span class="special">,</span> <span class="number">1000000</span><span class="identifier">th</span> <span class="identifier">root</span> 162 <span class="special">=</span> <span class="number">1.00000069314783</span></code>. 163 Use of the the <code class="computeroutput"><span class="identifier">pow</span></code> function 164 is more sensible for this unusual need. 165 </p></td></tr> 166</table></div> 167<p> 168 Full code of this example is at <a href="../../../../example/root_finding_n_example.cpp" target="_top">root_finding_n_example.cpp</a>. 169 </p> 170</div> 171<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 172<td align="left"></td> 173<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 174 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 175 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 176 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 177 Daryle Walker and Xiaogang Zhang<p> 178 Distributed under the Boost Software License, Version 1.0. (See accompanying 179 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>) 180 </p> 181</div></td> 182</tr></table> 183<hr> 184<div class="spirit-nav"> 185<a accesskey="p" href="multiprecision_root.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="elliptic_eg.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a> 186</div> 187</body> 188</html> 189
