ceph/src/boost/libs/math/doc/html/math_toolkit/roots/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.77.1">
   7 <link rel="home" href="../../../index.html" title="Math Toolkit 2.5.1">
   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><h4 class="title">
  27 <a name="math_toolkit.roots.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 </h4></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 template parameter <code class="computeroutput"><span class="identifier">N</span></code> is an integer and a compile time constant.
  34           Our functor and function now have an additional template parameter <code class="computeroutput"><span class="identifier">N</span></code>, for the root required.
  35         </p>
  36 <div class="note"><table border="0" summary="Note">
  37 <tr>
  38 <td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../../../doc/src/images/note.png"></td>
  39 <th align="left">Note</th>
  40 </tr>
  41 <tr><td align="left" valign="top"><p>
  42             Since the powers and derivatives are fixed at compile time, the resulting
  43             code is as efficient as as if hand-coded as the cube and fifth-root examples
  44             above. A good compiler should also optimise any repeated multiplications.
  45           </p></td></tr>
  46 </table></div>
  47 <p>
  48           Our <span class="emphasis"><em>n</em></span>th root functor is
  49         </p>
  50 <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>
  51 <span class="keyword">struct</span> <span class="identifier">nth_functor_2deriv</span>
  52 <span class="special">{</span> <span class="comment">// Functor returning both 1st and 2nd derivatives.</span>
  53   <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>
  54   <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>
  55
  56   <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>
  57   <span class="special">{</span> <span class="comment">/* Constructor stores value a to find root of, for example: */</span> <span class="special">}</span>
  58
  59   <span class="comment">// using boost::math::tuple; // to return three values.</span>
  60   <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>
  61   <span class="special">{</span>
  62     <span class="comment">// Return f(x), f'(x) and f''(x).</span>
  63     <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>
  64     <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>
  65     <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>
  66     <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>
  67
  68     <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>
  69   <span class="special">}</span>
  70 <span class="keyword">private</span><span class="special">:</span>
  71   <span class="identifier">T</span> <span class="identifier">a</span><span class="special">;</span>                                     <span class="comment">// to be 'nth_rooted'.</span>
  72 <span class="special">};</span>
  73 </pre>
  74 <p>
  75           and our <span class="emphasis"><em>n</em></span>th root function is
  76         </p>
  77 <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>
  78 <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>
  79 <span class="special">{</span> <span class="comment">// return nth root of x using 1st and 2nd derivatives and Halley.</span>
  80
  81   <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>
  82   <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>
  83
  84   <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>
  85   <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>
  86   <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>
  87
  88   <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>
  89
  90   <span class="keyword">int</span> <span class="identifier">exponent</span><span class="special">;</span>
  91   <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>
  92   <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>
  93   <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>
  94   <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>
  95
  96   <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>
  97                                                                 <span class="comment">// slightly more than one third of the digits are correct.</span>
  98   <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>
  99   <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>
 100   <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>
 101   <span class="keyword">return</span> <span class="identifier">result</span><span class="special">;</span>
 102 <span class="special">}</span>
 103 </pre>
 104 <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>
 105     <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>
 106 <span class="preprocessor">#ifndef</span> <span class="identifier">_MSC_VER</span>  <span class="comment">// float128 is not supported by Microsoft compiler 2013.</span>
 107     <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>
 108 <span class="preprocessor">#endif</span>
 109     <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>
 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">cpp_bin_float_50</span><span class="special">&gt;(</span><span class="number">2</span><span class="special">);</span> <span class="comment">// bin</span>
 111 </pre>
 112 <p>
 113           produces an output similar to this
 114         </p>
 115 <pre class="programlisting"> <span class="identifier">Using</span> <span class="identifier">MSVC</span> <span class="number">2013</span>
 116
 117 <span class="identifier">nth</span> <span class="identifier">Root</span> <span class="identifier">finding</span> <span class="identifier">Example</span><span class="special">.</span>
 118 <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>
 119 <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>
 120 <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>
 121   <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>
 122 <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>
 123   <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>
 124 </pre>
 125 <div class="tip"><table border="0" summary="Tip">
 126 <tr>
 127 <td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../../../doc/src/images/tip.png"></td>
 128 <th align="left">Tip</th>
 129 </tr>
 130 <tr><td align="left" valign="top">
 131 <p>
 132             Take care with the type passed to the function. It is best to pass a
 133             <code class="computeroutput"><span class="keyword">double</span></code> or greater-precision
 134             floating-point type.
 135           </p>
 136 <p>
 137             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
 138             be 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>
 139             <span class="keyword">double</span><span class="special">&gt;(</span><span class="number">2</span><span class="special">)</span></code> converts
 140             the integer to <code class="computeroutput"><span class="keyword">double</span></code>.
 141           </p>
 142 <p>
 143             Avoid passing a <code class="computeroutput"><span class="keyword">float</span></code> value
 144             that will provoke warnings (actually spurious) from the compiler about
 145             potential loss of data, as noted above.
 146           </p>
 147 </td></tr>
 148 </table></div>
 149 <div class="warning"><table border="0" summary="Warning">
 150 <tr>
 151 <td rowspan="2" align="center" valign="top" width="25"><img alt="[Warning]" src="../../../../../../../doc/src/images/warning.png"></td>
 152 <th align="left">Warning</th>
 153 </tr>
 154 <tr><td align="left" valign="top"><p>
 155             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> may lead to <a href="http://en.wikipedia.org/wiki/Loss_of_significance" target="_top">Loss
 156             of significance</a> like <code class="computeroutput"><span class="identifier">Type</span>
 157             <span class="keyword">double</span> <span class="identifier">value</span>
 158             <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>
 159             <span class="special">=</span> <span class="number">1.00000069314783</span></code>.
 160             Use of the the <code class="computeroutput"><span class="identifier">pow</span></code> function
 161             is more sensible for this unusual need.
 162           </p></td></tr>
 163 </table></div>
 164 <p>
 165           Full code of this example is at <a href="../../../../../example/root_finding_n_example.cpp" target="_top">root_finding_n_example.cpp</a>.
 166         </p>
 167 </div>
 168 <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
 169 <td align="left"></td>
 170 <td align="right"><div class="copyright-footer">Copyright &#169; 2006-2010, 2012-2014 Nikhar Agrawal,
 171       Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert
 172       Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Johan R&#229;de, Gautam Sewani,
 173       Benjamin Sobotta, Thijs van den Berg, Daryle Walker and Xiaogang Zhang<p>
 174         Distributed under the Boost Software License, Version 1.0. (See accompanying
 175         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>)
 176       </p>
 177 </div></td>
 178 </tr></table>
 179 <hr>
 180 <div class="spirit-nav">
 181 <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>
 182 </div>
 183 </body>
 184 </html>