1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
|
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Root Finding With Derivatives: Newton-Raphson, Halley & Schroeder</title>
<link rel="stylesheet" href="../../math.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.77.1">
<link rel="home" href="../../index.html" title="Math Toolkit 2.2.0">
<link rel="up" href="../internals1.html" title="Utilities & internals">
<link rel="prev" href="rational.html" title="Polynomial and Rational Function Evaluation">
<link rel="next" href="roots2.html" title="Root Finding Without Derivatives: Bisection, Bracket and TOMS748">
</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">
<div class="titlepage"><div><div><h3 class="title">
<a name="math_toolkit.internals1.roots"></a><a class="link" href="roots.html" title="Root Finding With Derivatives: Newton-Raphson, Halley & Schroeder">Root Finding With Derivatives:
Newton-Raphson, Halley & Schroeder</a>
</h3></div></div></div>
<h5>
<a name="math_toolkit.internals1.roots.h0"></a>
<span class="phrase"><a name="math_toolkit.internals1.roots.synopsis"></a></span><a class="link" href="roots.html#math_toolkit.internals1.roots.synopsis">Synopsis</a>
</h5>
<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>
<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>
<h5>
<a name="math_toolkit.internals1.roots.h1"></a>
<span class="phrase"><a name="math_toolkit.internals1.roots.description"></a></span><a class="link" href="roots.html#math_toolkit.internals1.roots.description">Description</a>
</h5>
<p>
These functions all perform iterative root finding using derivatives:
</p>
<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
<li class="listitem">
<code class="computeroutput"><span class="identifier">newton_raphson_iterate</span></code>
performs second order <a class="link" href="roots.html#math_toolkit.internals1.roots.newton">Newton-Raphson
iteration</a>,
</li>
<li class="listitem">
<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#math_toolkit.internals1.roots.halley">Halley</a> and
<a class="link" href="roots.html#math_toolkit.internals1.roots.schroeder">Schroeder</a>
iteration.
</li>
</ul></div>
<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 class="variablelist">
<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 <a class="link" href="tuples.html" title="Tuples">boost::math::tuple</a>:
</p>
<p>
For the second order iterative methods (<a href="http://en.wikipedia.org/wiki/Newton_Raphson" target="_top">Newton
Raphson</a>) the <a class="link" href="tuples.html" title="Tuples">boost::math::tuple</a>
should have <span class="bold"><strong>two</strong></span> elements containing
the evaluation of the function and its first derivative.
</p>
<p>
For the third order methods (<a href="http://en.wikipedia.org/wiki/Halley%27s_method" target="_top">Halley</a>
and Schroeder) the <a class="link" href="tuples.html" title="Tuples">boost::math::tuple</a>
should have <span class="bold"><strong>three</strong></span> 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. A good guess is crucial to quick convergence!
</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. On exit this is
set to the actual number of iterations performed.
</p></dd>
</dl>
</div>
<p>
When using these functions you should note that:
</p>
<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
<li class="listitem">
Default max_iter = <code class="computeroutput"><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">>::</span><span class="identifier">max</span><span class="special">)()</span></code> is effectively 'iterate for ever'!.
</li>
<li class="listitem">
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 class="listitem">
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 functor F is defined to return an arbitrarily
small value <span class="emphasis"><em>of the correct sign</em></span> rather than zero.
</li>
<li class="listitem">
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 class="listitem">
If the function is 'Really Well Behaved' (monotonic and has only one
root) the bracket bounds min and max may as well be set to the widest
limits like zero and <code class="computeroutput"><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">max</span><span class="special">()</span></code>.
</li>
<li class="listitem">
But if the function more complex and may have more than one root or a
pole, the choice of bounds is protection against jumping out to seek
the 'wrong' root.
</li>
<li class="listitem">
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 class="listitem">
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 class="listitem">
To get the binary digits of accuracy, use policies::get_max_root_iterations<Policy>()).
</li>
<li class="listitem">
If you need some diagnostic output to see what is going on, you can
<code class="computeroutput"><span class="preprocessor">#define</span> <span class="identifier">BOOST_MATH_INSTRUMENT</span></code>
before the <code class="computeroutput"><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></code>, and also ensure that display of
all the possibly significant digits with <code class="computeroutput"> <span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</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="keyword">double</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">)</span></code>: but be warned, this may produce copious
output!
</li>
<li class="listitem">
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>
<h5>
<a name="math_toolkit.internals1.roots.h2"></a>
<span class="phrase"><a name="math_toolkit.internals1.roots.newton"></a></span><a class="link" href="roots.html#math_toolkit.internals1.roots.newton">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.svg"></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>
<h5>
<a name="math_toolkit.internals1.roots.h3"></a>
<span class="phrase"><a name="math_toolkit.internals1.roots.halley"></a></span><a class="link" href="roots.html#math_toolkit.internals1.roots.halley">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.svg"></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>
<h5>
<a name="math_toolkit.internals1.roots.h4"></a>
<span class="phrase"><a name="math_toolkit.internals1.roots.schroeder"></a></span><a class="link" href="roots.html#math_toolkit.internals1.roots.schroeder">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.svg"></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>
<h5>
<a name="math_toolkit.internals1.roots.h5"></a>
<span class="phrase"><a name="math_toolkit.internals1.roots.example"></a></span><a class="link" href="roots.html#math_toolkit.internals1.roots.example">Example</a>
</h5>
<p>
Let's 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.svg"></span>
</p>
<p>
To begin with lets solve the problem using Newton-Raphson iterations, we'll
begin by defining a function object (functor) that returns the evaluation
of the function to solve, along with its first derivative f'(x):
</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="special">{</span> <span class="comment">// Constructor stores value to be 'cube-rooted'.</span>
<span class="special">}</span>
<a class="link" href="tuples.html" title="Tuples">boost::math::tuple</a><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="comment">// z is estimate so far.</span>
<span class="keyword">return</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">make_tuple</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="identifier">z</span> <span class="special">-</span> <span class="identifier">a</span><span class="special">,</span> <span class="comment">// return both f(x)</span>
<span class="number">3</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="comment">// and f'(x)</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="comment">// to be 'cube-rooted'.</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="comment">// for frexp, ldexp, numeric_limits.</span>
<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="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="comment">// Get exponent of z (ignore mantissa).</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="comment">// Rough guess is to divide the exponent by three.</span>
<span class="keyword">int</span> <span class="identifier">digits</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">digits</span><span class="special">;</span> <span class="comment">// Maximum possible binary digits accuracy for type T.</span>
<span class="keyword">return</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 <code class="computeroutput"><span class="identifier">libs</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">test</span><span class="special">/</span><span class="identifier">cbrt_test</span><span class="special">.</span><span class="identifier">cpp</span></code> 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 it is 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 a probably optimization by computing
z², and reusing it, omits error handling, and does not handle negative values
of z correctly. (These are left as an exercise for the reader!)
</p>
<p>
The <code class="computeroutput"><span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cbrt</span></code> function also includes these and other
improvements.
</p>
<p>
Now let's 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>
<a class="link" href="tuples.html" title="Tuples">boost::math::tuple</a><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="keyword">return</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">make_tuple</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="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">z</span><span class="special">*</span><span class="identifier">z</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">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="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">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 <code class="computeroutput"><span class="identifier">cbrt</span></code>
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.
</p>
<p>
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-2010, 2012-2014 Nikhar Agrawal,
Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert
Holin, Bruno Lalande, John Maddock, Johan Råde, Gautam Sewani, Benjamin Sobotta,
Thijs van den Berg, Daryle Walker and Xiaogang Zhang<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>
|
