<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Root Finding With Derivatives</title>
<link rel="stylesheet" href="../../../../../../../../doc/src/boostbook.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.74.0">
<link rel="home" href="../../../index.html" title="Math Toolkit">
<link rel="up" href="../internals1.html" title="Reused Utilities">
<link rel="prev" href="rational.html" title="Polynomial and Rational Function Evaluation">
<link rel="next" href="roots2.html" title="Root Finding Without Derivatives">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../../boost.png"></td>
<td align="center"><a href="../../../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="rational.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../internals1.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="roots2.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section" lang="en">
<div class="titlepage"><div><div><h4 class="title">
<a name="math_toolkit.toolkit.internals1.roots"></a><a class="link" href="roots.html" title="Root Finding With Derivatives"> Root Finding
With Derivatives</a>
</h4></div></div></div>
<a name="math_toolkit.toolkit.internals1.roots.synopsis"></a><h5>
<a name="id1205005"></a>
<a class="link" href="roots.html#math_toolkit.toolkit.internals1.roots.synopsis">Synopsis</a>
</h5>
<p>
</p>
<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>
</pre>
<p>
</p>
<pre class="programlisting"><span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span><span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">tools</span><span class="special">{</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span>
<span class="identifier">T</span> <span class="identifier">newton_raphson_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">);</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span>
<span class="identifier">T</span> <span class="identifier">newton_raphson_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</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_iter</span><span class="special">);</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span>
<span class="identifier">T</span> <span class="identifier">halley_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">);</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span>
<span class="identifier">T</span> <span class="identifier">halley_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</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_iter</span><span class="special">);</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span>
<span class="identifier">T</span> <span class="identifier">schroeder_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">);</span>
<span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span>
<span class="identifier">T</span> <span class="identifier">schroeder_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</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_iter</span><span class="special">);</span>
<span class="special">}}}</span> <span class="comment">// namespaces
</span></pre>
<a name="math_toolkit.toolkit.internals1.roots.description"></a><h5>
<a name="id1205872"></a>
<a class="link" href="roots.html#math_toolkit.toolkit.internals1.roots.description">Description</a>
</h5>
<p>
These functions all perform iterative root finding: <code class="computeroutput"><span class="identifier">newton_raphson_iterate</span></code>
performs second order <a class="link" href="roots.html#newton">Newton Raphson iteration</a>,
while <code class="computeroutput"><span class="identifier">halley_iterate</span></code> and
<code class="computeroutput"><span class="identifier">schroeder_iterate</span></code> perform
third order <a class="link" href="roots.html#halley">Halley</a> and <a class="link" href="roots.html#schroeder">Schroeder</a>
iteration respectively.
</p>
<p>
The functions all take the same parameters:
</p>
<div class="variablelist">
<p class="title"><b>Parameters of the root finding functions</b></p>
<dl>
<dt><span class="term">F f</span></dt>
<dd>
<p>
Type F must be a callable function object that accepts one parameter
and returns a tuple:
</p>
<p>
For the second order iterative methods (Newton Raphson) the tuple
should have two elements containing the evaluation of the function
and it's first derivative.
</p>
<p>
For the third order methods (Halley and Schroeder) the tuple should
have three elements containing the evaluation of the function and
its first and second derivatives.
</p>
</dd>
<dt><span class="term">T guess</span></dt>
<dd><p>
The initial starting value.
</p></dd>
<dt><span class="term">T min</span></dt>
<dd><p>
The minimum possible value for the result, this is used as an initial
lower bracket.
</p></dd>
<dt><span class="term">T max</span></dt>
<dd><p>
The maximum possible value for the result, this is used as an initial
upper bracket.
</p></dd>
<dt><span class="term">int digits</span></dt>
<dd><p>
The desired number of binary digits.
</p></dd>
<dt><span class="term">uintmax_t max_iter</span></dt>
<dd><p>
An optional maximum number of iterations to perform.
</p></dd>
</dl>
</div>
<p>
When using these functions you should note that:
</p>
<div class="itemizedlist"><ul type="disc">
<li>
They may be very sensitive to the initial guess, typically they converge
very rapidly if the initial guess has two or three decimal digits correct.
However convergence can be no better than bisection, or in some rare
cases even worse than bisection if the initial guess is a long way
from the correct value and the derivatives are close to zero.
</li>
<li>
These functions include special cases to handle zero first (and second
where appropriate) derivatives, and fall back to bisection in this
case. However, it is helpful if F is defined to return an arbitrarily
small value <span class="emphasis"><em>of the correct sign</em></span> rather than zero.
</li>
<li>
If the derivative at the current best guess for the result is infinite
(or very close to being infinite) then these functions may terminate
prematurely. A large first derivative leads to a very small next step,
triggering the termination condition. Derivative based iteration may
not be appropriate in such cases.
</li>
<li>
These functions fall back to bisection if the next computed step would
take the next value out of bounds. The bounds are updated after each
step to ensure this leads to convergence. However, a good initial guess
backed up by asymptotically-tight bounds will improve performance no
end rather than relying on bisection.
</li>
<li>
The value of <span class="emphasis"><em>digits</em></span> is crucial to good performance
of these functions, if it is set too high then at best you will get
one extra (unnecessary) iteration, and at worst the last few steps
will proceed by bisection. Remember that the returned value can never
be more accurate than f(x) can be evaluated, and that if f(x) suffers
from cancellation errors as it tends to zero then the computed steps
will be effectively random. The value of <span class="emphasis"><em>digits</em></span>
should be set so that iteration terminates before this point: remember
that for second and third order methods the number of correct digits
in the result is increasing quite substantially with each iteration,
<span class="emphasis"><em>digits</em></span> should be set by experiment so that the
final iteration just takes the next value into the zone where f(x)
becomes inaccurate.
</li>
<li>
Finally: you may well be able to do better than these functions by
hand-coding the heuristics used so that they are tailored to a specific
function. You may also be able to compute the ratio of derivatives
used by these methods more efficiently than computing the derivatives
themselves. As ever, algebraic simplification can be a big win.
</li>
</ul></div>
<a name="newton"></a><p>
</p>
<a name="math_toolkit.toolkit.internals1.roots.newton_raphson_method"></a><h5>
<a name="id1206381"></a>
<a class="link" href="roots.html#math_toolkit.toolkit.internals1.roots.newton_raphson_method">Newton
Raphson Method</a>
</h5>
<p>
Given an initial guess x0 the subsequent values are computed using:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/roots1.png"></span>
</p>
<p>
Out of bounds steps revert to bisection of the current bounds.
</p>
<p>
Under ideal conditions, the number of correct digits doubles with each
iteration.
</p>
<a name="halley"></a><p>
</p>
<a name="math_toolkit.toolkit.internals1.roots.halley_s_method"></a><h5>
<a name="id1206438"></a>
<a class="link" href="roots.html#math_toolkit.toolkit.internals1.roots.halley_s_method">Halley's
Method</a>
</h5>
<p>
Given an initial guess x0 the subsequent values are computed using:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/roots2.png"></span>
</p>
<p>
Over-compensation by the second derivative (one which would proceed in
the wrong direction) causes the method to revert to a Newton-Raphson step.
</p>
<p>
Out of bounds steps revert to bisection of the current bounds.
</p>
<p>
Under ideal conditions, the number of correct digits trebles with each
iteration.
</p>
<a name="schroeder"></a><p>
</p>
<a name="math_toolkit.toolkit.internals1.roots.schroeder_s_method"></a><h5>
<a name="id1206499"></a>
<a class="link" href="roots.html#math_toolkit.toolkit.internals1.roots.schroeder_s_method">Schroeder's
Method</a>
</h5>
<p>
Given an initial guess x0 the subsequent values are computed using:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/roots3.png"></span>
</p>
<p>
Over-compensation by the second derivative (one which would proceed in
the wrong direction) causes the method to revert to a Newton-Raphson step.
Likewise a Newton step is used whenever that Newton step would change the
next value by more than 10%.
</p>
<p>
Out of bounds steps revert to bisection of the current bounds.
</p>
<p>
Under ideal conditions, the number of correct digits trebles with each
iteration.
</p>
<a name="math_toolkit.toolkit.internals1.roots.example"></a><h5>
<a name="id1206553"></a>
<a class="link" href="roots.html#math_toolkit.toolkit.internals1.roots.example">Example</a>
</h5>
<p>
Lets suppose we want to find the cube root of a number, the equation we
want to solve along with its derivatives are:
</p>
<p>
<span class="inlinemediaobject"><img src="../../../../equations/roots4.png"></span>
</p>
<p>
To begin with lets solve the problem using Newton Raphson iterations, we'll
begin be defining a function object that returns the evaluation of the
function to solve, along with its first derivative:
</p>
<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>
<span class="keyword">struct</span> <span class="identifier">cbrt_functor</span>
<span class="special">{</span>
<span class="identifier">cbrt_functor</span><span class="special">(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">target</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">a</span><span class="special">(</span><span class="identifier">target</span><span class="special">){}</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">tr1</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="keyword">operator</span><span class="special">()(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">z</span><span class="special">)</span>
<span class="special">{</span>
<span class="identifier">T</span> <span class="identifier">sqr</span> <span class="special">=</span> <span class="identifier">z</span> <span class="special">*</span> <span class="identifier">z</span><span class="special">;</span>
<span class="keyword">return</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">tr1</span><span class="special">::</span><span class="identifier">make_tuple</span><span class="special">(</span><span class="identifier">sqr</span> <span class="special">*</span> <span class="identifier">z</span> <span class="special">-</span> <span class="identifier">a</span><span class="special">,</span> <span class="number">3</span> <span class="special">*</span> <span class="identifier">sqr</span><span class="special">);</span>
<span class="special">}</span>
<span class="keyword">private</span><span class="special">:</span>
<span class="identifier">T</span> <span class="identifier">a</span><span class="special">;</span>
<span class="special">};</span>
</pre>
<p>
Implementing the cube root is fairly trivial now, the hardest part is finding
a good approximation to begin with: in this case we'll just divide the
exponent by three:
</p>
<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>
<span class="identifier">T</span> <span class="identifier">cbrt</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">)</span>
<span class="special">{</span>
<span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">std</span><span class="special">;</span>
<span class="keyword">int</span> <span class="identifier">exp</span><span class="special">;</span>
<span class="identifier">frexp</span><span class="special">(</span><span class="identifier">z</span><span class="special">,</span> <span class="special">&</span><span class="identifier">exp</span><span class="special">);</span>
<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">exp</span><span class="special">/</span><span class="number">3</span><span class="special">);</span>
<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.0</span><span class="special">,</span> <span class="identifier">exp</span><span class="special">/</span><span class="number">3</span><span class="special">);</span>
<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.0</span><span class="special">,</span> <span class="identifier">exp</span><span class="special">/</span><span class="number">3</span><span class="special">);</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="keyword">return</span> <span class="identifier">tools</span><span class="special">::</span><span class="identifier">newton_raphson_iterate</span><span class="special">(</span><span class="identifier">detail</span><span class="special">::</span><span class="identifier">cbrt_functor</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="identifier">z</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="special">}</span>
</pre>
<p>
Using the test data in libs/math/test/cbrt_test.cpp this found the cube
root exact to the last digit in every case, and in no more than 6 iterations
at double precision. However, you will note that a high precision was used
in this example, exactly what was warned against earlier on in these docs!
In this particular case its possible to compute f(x) exactly and without
undue cancellation error, so a high limit is not too much of an issue.
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> gave
full precision in all but one of the test cases (and that one was out by
just one bit). The maximum number of iterations remained 6, but in most
cases was reduced by one.
</p>
<p>
Note also that the above code omits error handling, and does not handle
negative values of z correctly. That will be left as an exercise for the
reader!
</p>
<p>
Now lets adapt the functor slightly to return the second derivative as
well:
</p>
<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>
<span class="keyword">struct</span> <span class="identifier">cbrt_functor</span>
<span class="special">{</span>
<span class="identifier">cbrt_functor</span><span class="special">(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">target</span><span class="special">)</span> <span class="special">:</span> <span class="identifier">a</span><span class="special">(</span><span class="identifier">target</span><span class="special">){}</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">tr1</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">z</span><span class="special">)</span>
<span class="special">{</span>
<span class="identifier">T</span> <span class="identifier">sqr</span> <span class="special">=</span> <span class="identifier">z</span> <span class="special">*</span> <span class="identifier">z</span><span class="special">;</span>
<span class="keyword">return</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">tr1</span><span class="special">::</span><span class="identifier">make_tuple</span><span class="special">(</span><span class="identifier">sqr</span> <span class="special">*</span> <span class="identifier">z</span> <span class="special">-</span> <span class="identifier">a</span><span class="special">,</span> <span class="number">3</span> <span class="special">*</span> <span class="identifier">sqr</span><span class="special">,</span> <span class="number">6</span> <span class="special">*</span> <span class="identifier">z</span><span class="special">);</span>
<span class="special">}</span>
<span class="keyword">private</span><span class="special">:</span>
<span class="identifier">T</span> <span class="identifier">a</span><span class="special">;</span>
<span class="special">};</span>
</pre>
<p>
And then adapt the <code class="computeroutput"><span class="identifier">cbrt</span></code>
function to use Halley iterations:
</p>
<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>
<span class="identifier">T</span> <span class="identifier">cbrt</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">z</span><span class="special">)</span>
<span class="special">{</span>
<span class="keyword">using</span> <span class="keyword">namespace</span> <span class="identifier">std</span><span class="special">;</span>
<span class="keyword">int</span> <span class="identifier">exp</span><span class="special">;</span>
<span class="identifier">frexp</span><span class="special">(</span><span class="identifier">z</span><span class="special">,</span> <span class="special">&</span><span class="identifier">exp</span><span class="special">);</span>
<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">exp</span><span class="special">/</span><span class="number">3</span><span class="special">);</span>
<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.0</span><span class="special">,</span> <span class="identifier">exp</span><span class="special">/</span><span class="number">3</span><span class="special">);</span>
<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.0</span><span class="special">,</span> <span class="identifier">exp</span><span class="special">/</span><span class="number">3</span><span class="special">);</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="number">2</span><span class="special">;</span>
<span class="keyword">return</span> <span class="identifier">tools</span><span class="special">::</span><span class="identifier">halley_iterate</span><span class="special">(</span><span class="identifier">detail</span><span class="special">::</span><span class="identifier">cbrt_functor</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="identifier">z</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="special">}</span>
</pre>
<p>
Note that the iterations are set to stop at just one-half of full precision,
and yet even so not one of the test cases had a single bit wrong. What's
more, the maximum number of iterations was now just 4.
</p>
<p>
Just to complete the picture, we could have called <code class="computeroutput"><span class="identifier">schroeder_iterate</span></code>
in the last example: and in fact it makes no difference to the accuracy
or number of iterations in this particular case. However, the relative
performance of these two methods may vary depending upon the nature of
f(x), and the accuracy to which the initial guess can be computed. There
appear to be no generalisations that can be made except "try them
and see".
</p>
<p>
Finally, had we called cbrt with <a href="http://shoup.net/ntl/doc/RR.txt" target="_top">NTL::RR</a>
set to 1000 bit precision, then full precision can be obtained with just
7 iterations. To put that in perspective an increase in precision by a
factor of 20, has less than doubled the number of iterations. That just
goes to emphasise that most of the iterations are used up getting the first
few digits correct: after that these methods can churn out further digits
with remarkable efficiency. Or to put it another way: <span class="emphasis"><em>nothing
beats a really good initial guess!</em></span>
</p>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2006 , 2007, 2008, 2009 John Maddock, Paul A. Bristow,
Hubert Holin, Xiaogang Zhang, Bruno Lalande, Johan Råde, Gautam Sewani
and Thijs van den Berg<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
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>)
</p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="rational.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../internals1.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="roots2.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
