libs/math/doc/html/math_toolkit/root_finding_examples/nth_root.html

   1 <html>
   2 <head>
   3 <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
   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.11.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">&lt;</span><span class="identifier">N</span><span class="special">&gt;</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">&lt;</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">&gt;</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">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</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">&gt;</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 &gt; 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">&amp;</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">&lt;</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">&gt;</span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&amp;</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">&lt;</span><span class="identifier">N</span><span class="special">&gt;(</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">&lt;</span><span class="identifier">N</span> <span class="special">-</span> <span class="number">1</span><span class="special">&gt;(</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">&lt;</span><span class="identifier">N</span> <span class="special">-</span> <span class="number">2</span> <span class="special">&gt;(</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">&lt;</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">&gt;</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">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</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">&gt;</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 &gt; 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">&gt;</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">&lt;</span><span class="identifier">guess_type</span><span class="special">&gt;(</span><span class="identifier">x</span><span class="special">),</span> <span class="special">&amp;</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">&lt;</span><span class="identifier">guess_type</span><span class="special">&gt;(</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">&lt;</span><span class="identifier">guess_type</span><span class="special">&gt;(</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">&lt;</span><span class="identifier">guess_type</span><span class="special">&gt;(</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">&lt;</span><span class="identifier">T</span><span class="special">&gt;::</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">&lt;</span><span class="identifier">N</span><span class="special">,</span> <span class="identifier">T</span><span class="special">&gt;(</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">&lt;</span><span class="number">5</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;(</span><span class="number">2.</span><span class="special">);</span>
 108     <span class="identifier">show_nth_root</span><span class="special">&lt;</span><span class="number">5</span><span class="special">,</span> <span class="keyword">long</span> <span class="keyword">double</span><span class="special">&gt;(</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">&lt;</span><span class="number">5</span><span class="special">,</span> <span class="identifier">float128</span><span class="special">&gt;(</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">&lt;</span><span class="number">5</span><span class="special">,</span> <span class="identifier">cpp_dec_float_50</span><span class="special">&gt;(</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">&lt;</span><span class="number">5</span><span class="special">,</span> <span class="identifier">cpp_bin_float_50</span><span class="special">&gt;(</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">&lt;</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">&lt;</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">&gt;,</span><span class="number">1</span><span class="special">&gt;</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">&lt;</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">&lt;</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">&gt;,</span><span class="number">0</span><span class="special">&gt;</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">&lt;</span><span class="number">5</span><span class="special">&gt;(</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">&lt;</span><span class="number">5</span><span class="special">,</span>
 141           <span class="keyword">double</span><span class="special">&gt;(</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">&lt;</span><span class="number">1000000</span><span class="special">&gt;(</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 &#169; 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&#229;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>