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
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
|
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>Locating Function Minima using Brent's algorithm</title>
<link rel="stylesheet" href="../math.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
<link rel="home" href="../index.html" title="Math Toolkit 4.2.1">
<link rel="up" href="../root_finding.html" title="Chapter 10. Root Finding & Minimization Algorithms">
<link rel="prev" href="bad_roots.html" title="Examples Where Root Finding Goes Wrong">
<link rel="next" href="root_comparison.html" title="Comparison of Root Finding Algorithms">
<meta name="viewport" content="width=device-width, initial-scale=1">
</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="bad_roots.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../root_finding.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="root_comparison.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h2 class="title" style="clear: both">
<a name="math_toolkit.brent_minima"></a><a class="link" href="brent_minima.html" title="Locating Function Minima using Brent's algorithm">Locating Function Minima using
Brent's algorithm</a>
</h2></div></div></div>
<h5>
<a name="math_toolkit.brent_minima.h0"></a>
<span class="phrase"><a name="math_toolkit.brent_minima.synopsis"></a></span><a class="link" href="brent_minima.html#math_toolkit.brent_minima.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">minima</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
</pre>
<pre class="programlisting"><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">std</span><span class="special">::</span><span class="identifier">pair</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">brent_find_minima</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">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">bits</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">std</span><span class="special">::</span><span class="identifier">pair</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">brent_find_minima</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">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">bits</span><span class="special">,</span> <span class="identifier">std</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>
</pre>
<h5>
<a name="math_toolkit.brent_minima.h1"></a>
<span class="phrase"><a name="math_toolkit.brent_minima.description"></a></span><a class="link" href="brent_minima.html#math_toolkit.brent_minima.description">Description</a>
</h5>
<p>
These two functions locate the minima of the continuous function <span class="emphasis"><em>f</em></span>
using <a href="http://en.wikipedia.org/wiki/Brent%27s_method" target="_top">Brent's method</a>:
specifically it uses quadratic interpolation to locate the minima, or if that
fails, falls back to a <a href="http://en.wikipedia.org/wiki/Golden_section_search" target="_top">golden-section
search</a>.
</p>
<p>
<span class="bold"><strong>Parameters</strong></span>
</p>
<div class="variablelist">
<p class="title"><b></b></p>
<dl class="variablelist">
<dt><span class="term">f</span></dt>
<dd><p>
The function to minimise: a function object (or C++ lambda) that should
be smooth over the range <span class="emphasis"><em>[min, max]</em></span>, with no maxima
occurring in that interval.
</p></dd>
<dt><span class="term">min</span></dt>
<dd><p>
The lower endpoint of the range in which to search for the minima.
</p></dd>
<dt><span class="term">max</span></dt>
<dd><p>
The upper endpoint of the range in which to search for the minima.
</p></dd>
<dt><span class="term">bits</span></dt>
<dd><p>
The number of bits precision to which the minima should be found.<br>
Note that in principle, the minima can not be located to greater accuracy
than the square root of machine epsilon (for 64-bit double, sqrt(1e-16)≅1e-8),
therefore the value of <span class="emphasis"><em>bits</em></span> will be ignored if it's
greater than half the number of bits in the mantissa of T.
</p></dd>
<dt><span class="term">max_iter</span></dt>
<dd><p>
The maximum number of iterations to use in the algorithm, if not provided
the algorithm will just keep on going until the minima is found.
</p></dd>
</dl>
</div>
<p>
<span class="bold"><strong>Returns:</strong></span>
</p>
<p>
A <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span></code> of type T containing the value of the
abscissa (x) at the minima and the value of <span class="emphasis"><em>f(x)</em></span> at the
minima.
</p>
<div class="tip"><table border="0" summary="Tip">
<tr>
<td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../doc/src/images/tip.png"></td>
<th align="left">Tip</th>
</tr>
<tr><td align="left" valign="top">
<p>
Defining BOOST_MATH_INSTRUMENT will show some parameters, for example:
</p>
<pre class="programlisting"><span class="identifier">Type</span> <span class="identifier">T</span> <span class="identifier">is</span> <span class="keyword">double</span>
<span class="identifier">bits</span> <span class="special">=</span> <span class="number">24</span><span class="special">,</span> <span class="identifier">maximum</span> <span class="number">26</span>
<span class="identifier">tolerance</span> <span class="special">=</span> <span class="number">1.19209289550781e-007</span>
<span class="identifier">seeking</span> <span class="identifier">minimum</span> <span class="identifier">in</span> <span class="identifier">range</span> <span class="identifier">min</span><span class="special">-</span><span class="number">4</span> <span class="identifier">to</span> <span class="number">1.33333333333333</span>
<span class="identifier">maximum</span> <span class="identifier">iterations</span> <span class="number">18446744073709551615</span>
<span class="number">10</span> <span class="identifier">iterations</span><span class="special">.</span>
</pre>
</td></tr>
</table></div>
<h5>
<a name="math_toolkit.brent_minima.h2"></a>
<span class="phrase"><a name="math_toolkit.brent_minima.example"></a></span><a class="link" href="brent_minima.html#math_toolkit.brent_minima.example">Brent
Minimisation Example</a>
</h5>
<p>
As a demonstration, we replicate this <a href="http://en.wikipedia.org/wiki/Brent%27s_method#Example" target="_top">Wikipedia
example</a> minimising the function <span class="emphasis"><em>y= (x+3)(x-1)<sup>2</sup></em></span>.
</p>
<p>
It is obvious from the equation and the plot that there is a minimum at exactly
unity (x = 1) and the value of the function at one is exactly zero (y = 0).
</p>
<div class="tip"><table border="0" summary="Tip">
<tr>
<td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../doc/src/images/tip.png"></td>
<th align="left">Tip</th>
</tr>
<tr><td align="left" valign="top"><p>
This observation shows that an analytical or <a href="http://en.wikipedia.org/wiki/Closed-form_expression" target="_top">Closed-form
expression</a> solution always beats brute-force hands-down for both
speed and precision.
</p></td></tr>
</table></div>
<div class="blockquote"><blockquote class="blockquote"><p>
<span class="inlinemediaobject"><img src="../../graphs/brent_test_function_1.svg" align="middle"></span>
</p></blockquote></div>
<p>
First an include is needed:
</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">minima</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
</pre>
<p>
This function is encoded in C++ as function object (functor) using <code class="computeroutput"><span class="keyword">double</span></code> precision thus:
</p>
<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">funcdouble</span>
<span class="special">{</span>
<span class="keyword">double</span> <span class="keyword">operator</span><span class="special">()(</span><span class="keyword">double</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">x</span><span class="special">)</span>
<span class="special">{</span>
<span class="keyword">return</span> <span class="special">(</span><span class="identifier">x</span> <span class="special">+</span> <span class="number">3</span><span class="special">)</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">x</span> <span class="special">-</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">x</span> <span class="special">-</span> <span class="number">1</span><span class="special">);</span> <span class="comment">// (x + 3)(x - 1)^2</span>
<span class="special">}</span>
<span class="special">};</span>
</pre>
<p>
The Brent function is conveniently accessed through a <code class="computeroutput"><span class="keyword">using</span></code>
statement (noting sub-namespace <code class="computeroutput"><span class="special">::</span><span class="identifier">tools</span></code>).
</p>
<pre class="programlisting"><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">tools</span><span class="special">::</span><span class="identifier">brent_find_minima</span><span class="special">;</span>
</pre>
<p>
The search minimum and maximum are chosen as -4 to 4/3 (as in the Wikipedia
example).
</p>
<div class="tip"><table border="0" summary="Tip">
<tr>
<td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../doc/src/images/tip.png"></td>
<th align="left">Tip</th>
</tr>
<tr><td align="left" valign="top"><p>
S A Stage (reference 6) reports that the Brent algorithm is <span class="emphasis"><em>slow
to start, but fast to converge</em></span>, so choosing a tight min-max range
is good.
</p></td></tr>
</table></div>
<p>
For simplicity, we set the precision parameter <code class="computeroutput"><span class="identifier">bits</span></code>
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="keyword">double</span><span class="special">>::</span><span class="identifier">digits</span></code>, which is effectively the maximum
possible <code class="computeroutput"><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">digits</span><span class="special">/</span><span class="number">2</span></code>.
</p>
<p>
Nor do we provide a value for maximum iterations parameter <code class="computeroutput"><span class="identifier">max_iter</span></code>,
(probably unwisely), so the function will iterate until it finds a minimum.
</p>
<pre class="programlisting"><span class="keyword">const</span> <span class="keyword">int</span> <span class="identifier">double_bits</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">digits</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">brent_find_minima</span><span class="special">(</span><span class="identifier">funcdouble</span><span class="special">(),</span> <span class="special">-</span><span class="number">4.</span><span class="special">,</span> <span class="number">4.</span> <span class="special">/</span> <span class="number">3</span><span class="special">,</span> <span class="identifier">double_bits</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">streamsize</span> <span class="identifier">precision_1</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><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">digits10</span><span class="special">);</span>
<span class="comment">// Show all double precision decimal digits and trailing zeros.</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"x at minimum = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span>
<span class="special"><<</span> <span class="string">", f("</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
</pre>
<p>
The resulting <a href="http://en.cppreference.com/w/cpp/utility/pair" target="_top">std::pair</a>
contains the minimum close to one, and the minimum value close to zero.
</p>
<pre class="programlisting"><span class="identifier">x</span> <span class="identifier">at</span> <span class="identifier">minimum</span> <span class="special">=</span> <span class="number">1.00000000112345</span><span class="special">,</span>
<span class="identifier">f</span><span class="special">(</span><span class="number">1.00000000112345</span><span class="special">)</span> <span class="special">=</span> <span class="number">5.04852568272458e-018</span>
</pre>
<p>
The differences from the expected <span class="emphasis"><em>one</em></span> and <span class="emphasis"><em>zero</em></span>
are less than the uncertainty, for <code class="computeroutput"><span class="keyword">double</span></code>
1.5e-008, calculated from <code class="computeroutput"><span class="identifier">sqrt</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">epsilon</span><span class="special">())</span></code>.
</p>
<p>
We can use this uncertainty to check that the two values are close-enough to
those expected,
</p>
<pre class="programlisting"><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">fpc</span><span class="special">::</span><span class="identifier">close_at_tolerance</span><span class="special">;</span>
<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">fpc</span><span class="special">::</span><span class="identifier">is_small</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"x = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">", f(x) = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">boolalpha</span> <span class="special"><<</span> <span class="string">"x == 1 (compared to uncertainty "</span>
<span class="special"><<</span> <span class="identifier">uncertainty</span> <span class="special"><<</span> <span class="string">") is "</span> <span class="special"><<</span> <span class="identifier">is_close</span><span class="special">(</span><span class="number">1.</span><span class="special">,</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span><span class="special">,</span> <span class="identifier">uncertainty</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// true</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">boolalpha</span> <span class="special"><<</span> <span class="string">"f(x) == 0 (compared to uncertainty "</span>
<span class="special"><<</span> <span class="identifier">uncertainty</span> <span class="special"><<</span> <span class="string">") is "</span> <span class="special"><<</span> <span class="identifier">is_close</span><span class="special">(</span><span class="number">0.</span><span class="special">,</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span><span class="special">,</span> <span class="identifier">uncertainty</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> <span class="comment">// true</span>
</pre>
<p>
that outputs
</p>
<pre class="programlisting"><span class="identifier">x</span> <span class="special">==</span> <span class="number">1</span> <span class="special">(</span><span class="identifier">compared</span> <span class="identifier">to</span> <span class="identifier">uncertainty</span> <span class="number">1.49011611938477e-08</span><span class="special">)</span> <span class="identifier">is</span> <span class="keyword">true</span>
<span class="identifier">f</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special">==</span> <span class="number">0</span> <span class="special">(</span><span class="identifier">compared</span> <span class="identifier">to</span> <span class="identifier">uncertainty</span> <span class="number">1.49011611938477e-08</span><span class="special">)</span> <span class="identifier">is</span> <span class="keyword">true</span>
</pre>
<div class="note"><table border="0" summary="Note">
<tr>
<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../doc/src/images/note.png"></td>
<th align="left">Note</th>
</tr>
<tr><td align="left" valign="top"><p>
The function <code class="computeroutput"><span class="identifier">close_at_tolerance</span></code>
is ineffective for testing if a value is small or zero (because it may divide
by small or zero and cause overflow). Function <code class="computeroutput"><span class="identifier">is_small</span></code>
does this job.
</p></td></tr>
</table></div>
<p>
It is possible to make this comparison more generally with a templated function,
<code class="computeroutput"><span class="identifier">is_close</span></code>, first checking <code class="computeroutput"><span class="identifier">is_small</span></code> before checking <code class="computeroutput"><span class="identifier">close_at_tolerance</span></code>,
returning <code class="computeroutput"><span class="keyword">true</span></code> when small or close,
for example:
</p>
<pre class="programlisting"><span class="comment">//! Compare if value got is close to expected,</span>
<span class="comment">//! checking first if expected is very small</span>
<span class="comment">//! (to avoid divide by tiny or zero during comparison)</span>
<span class="comment">//! before comparing expect with value got.</span>
<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">bool</span> <span class="identifier">is_close</span><span class="special">(</span><span class="identifier">T</span> <span class="identifier">expect</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">got</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">tolerance</span><span class="special">)</span>
<span class="special">{</span>
<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">fpc</span><span class="special">::</span><span class="identifier">close_at_tolerance</span><span class="special">;</span>
<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">fpc</span><span class="special">::</span><span class="identifier">is_small</span><span class="special">;</span>
<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">fpc</span><span class="special">::</span><span class="identifier">FPC_STRONG</span><span class="special">;</span>
<span class="keyword">if</span> <span class="special">(</span><span class="identifier">is_small</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="identifier">expect</span><span class="special">,</span> <span class="identifier">tolerance</span><span class="special">))</span>
<span class="special">{</span>
<span class="keyword">return</span> <span class="identifier">is_small</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="identifier">got</span><span class="special">,</span> <span class="identifier">tolerance</span><span class="special">);</span>
<span class="special">}</span>
<span class="keyword">return</span> <span class="identifier">close_at_tolerance</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="identifier">tolerance</span><span class="special">,</span> <span class="identifier">FPC_STRONG</span><span class="special">)</span> <span class="special">(</span><span class="identifier">expect</span><span class="special">,</span> <span class="identifier">got</span><span class="special">);</span>
<span class="special">}</span> <span class="comment">// bool is_close(T expect, T got, T tolerance)</span>
</pre>
<p>
In practical applications, we might want to know how many iterations, and maybe
to limit iterations (in case the function proves ill-behaved), and/or perhaps
to trade some loss of precision for speed, for example:
</p>
<pre class="programlisting"><span class="keyword">const</span> <span class="identifier">std</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>
<span class="identifier">std</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>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">brent_find_minima</span><span class="special">(</span><span class="identifier">funcdouble</span><span class="special">(),</span> <span class="special">-</span><span class="number">4.</span><span class="special">,</span> <span class="number">4.</span> <span class="special">/</span> <span class="number">3</span><span class="special">,</span> <span class="identifier">bits</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"x at minimum = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">", f("</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span>
<span class="special"><<</span> <span class="string">" after "</span> <span class="special"><<</span> <span class="identifier">it</span> <span class="special"><<</span> <span class="string">" iterations. "</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
</pre>
<p>
limits to a maximum of 20 iterations (a reasonable estimate for this example
function, even for much higher precision shown later).
</p>
<p>
The parameter <code class="computeroutput"><span class="identifier">it</span></code> is updated
to return the actual number of iterations (so it may be useful to also keep
a record of the limit set in <code class="computeroutput"><span class="keyword">const</span> <span class="identifier">maxit</span></code>).
</p>
<p>
It is neat to avoid showing insignificant digits by computing the number of
decimal digits to display.
</p>
<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">streamsize</span> <span class="identifier">prec</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">int</span><span class="special">>(</span><span class="number">2</span> <span class="special">+</span> <span class="identifier">sqrt</span><span class="special">((</span><span class="keyword">double</span><span class="special">)</span><span class="identifier">bits</span><span class="special">));</span> <span class="comment">// Number of significant decimal digits.</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">streamsize</span> <span class="identifier">precision_3</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">prec</span><span class="special">);</span> <span class="comment">// Save and set new precision.</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Showing "</span> <span class="special"><<</span> <span class="identifier">bits</span> <span class="special"><<</span> <span class="string">" bits "</span>
<span class="string">"precision with "</span> <span class="special"><<</span> <span class="identifier">prec</span>
<span class="special"><<</span> <span class="string">" decimal digits from tolerance "</span> <span class="special"><<</span> <span class="identifier">sqrt</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">epsilon</span><span class="special">())</span>
<span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"x at minimum = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span>
<span class="special"><<</span> <span class="string">", f("</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span>
<span class="special"><<</span> <span class="string">" after "</span> <span class="special"><<</span> <span class="identifier">it</span> <span class="special"><<</span> <span class="string">" iterations. "</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">precision_3</span><span class="special">);</span> <span class="comment">// Restore.</span>
</pre>
<pre class="programlisting"><span class="identifier">Showing</span> <span class="number">53</span> <span class="identifier">bits</span> <span class="identifier">precision</span> <span class="identifier">with</span> <span class="number">9</span> <span class="identifier">decimal</span> <span class="identifier">digits</span> <span class="identifier">from</span> <span class="identifier">tolerance</span> <span class="number">1.49011611938477e-008</span>
<span class="identifier">x</span> <span class="identifier">at</span> <span class="identifier">minimum</span> <span class="special">=</span> <span class="number">1</span><span class="special">,</span> <span class="identifier">f</span><span class="special">(</span><span class="number">1</span><span class="special">)</span> <span class="special">=</span> <span class="number">5.04852568e-018</span>
</pre>
<p>
We can also half the number of precision bits from 52 to 26:
</p>
<pre class="programlisting"><span class="keyword">const</span> <span class="keyword">int</span> <span class="identifier">bits_div_2</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">digits</span> <span class="special">/</span> <span class="number">2</span><span class="special">;</span> <span class="comment">// Half digits precision (effective maximum).</span>
<span class="keyword">double</span> <span class="identifier">epsilon_2</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">pow</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">digits</span><span class="special">/</span><span class="number">2</span> <span class="special">-</span> <span class="number">1</span><span class="special">),</span> <span class="keyword">double</span><span class="special">>(</span><span class="number">2</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">streamsize</span> <span class="identifier">prec</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">int</span><span class="special">>(</span><span class="number">2</span> <span class="special">+</span> <span class="identifier">sqrt</span><span class="special">((</span><span class="keyword">double</span><span class="special">)</span><span class="identifier">bits_div_2</span><span class="special">));</span> <span class="comment">// Number of significant decimal digits.</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Showing "</span> <span class="special"><<</span> <span class="identifier">bits_div_2</span> <span class="special"><<</span> <span class="string">" bits precision with "</span> <span class="special"><<</span> <span class="identifier">prec</span>
<span class="special"><<</span> <span class="string">" decimal digits from tolerance "</span> <span class="special"><<</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">epsilon_2</span><span class="special">)</span>
<span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">streamsize</span> <span class="identifier">precision_4</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">prec</span><span class="special">);</span> <span class="comment">// Save.</span>
<span class="keyword">const</span> <span class="identifier">std</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>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it_4</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">brent_find_minima</span><span class="special">(</span><span class="identifier">funcdouble</span><span class="special">(),</span> <span class="special">-</span><span class="number">4.</span><span class="special">,</span> <span class="number">4.</span> <span class="special">/</span> <span class="number">3</span><span class="special">,</span> <span class="identifier">bits_div_2</span><span class="special">,</span> <span class="identifier">it_4</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"x at minimum = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">", f("</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">it_4</span> <span class="special"><<</span> <span class="string">" iterations. "</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">precision_4</span><span class="special">);</span> <span class="comment">// Restore.</span>
</pre>
<p>
showing no change in the result and no change in the number of iterations,
as expected.
</p>
<p>
It is only if we reduce the precision to a quarter, specifying only 13 precision
bits
</p>
<pre class="programlisting"><span class="keyword">const</span> <span class="keyword">int</span> <span class="identifier">bits_div_4</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">digits</span> <span class="special">/</span> <span class="number">4</span><span class="special">;</span> <span class="comment">// Quarter precision.</span>
<span class="keyword">double</span> <span class="identifier">epsilon_4</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">pow</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">digits</span> <span class="special">/</span> <span class="number">4</span> <span class="special">-</span> <span class="number">1</span><span class="special">),</span> <span class="keyword">double</span><span class="special">>(</span><span class="number">2</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">streamsize</span> <span class="identifier">prec</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">int</span><span class="special">>(</span><span class="number">2</span> <span class="special">+</span> <span class="identifier">sqrt</span><span class="special">((</span><span class="keyword">double</span><span class="special">)</span><span class="identifier">bits_div_4</span><span class="special">));</span> <span class="comment">// Number of significant decimal digits.</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Showing "</span> <span class="special"><<</span> <span class="identifier">bits_div_4</span> <span class="special"><<</span> <span class="string">" bits precision with "</span> <span class="special"><<</span> <span class="identifier">prec</span>
<span class="special"><<</span> <span class="string">" decimal digits from tolerance "</span> <span class="special"><<</span> <span class="identifier">sqrt</span><span class="special">(</span><span class="identifier">epsilon_4</span><span class="special">)</span>
<span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">streamsize</span> <span class="identifier">precision_5</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">prec</span><span class="special">);</span> <span class="comment">// Save & set.</span>
<span class="keyword">const</span> <span class="identifier">std</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>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it_5</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">brent_find_minima</span><span class="special">(</span><span class="identifier">funcdouble</span><span class="special">(),</span> <span class="special">-</span><span class="number">4.</span><span class="special">,</span> <span class="number">4.</span> <span class="special">/</span> <span class="number">3</span><span class="special">,</span> <span class="identifier">bits_div_4</span><span class="special">,</span> <span class="identifier">it_5</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"x at minimum = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">", f("</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span>
<span class="special"><<</span> <span class="string">", after "</span> <span class="special"><<</span> <span class="identifier">it_5</span> <span class="special"><<</span> <span class="string">" iterations. "</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">precision_5</span><span class="special">);</span> <span class="comment">// Restore.</span>
</pre>
<p>
that we reduce the number of iterations from 10 to 7 that the result slightly
differs from <span class="emphasis"><em>one</em></span> and <span class="emphasis"><em>zero</em></span>.
</p>
<pre class="programlisting"><span class="identifier">Showing</span> <span class="number">13</span> <span class="identifier">bits</span> <span class="identifier">precision</span> <span class="identifier">with</span> <span class="number">9</span> <span class="identifier">decimal</span> <span class="identifier">digits</span> <span class="identifier">from</span> <span class="identifier">tolerance</span> <span class="number">0.015625</span>
<span class="identifier">x</span> <span class="identifier">at</span> <span class="identifier">minimum</span> <span class="special">=</span> <span class="number">0.9999776</span><span class="special">,</span> <span class="identifier">f</span><span class="special">(</span><span class="number">0.9999776</span><span class="special">)</span> <span class="special">=</span> <span class="number">2.0069572e-009</span> <span class="identifier">after</span> <span class="number">7</span> <span class="identifier">iterations</span><span class="special">.</span>
</pre>
<h6>
<a name="math_toolkit.brent_minima.h3"></a>
<span class="phrase"><a name="math_toolkit.brent_minima.template"></a></span><a class="link" href="brent_minima.html#math_toolkit.brent_minima.template">Templating
on floating-point type</a>
</h6>
<p>
If we want to switch the floating-point type, then the functor must be revised.
Since the functor is stateless, the easiest option is to simply make <code class="computeroutput"><span class="keyword">operator</span><span class="special">()</span></code> a
template member function:
</p>
<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">func</span>
<span class="special">{</span>
<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="keyword">operator</span><span class="special">()(</span><span class="identifier">T</span> <span class="keyword">const</span><span class="special">&</span> <span class="identifier">x</span><span class="special">)</span>
<span class="special">{</span>
<span class="keyword">return</span> <span class="special">(</span><span class="identifier">x</span> <span class="special">+</span> <span class="number">3</span><span class="special">)</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">x</span> <span class="special">-</span> <span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="special">(</span><span class="identifier">x</span> <span class="special">-</span> <span class="number">1</span><span class="special">);</span> <span class="comment">// (x + 3)(x - 1)^2</span>
<span class="special">}</span>
<span class="special">};</span>
</pre>
<p>
The <code class="computeroutput"><span class="identifier">brent_find_minima</span></code> function
can now be used in template form, for example using <code class="computeroutput"><span class="keyword">long</span>
<span class="keyword">double</span></code>:
</p>
<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">streamsize</span> <span class="identifier">precision_t1</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><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">long</span> <span class="keyword">double</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span> <span class="comment">// Save & set.</span>
<span class="keyword">long</span> <span class="keyword">double</span> <span class="identifier">bracket_min</span> <span class="special">=</span> <span class="special">-</span><span class="number">4.</span><span class="special">;</span>
<span class="keyword">long</span> <span class="keyword">double</span> <span class="identifier">bracket_max</span> <span class="special">=</span> <span class="number">4.</span> <span class="special">/</span> <span class="number">3</span><span class="special">;</span>
<span class="keyword">const</span> <span class="keyword">int</span> <span class="identifier">bits</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">long</span> <span class="keyword">double</span><span class="special">>::</span><span class="identifier">digits</span><span class="special">;</span>
<span class="keyword">const</span> <span class="identifier">std</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>
<span class="identifier">std</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>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="keyword">long</span> <span class="keyword">double</span><span class="special">,</span> <span class="keyword">long</span> <span class="keyword">double</span><span class="special">></span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">brent_find_minima</span><span class="special">(</span><span class="identifier">func</span><span class="special">(),</span> <span class="identifier">bracket_min</span><span class="special">,</span> <span class="identifier">bracket_max</span><span class="special">,</span> <span class="identifier">bits</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"x at minimum = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">", f("</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span>
<span class="special"><<</span> <span class="string">", after "</span> <span class="special"><<</span> <span class="identifier">it</span> <span class="special"><<</span> <span class="string">" iterations. "</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">precision_t1</span><span class="special">);</span> <span class="comment">// Restore.</span>
</pre>
<p>
The form shown uses the floating-point type <code class="computeroutput"><span class="keyword">long</span>
<span class="keyword">double</span></code> by deduction, but it is also
possible to be more explicit, for example:
</p>
<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="keyword">long</span> <span class="keyword">double</span><span class="special">,</span> <span class="keyword">long</span> <span class="keyword">double</span><span class="special">></span> <span class="identifier">r</span> <span class="special">=</span> <span class="identifier">brent_find_minima</span><span class="special"><</span><span class="identifier">func</span><span class="special">,</span> <span class="keyword">long</span> <span class="keyword">double</span><span class="special">></span>
<span class="special">(</span><span class="identifier">func</span><span class="special">(),</span> <span class="identifier">bracket_min</span><span class="special">,</span> <span class="identifier">bracket_max</span><span class="special">,</span> <span class="identifier">bits</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span>
</pre>
<p>
In order to show the use of multiprecision below (or other user-defined types),
it may be convenient to write a templated function to use this:
</p>
<pre class="programlisting"><span class="comment">//! Example template function to find and show minima.</span>
<span class="comment">//! \tparam T floating-point or fixed_point type.</span>
<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">void</span> <span class="identifier">show_minima</span><span class="special">()</span>
<span class="special">{</span>
<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">tools</span><span class="special">::</span><span class="identifier">brent_find_minima</span><span class="special">;</span>
<span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">sqrt</span><span class="special">;</span>
<span class="keyword">try</span>
<span class="special">{</span> <span class="comment">// Always use try'n'catch blocks with Boost.Math to ensure you get any error messages.</span>
<span class="keyword">int</span> <span class="identifier">bits</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="comment">// Maximum is digits/2;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">streamsize</span> <span class="identifier">prec</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="keyword">int</span><span class="special">>(</span><span class="number">2</span> <span class="special">+</span> <span class="identifier">sqrt</span><span class="special">((</span><span class="keyword">double</span><span class="special">)</span><span class="identifier">bits</span><span class="special">));</span> <span class="comment">// Number of significant decimal digits.</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">streamsize</span> <span class="identifier">precision</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">prec</span><span class="special">);</span> <span class="comment">// Save and set.</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"\n\nFor type: "</span> <span class="special"><<</span> <span class="keyword">typeid</span><span class="special">(</span><span class="identifier">T</span><span class="special">).</span><span class="identifier">name</span><span class="special">()</span>
<span class="special"><<</span> <span class="string">",\n epsilon = "</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">epsilon</span><span class="special">()</span>
<span class="comment">// << ", precision of " << bits << " bits"</span>
<span class="special"><<</span> <span class="string">",\n the maximum theoretical precision from Brent's minimization is "</span>
<span class="special"><<</span> <span class="identifier">sqrt</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">epsilon</span><span class="special">())</span>
<span class="special"><<</span> <span class="string">"\n Displaying to std::numeric_limits<T>::digits10 "</span> <span class="special"><<</span> <span class="identifier">prec</span> <span class="special"><<</span> <span class="string">", significant decimal digits."</span>
<span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="keyword">const</span> <span class="identifier">std</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>
<span class="identifier">std</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>
<span class="comment">// Construct using string, not double, avoids loss of precision.</span>
<span class="comment">//T bracket_min = static_cast<T>("-4");</span>
<span class="comment">//T bracket_max = static_cast<T>("1.3333333333333333333333333333333333333333333333333");</span>
<span class="comment">// Construction from double may cause loss of precision for multiprecision types like cpp_bin_float,</span>
<span class="comment">// but brackets values are good enough for using Brent minimization.</span>
<span class="identifier">T</span> <span class="identifier">bracket_min</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(-</span><span class="number">4</span><span class="special">);</span>
<span class="identifier">T</span> <span class="identifier">bracket_max</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="number">1.3333333333333333333333333333333333333333333333333</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</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">r</span> <span class="special">=</span> <span class="identifier">brent_find_minima</span><span class="special"><</span><span class="identifier">func</span><span class="special">,</span> <span class="identifier">T</span><span class="special">>(</span><span class="identifier">func</span><span class="special">(),</span> <span class="identifier">bracket_min</span><span class="special">,</span> <span class="identifier">bracket_max</span><span class="special">,</span> <span class="identifier">bits</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">" x at minimum = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">", f("</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span><span class="special">;</span>
<span class="keyword">if</span> <span class="special">(</span><span class="identifier">it</span> <span class="special"><</span> <span class="identifier">maxit</span><span class="special">)</span>
<span class="special">{</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">",\n met "</span> <span class="special"><<</span> <span class="identifier">bits</span> <span class="special"><<</span> <span class="string">" bits precision"</span> <span class="special"><<</span> <span class="string">", after "</span> <span class="special"><<</span> <span class="identifier">it</span> <span class="special"><<</span> <span class="string">" iterations."</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="special">}</span>
<span class="keyword">else</span>
<span class="special">{</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">",\n did NOT meet "</span> <span class="special"><<</span> <span class="identifier">bits</span> <span class="special"><<</span> <span class="string">" bits precision"</span> <span class="special"><<</span> <span class="string">" after "</span> <span class="special"><<</span> <span class="identifier">it</span> <span class="special"><<</span> <span class="string">" iterations!"</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="special">}</span>
<span class="comment">// Check that result is that expected (compared to theoretical uncertainty).</span>
<span class="identifier">T</span> <span class="identifier">uncertainty</span> <span class="special">=</span> <span class="identifier">sqrt</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">epsilon</span><span class="special">());</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">boolalpha</span> <span class="special"><<</span> <span class="string">"x == 1 (compared to uncertainty "</span> <span class="special"><<</span> <span class="identifier">uncertainty</span> <span class="special"><<</span> <span class="string">") is "</span>
<span class="special"><<</span> <span class="identifier">is_close</span><span class="special">(</span><span class="keyword">static_cast</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="number">1</span><span class="special">),</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span><span class="special">,</span> <span class="identifier">uncertainty</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">boolalpha</span> <span class="special"><<</span> <span class="string">"f(x) == (0 compared to uncertainty "</span> <span class="special"><<</span> <span class="identifier">uncertainty</span> <span class="special"><<</span> <span class="string">") is "</span>
<span class="special"><<</span> <span class="identifier">is_close</span><span class="special">(</span><span class="keyword">static_cast</span><span class="special"><</span><span class="identifier">T</span><span class="special">>(</span><span class="number">0</span><span class="special">),</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span><span class="special">,</span> <span class="identifier">uncertainty</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="comment">// Problems with this using multiprecision with expression template on?</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">precision</span><span class="special">);</span> <span class="comment">// Restore.</span>
<span class="special">}</span>
<span class="keyword">catch</span> <span class="special">(</span><span class="keyword">const</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">exception</span><span class="special">&</span> <span class="identifier">e</span><span class="special">)</span>
<span class="special">{</span> <span class="comment">// Always useful to include try & catch blocks because default policies</span>
<span class="comment">// are to throw exceptions on arguments that cause errors like underflow, overflow.</span>
<span class="comment">// Lacking try & catch blocks, the program will abort without a message below,</span>
<span class="comment">// which may give some helpful clues as to the cause of the exception.</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span>
<span class="string">"\n"</span><span class="string">"Message from thrown exception was:\n "</span> <span class="special"><<</span> <span class="identifier">e</span><span class="special">.</span><span class="identifier">what</span><span class="special">()</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="special">}</span>
<span class="special">}</span> <span class="comment">// void show_minima()</span>
</pre>
<div class="note"><table border="0" summary="Note">
<tr>
<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../doc/src/images/note.png"></td>
<th align="left">Note</th>
</tr>
<tr><td align="left" valign="top"><p>
the prudent addition of <code class="computeroutput"><span class="keyword">try</span><span class="char">'n'</span><span class="keyword">catch</span></code> blocks
within the function. This ensures that any malfunction will provide a Boost.Math
error message rather than uncommunicatively calling <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">abort</span></code>.
</p></td></tr>
</table></div>
<p>
We can use this with all built-in floating-point types, for example
</p>
<pre class="programlisting"><span class="identifier">show_minima</span><span class="special"><</span><span class="keyword">float</span><span class="special">>();</span>
<span class="identifier">show_minima</span><span class="special"><</span><span class="keyword">double</span><span class="special">>();</span>
<span class="identifier">show_minima</span><span class="special"><</span><span class="keyword">long</span> <span class="keyword">double</span><span class="special">>();</span>
</pre>
<h6>
<a name="math_toolkit.brent_minima.h4"></a>
<span class="phrase"><a name="math_toolkit.brent_minima.quad_precision"></a></span><a class="link" href="brent_minima.html#math_toolkit.brent_minima.quad_precision">Quad
128-bit precision</a>
</h6>
<p>
On platforms that provide it, a <a href="http://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format" target="_top">128-bit
quad</a> type can be used. (See <a href="../../../../../libs/multiprecision/doc/html/boost_multiprecision/tut/floats/float128.html" target="_top">float128</a>).
</p>
<p>
If this type is available, the build should define the macro BOOST_HAVE_QUADMATH:
</p>
<pre class="programlisting"><span class="preprocessor">#ifdef</span> <span class="identifier">BOOST_HAVE_QUADMATH</span> <span class="comment">// Defined only if GCC or Intel and have quadmath.lib or .dll library available.</span>
<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">float128</span><span class="special">;</span>
<span class="preprocessor">#endif</span>
</pre>
<p>
and can be (conditionally) used this:
</p>
<pre class="programlisting"><span class="preprocessor">#ifdef</span> <span class="identifier">BOOST_HAVE_QUADMATH</span> <span class="comment">// Defined only if GCC or Intel and have quadmath.lib or .dll library available.</span>
<span class="identifier">show_minima</span><span class="special"><</span><span class="identifier">float128</span><span class="special">>();</span> <span class="comment">// Needs quadmath_snprintf, sqrtQ, fabsq that are in in quadmath library.</span>
<span class="preprocessor">#endif</span>
</pre>
<h6>
<a name="math_toolkit.brent_minima.h5"></a>
<span class="phrase"><a name="math_toolkit.brent_minima.multiprecision"></a></span><a class="link" href="brent_minima.html#math_toolkit.brent_minima.multiprecision">Multiprecision</a>
</h6>
<p>
If a higher precision than <code class="computeroutput"><span class="keyword">double</span></code>
(or <code class="computeroutput"><span class="keyword">long</span> <span class="keyword">double</span></code>
if that is more precise) is required, then this is easily achieved using <a href="../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a>
with some includes, for example:
</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">multiprecision</span><span class="special">/</span><span class="identifier">cpp_dec_float</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> <span class="comment">// For decimal boost::multiprecision::cpp_dec_float_50.</span>
<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">multiprecision</span><span class="special">/</span><span class="identifier">cpp_bin_float</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> <span class="comment">// For binary boost::multiprecision::cpp_bin_float_50;</span>
</pre>
<p>
and some <code class="computeroutput"><span class="identifier">typdef</span></code>s.
</p>
<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_bin_float_50</span><span class="special">;</span> <span class="comment">// binary multiprecision typedef.</span>
<span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_dec_float_50</span><span class="special">;</span> <span class="comment">// decimal multiprecision typedef.</span>
<span class="comment">// One might also need typedefs like these to switch expression templates off and on (default is on).</span>
<span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">number</span><span class="special"><</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_bin_float</span><span class="special"><</span><span class="number">50</span><span class="special">>,</span>
<span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">et_on</span><span class="special">></span>
<span class="identifier">cpp_bin_float_50_et_on</span><span class="special">;</span> <span class="comment">// et_on is default so is same as cpp_bin_float_50.</span>
<span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">number</span><span class="special"><</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_bin_float</span><span class="special"><</span><span class="number">50</span><span class="special">>,</span>
<span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">et_off</span><span class="special">></span>
<span class="identifier">cpp_bin_float_50_et_off</span><span class="special">;</span>
<span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">number</span><span class="special"><</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_dec_float</span><span class="special"><</span><span class="number">50</span><span class="special">>,</span>
<span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">et_on</span><span class="special">></span> <span class="comment">// et_on is default so is same as cpp_dec_float_50.</span>
<span class="identifier">cpp_dec_float_50_et_on</span><span class="special">;</span>
<span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">number</span><span class="special"><</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">cpp_dec_float</span><span class="special"><</span><span class="number">50</span><span class="special">>,</span>
<span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">et_off</span><span class="special">></span>
<span class="identifier">cpp_dec_float_50_et_off</span><span class="special">;</span>
</pre>
<p>
Used thus
</p>
<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><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="identifier">cpp_bin_float_50</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span>
<span class="keyword">int</span> <span class="identifier">bits</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">cpp_bin_float_50</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">2</span><span class="special">;</span>
<span class="identifier">cpp_bin_float_50</span> <span class="identifier">bracket_min</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="identifier">cpp_bin_float_50</span><span class="special">>(</span><span class="string">"-4"</span><span class="special">);</span>
<span class="identifier">cpp_bin_float_50</span> <span class="identifier">bracket_max</span> <span class="special">=</span> <span class="keyword">static_cast</span><span class="special"><</span><span class="identifier">cpp_bin_float_50</span><span class="special">>(</span><span class="string">"1.3333333333333333333333333333333333333333333333333"</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Bracketing "</span> <span class="special"><<</span> <span class="identifier">bracket_min</span> <span class="special"><<</span> <span class="string">" to "</span> <span class="special"><<</span> <span class="identifier">bracket_max</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="keyword">const</span> <span class="identifier">std</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>
<span class="identifier">std</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> <span class="comment">// Will be updated with actual iteration count.</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="identifier">cpp_bin_float_50</span><span class="special">,</span> <span class="identifier">cpp_bin_float_50</span><span class="special">></span> <span class="identifier">r</span>
<span class="special">=</span> <span class="identifier">brent_find_minima</span><span class="special">(</span><span class="identifier">func</span><span class="special">(),</span> <span class="identifier">bracket_min</span><span class="special">,</span> <span class="identifier">bracket_max</span><span class="special">,</span> <span class="identifier">bits</span><span class="special">,</span> <span class="identifier">it</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"x at minimum = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">",\n f("</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span>
<span class="comment">// x at minimum = 1, f(1) = 5.04853e-018</span>
<span class="special"><<</span> <span class="string">", after "</span> <span class="special"><<</span> <span class="identifier">it</span> <span class="special"><<</span> <span class="string">" iterations. "</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">is_close_to</span><span class="special">(</span><span class="keyword">static_cast</span><span class="special"><</span><span class="identifier">cpp_bin_float_50</span><span class="special">>(</span><span class="string">"1"</span><span class="special">),</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">first</span><span class="special">,</span> <span class="identifier">sqrt</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">cpp_bin_float_50</span><span class="special">>::</span><span class="identifier">epsilon</span><span class="special">()));</span>
<span class="identifier">is_close_to</span><span class="special">(</span><span class="keyword">static_cast</span><span class="special"><</span><span class="identifier">cpp_bin_float_50</span><span class="special">>(</span><span class="string">"0"</span><span class="special">),</span> <span class="identifier">r</span><span class="special">.</span><span class="identifier">second</span><span class="special">,</span> <span class="identifier">sqrt</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">cpp_bin_float_50</span><span class="special">>::</span><span class="identifier">epsilon</span><span class="special">()));</span>
</pre>
<p>
and with our <code class="computeroutput"><span class="identifier">show_minima</span></code> function
</p>
<pre class="programlisting"><span class="identifier">show_minima</span><span class="special"><</span><span class="identifier">cpp_bin_float_50_et_on</span><span class="special">>();</span> <span class="comment">//</span>
</pre>
<pre class="programlisting"><span class="identifier">For</span> <span class="identifier">type</span> <span class="keyword">class</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">number</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">backends</span><span class="special">::</span><span class="identifier">cpp_bin_float</span><span class="special"><</span><span class="number">50</span><span class="special">,</span><span class="number">10</span><span class="special">,</span><span class="keyword">void</span><span class="special">,</span><span class="keyword">int</span><span class="special">,</span><span class="number">0</span><span class="special">,</span><span class="number">0</span><span class="special">>,</span><span class="number">1</span><span class="special">>,</span>
<span class="identifier">epsilon</span> <span class="special">=</span> <span class="number">5.3455294202e-51</span><span class="special">,</span>
<span class="identifier">the</span> <span class="identifier">maximum</span> <span class="identifier">theoretical</span> <span class="identifier">precision</span> <span class="identifier">from</span> <span class="identifier">Brent</span> <span class="identifier">minimization</span> <span class="identifier">is</span> <span class="number">7.311312755e-26</span>
<span class="identifier">Displaying</span> <span class="identifier">to</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">digits10</span> <span class="number">11</span> <span class="identifier">significant</span> <span class="identifier">decimal</span> <span class="identifier">digits</span><span class="special">.</span>
<span class="identifier">x</span> <span class="identifier">at</span> <span class="identifier">minimum</span> <span class="special">=</span> <span class="number">1</span><span class="special">,</span> <span class="identifier">f</span><span class="special">(</span><span class="number">1</span><span class="special">)</span> <span class="special">=</span> <span class="number">5.6273022713e-58</span><span class="special">,</span>
<span class="identifier">met</span> <span class="number">84</span> <span class="identifier">bits</span> <span class="identifier">precision</span><span class="special">,</span> <span class="identifier">after</span> <span class="number">14</span> <span class="identifier">iterations</span><span class="special">.</span>
<span class="identifier">x</span> <span class="special">==</span> <span class="number">1</span> <span class="special">(</span><span class="identifier">compared</span> <span class="identifier">to</span> <span class="identifier">uncertainty</span> <span class="number">7.311312755e-26</span><span class="special">)</span> <span class="identifier">is</span> <span class="keyword">true</span>
<span class="identifier">f</span><span class="special">(</span><span class="identifier">x</span><span class="special">)</span> <span class="special">==</span> <span class="special">(</span><span class="number">0</span> <span class="identifier">compared</span> <span class="identifier">to</span> <span class="identifier">uncertainty</span> <span class="number">7.311312755e-26</span><span class="special">)</span> <span class="identifier">is</span> <span class="keyword">true</span>
<span class="special">-</span><span class="number">4</span> <span class="number">1.3333333333333333333333333333333333333333333333333</span>
<span class="identifier">x</span> <span class="identifier">at</span> <span class="identifier">minimum</span> <span class="special">=</span> <span class="number">0.99999999999999999999999999998813903221565569205253</span><span class="special">,</span>
<span class="identifier">f</span><span class="special">(</span><span class="number">0.99999999999999999999999999998813903221565569205253</span><span class="special">)</span> <span class="special">=</span>
<span class="number">5.6273022712501408640665300316078046703496236636624e-58</span>
<span class="number">14</span> <span class="identifier">iterations</span>
</pre>
<pre class="programlisting"><span class="identifier">For</span> <span class="identifier">type</span> <span class="keyword">class</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">number</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">backends</span><span class="special">::</span><span class="identifier">cpp_bin_float</span><span class="special"><</span><span class="number">50</span><span class="special">,</span> <span class="number">10</span><span class="special">,</span> <span class="keyword">void</span><span class="special">,</span> <span class="keyword">int</span><span class="special">,</span> <span class="number">0</span><span class="special">,</span> <span class="number">0</span><span class="special">>,</span> <span class="number">1</span><span class="special">>,</span>
</pre>
<div class="tip"><table border="0" summary="Tip">
<tr>
<td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../doc/src/images/tip.png"></td>
<th align="left">Tip</th>
</tr>
<tr><td align="left" valign="top"><p>
One can usually rely on template argument deduction to avoid specifying the
verbose multiprecision types, but great care in needed with the <span class="emphasis"><em>type
of the values</em></span> provided to avoid confusing the compiler.
</p></td></tr>
</table></div>
<div class="tip"><table border="0" summary="Tip">
<tr>
<td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../doc/src/images/tip.png"></td>
<th align="left">Tip</th>
</tr>
<tr><td align="left" valign="top"><p>
Using <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><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="identifier">T</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span></code>
or <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><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="identifier">T</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">);</span></code>
during debugging may be wise because it gives some warning if construction
of multiprecision values involves unintended conversion from <code class="computeroutput"><span class="keyword">double</span></code> by showing trailing zero or random
digits after <a href="http://en.cppreference.com/w/cpp/types/numeric_limits/max_digits10" target="_top">max_digits10</a>,
that is 17 for <code class="computeroutput"><span class="keyword">double</span></code>, digit
18... may be just noise.
</p></td></tr>
</table></div>
<p>
The complete example code is at <a href="../../../example/brent_minimise_example.cpp" target="_top">brent_minimise_example.cpp</a>.
</p>
<h5>
<a name="math_toolkit.brent_minima.h6"></a>
<span class="phrase"><a name="math_toolkit.brent_minima.implementation"></a></span><a class="link" href="brent_minima.html#math_toolkit.brent_minima.implementation">Implementation</a>
</h5>
<p>
This is a reasonably faithful implementation of Brent's algorithm.
</p>
<h5>
<a name="math_toolkit.brent_minima.h7"></a>
<span class="phrase"><a name="math_toolkit.brent_minima.references"></a></span><a class="link" href="brent_minima.html#math_toolkit.brent_minima.references">References</a>
</h5>
<div class="orderedlist"><ol class="orderedlist" type="1">
<li class="listitem">
Brent, R.P. 1973, Algorithms for Minimization without Derivatives, (Englewood
Cliffs, NJ: Prentice-Hall), Chapter 5.
</li>
<li class="listitem">
Numerical Recipes in C, The Art of Scientific Computing, Second Edition,
William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P.
Flannery. Cambridge University Press. 1988, 1992.
</li>
<li class="listitem">
An algorithm with guaranteed convergence for finding a zero of a function,
R. P. Brent, The Computer Journal, Vol 44, 1971.
</li>
<li class="listitem">
<a href="http://en.wikipedia.org/wiki/Brent%27s_method" target="_top">Brent's method
in Wikipedia.</a>
</li>
<li class="listitem">
Z. Zhang, An Improvement to the Brent's Method, IJEA, vol. 2, pp. 2 to
26, May 31, 2011. <a href="http://www.cscjournals.org/manuscript/Journals/IJEA/volume2/Issue1/IJEA-7.pdf" target="_top">http://www.cscjournals.org/manuscript/Journals/IJEA/volume2/Issue1/IJEA-7.pdf</a>
</li>
<li class="listitem">
Steven A. Stage, Comments on An Improvement to the Brent's Method (and
comparison of various algorithms) <a href="http://www.cscjournals.org/manuscript/Journals/IJEA/volume4/Issue1/IJEA-33.pdf" target="_top">http://www.cscjournals.org/manuscript/Journals/IJEA/volume4/Issue1/IJEA-33.pdf</a>
Stage concludes that Brent's algorithm is slow to start, but fast to finish
convergence, and has good accuracy.
</li>
</ol></div>
</div>
<div class="copyright-footer">Copyright © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland,
Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno
Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Råde,
Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, 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>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="bad_roots.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../root_finding.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="root_comparison.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
|
