]>
Commit | Line | Data |
---|---|---|
1 | <html> | |
2 | <head> | |
3 | <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII"> | |
4 | <title>Locating Function Minima using Brent's algorithm</title> | |
5 | <link rel="stylesheet" href="../../math.css" type="text/css"> | |
6 | <meta name="generator" content="DocBook XSL Stylesheets V1.77.1"> | |
7 | <link rel="home" href="../../index.html" title="Math Toolkit 2.5.1"> | |
8 | <link rel="up" href="../roots.html" title="Root finding"> | |
9 | <link rel="prev" href="bad_roots.html" title="Examples Where Root Finding Goes Wrong"> | |
10 | <link rel="next" href="root_comparison.html" title="Comparison of Root Finding Algorithms"> | |
11 | </head> | |
12 | <body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> | |
13 | <table cellpadding="2" width="100%"><tr> | |
14 | <td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../boost.png"></td> | |
15 | <td align="center"><a href="../../../../../../index.html">Home</a></td> | |
16 | <td align="center"><a href="../../../../../../libs/libraries.htm">Libraries</a></td> | |
17 | <td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> | |
18 | <td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> | |
19 | <td align="center"><a href="../../../../../../more/index.htm">More</a></td> | |
20 | </tr></table> | |
21 | <hr> | |
22 | <div class="spirit-nav"> | |
23 | <a accesskey="p" href="bad_roots.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../roots.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> | |
24 | </div> | |
25 | <div class="section"> | |
26 | <div class="titlepage"><div><div><h3 class="title"> | |
27 | <a name="math_toolkit.roots.brent_minima"></a><a class="link" href="brent_minima.html" title="Locating Function Minima using Brent's algorithm">Locating Function Minima | |
28 | using Brent's algorithm</a> | |
29 | </h3></div></div></div> | |
30 | <h5> | |
31 | <a name="math_toolkit.roots.brent_minima.h0"></a> | |
32 | <span class="phrase"><a name="math_toolkit.roots.brent_minima.synopsis"></a></span><a class="link" href="brent_minima.html#math_toolkit.roots.brent_minima.synopsis">Synopsis</a> | |
33 | </h5> | |
34 | <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> | |
35 | </pre> | |
36 | <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> | |
37 | <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> | |
38 | ||
39 | <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> | |
40 | <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">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&</span> <span class="identifier">max_iter</span><span class="special">);</span> | |
41 | </pre> | |
42 | <h5> | |
43 | <a name="math_toolkit.roots.brent_minima.h1"></a> | |
44 | <span class="phrase"><a name="math_toolkit.roots.brent_minima.description"></a></span><a class="link" href="brent_minima.html#math_toolkit.roots.brent_minima.description">Description</a> | |
45 | </h5> | |
46 | <p> | |
47 | These two functions locate the minima of the continuous function <span class="emphasis"><em>f</em></span> | |
48 | using <a href="http://en.wikipedia.org/wiki/Brent%27s_method" target="_top">Brent's | |
49 | method</a>: specifically it uses quadratic interpolation to locate the | |
50 | minima, or if that fails, falls back to a <a href="http://en.wikipedia.org/wiki/Golden_section_search" target="_top">golden-section | |
51 | search</a>. | |
52 | </p> | |
53 | <p> | |
54 | <span class="bold"><strong>Parameters</strong></span> | |
55 | </p> | |
56 | <div class="variablelist"> | |
57 | <p class="title"><b></b></p> | |
58 | <dl class="variablelist"> | |
59 | <dt><span class="term">f</span></dt> | |
60 | <dd><p> | |
61 | The function to minimise: a function object (functor) that should be | |
62 | smooth over the range <span class="emphasis"><em>[min, max]</em></span>, with no maxima | |
63 | occurring in that interval. | |
64 | </p></dd> | |
65 | <dt><span class="term">min</span></dt> | |
66 | <dd><p> | |
67 | The lower endpoint of the range in which to search for the minima. | |
68 | </p></dd> | |
69 | <dt><span class="term">max</span></dt> | |
70 | <dd><p> | |
71 | The upper endpoint of the range in which to search for the minima. | |
72 | </p></dd> | |
73 | <dt><span class="term">bits</span></dt> | |
74 | <dd><p> | |
75 | The number of bits precision to which the minima should be found.<br> | |
76 | Note that in principle, the minima can not be located to greater accuracy | |
77 | than the square root of machine epsilon (for 64-bit double, sqrt(1e-16)≅1e-8), | |
78 | therefore the value of <span class="emphasis"><em>bits</em></span> will be ignored if | |
79 | it's greater than half the number of bits in the mantissa of T. | |
80 | </p></dd> | |
81 | <dt><span class="term">max_iter</span></dt> | |
82 | <dd><p> | |
83 | The maximum number of iterations to use in the algorithm, if not provided | |
84 | the algorithm will just keep on going until the minima is found. | |
85 | </p></dd> | |
86 | </dl> | |
87 | </div> | |
88 | <p> | |
89 | <span class="bold"><strong>Returns:</strong></span> | |
90 | </p> | |
91 | <p> | |
92 | A <code class="computeroutput"><span class="identifier">pair</span></code> of type T containing | |
93 | the value of the abscissa at the minima and the value of <span class="emphasis"><em>f(x)</em></span> | |
94 | at the minima. | |
95 | </p> | |
96 | <div class="tip"><table border="0" summary="Tip"> | |
97 | <tr> | |
98 | <td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../../doc/src/images/tip.png"></td> | |
99 | <th align="left">Tip</th> | |
100 | </tr> | |
101 | <tr><td align="left" valign="top"> | |
102 | <p> | |
103 | Defining BOOST_MATH_INSTRUMENT will show some parameters, for example: | |
104 | </p> | |
105 | <pre class="programlisting"><span class="identifier">Type</span> <span class="identifier">T</span> <span class="identifier">is</span> <span class="keyword">double</span> | |
106 | <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> | |
107 | <span class="identifier">tolerance</span> <span class="special">=</span> <span class="number">1.19209289550781e-007</span> | |
108 | <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> | |
109 | <span class="identifier">maximum</span> <span class="identifier">iterations</span> <span class="number">18446744073709551615</span> | |
110 | <span class="number">10</span> <span class="identifier">iterations</span><span class="special">.</span> | |
111 | </pre> | |
112 | </td></tr> | |
113 | </table></div> | |
114 | <h5> | |
115 | <a name="math_toolkit.roots.brent_minima.h2"></a> | |
116 | <span class="phrase"><a name="math_toolkit.roots.brent_minima.example"></a></span><a class="link" href="brent_minima.html#math_toolkit.roots.brent_minima.example">Brent | |
117 | Minimisation Example</a> | |
118 | </h5> | |
119 | <p> | |
120 | As a demonstration, we replicate this <a href="http://en.wikipedia.org/wiki/Brent%27s_method#Example" target="_top">Wikipedia | |
121 | example</a> minimising the function <span class="emphasis"><em>y= (x+3)(x-1)<sup>2</sup></em></span>. | |
122 | </p> | |
123 | <p> | |
124 | It is obvious from the equation and the plot that there is a minimum at exactly | |
125 | one and the value of the function at one is exactly zero. | |
126 | </p> | |
127 | <div class="tip"><table border="0" summary="Tip"> | |
128 | <tr> | |
129 | <td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../../doc/src/images/tip.png"></td> | |
130 | <th align="left">Tip</th> | |
131 | </tr> | |
132 | <tr><td align="left" valign="top"><p> | |
133 | This observation shows that an analytical or <a href="http://en.wikipedia.org/wiki/Closed-form_expression" target="_top">Closed-form | |
134 | expression</a> solution always beats brute-force hands-down for both | |
135 | speed and precision. | |
136 | </p></td></tr> | |
137 | </table></div> | |
138 | <p> | |
139 | <span class="inlinemediaobject"><img src="../../../graphs/brent_test_function_1.svg" align="middle"></span> | |
140 | </p> | |
141 | <p> | |
142 | First an include is needed: | |
143 | </p> | |
144 | <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> | |
145 | </pre> | |
146 | <p> | |
147 | This function is encoded in C++ as function object (functor) using <code class="computeroutput"><span class="keyword">double</span></code> precision thus: | |
148 | </p> | |
149 | <pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">funcdouble</span> | |
150 | <span class="special">{</span> | |
151 | <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> | |
152 | <span class="special">{</span> <span class="comment">//</span> | |
153 | <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> | |
154 | <span class="special">}</span> | |
155 | <span class="special">};</span> | |
156 | </pre> | |
157 | <p> | |
158 | The Brent function is conveniently accessed through a <code class="computeroutput"><span class="keyword">using</span></code> | |
159 | statement (noting sub-namespace <code class="computeroutput"><span class="special">::</span><span class="identifier">tools</span></code>). | |
160 | </p> | |
161 | <p> | |
162 | The search minimum and maximum are chosen as -4 to 4/3 (as in the Wikipedia | |
163 | example). | |
164 | </p> | |
165 | <div class="tip"><table border="0" summary="Tip"> | |
166 | <tr> | |
167 | <td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../../doc/src/images/tip.png"></td> | |
168 | <th align="left">Tip</th> | |
169 | </tr> | |
170 | <tr><td align="left" valign="top"><p> | |
171 | S A Stage (reference 6) reports that the Brent algorithm is <span class="emphasis"><em>slow | |
172 | to start, but fast to converge</em></span>, so choosing a tight min-max | |
173 | range is good. | |
174 | </p></td></tr> | |
175 | </table></div> | |
176 | <p> | |
177 | For simplicity, we set the precision parameter <code class="computeroutput"><span class="identifier">bits</span></code> | |
178 | 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 | |
179 | possible i.e. <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>/2. Nor do we provide a maximum iterations | |
180 | parameter <code class="computeroutput"><span class="identifier">max_iter</span></code>, (perhaps | |
181 | unwidely), so the function will iterate until it finds a minimum. | |
182 | </p> | |
183 | <pre class="programlisting"><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">double</span><span class="special">>::</span><span class="identifier">digits</span><span class="special">;</span> | |
184 | ||
185 | <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> | |
186 | ||
187 | <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> | |
188 | <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> | |
189 | <span class="comment">// x at minimum = 1.00000000112345, f(1.00000000112345) = 5.04852568272458e-018</span> | |
190 | </pre> | |
191 | <p> | |
192 | The resulting <a href="http://en.cppreference.com/w/cpp/utility/pair" target="_top">std::pair</a> | |
193 | contains the minimum close to one and the minimum value close to zero. | |
194 | </p> | |
195 | <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> | |
196 | </pre> | |
197 | <p> | |
198 | The differences from the expected <span class="emphasis"><em>one</em></span> and <span class="emphasis"><em>zero</em></span> | |
199 | are less than the uncertainty (for <code class="computeroutput"><span class="keyword">double</span></code>) | |
200 | 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">digits</span><span class="special">)</span> <span class="special">==</span> <span class="number">53</span></code>. | |
201 | </p> | |
202 | <p> | |
203 | We can use it like this to check that the two values are close-enough to | |
204 | those expected, | |
205 | </p> | |
206 | <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">is_close_to</span><span class="special">;</span> | |
207 | <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> | |
208 | ||
209 | <span class="keyword">double</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="keyword">double</span><span class="special">>::</span><span class="identifier">digits</span><span class="special">);</span> | |
210 | <span class="identifier">is_close_to</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> | |
211 | <span class="identifier">is_small</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> | |
212 | ||
213 | <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">0.00034527</span><span class="special">)</span> <span class="identifier">is</span> <span class="keyword">true</span> | |
214 | <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">0.00034527</span><span class="special">)</span> <span class="identifier">is</span> <span class="keyword">true</span> | |
215 | </pre> | |
216 | <p> | |
217 | It is possible to make this comparison more generally with a templated function, | |
218 | returning <code class="computeroutput"><span class="keyword">true</span></code> when this criterion | |
219 | is met, for example: | |
220 | </p> | |
221 | <pre class="programlisting"><span class="comment">//</span> | |
222 | <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">double</span><span class="special">></span> | |
223 | <span class="keyword">bool</span> <span class="identifier">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> | |
224 | <span class="special">{</span> | |
225 | <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_close_to</span><span class="special">;</span> | |
226 | <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> | |
227 | ||
228 | <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> | |
229 | <span class="special">{</span> | |
230 | <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> | |
231 | <span class="special">}</span> | |
232 | <span class="keyword">else</span> | |
233 | <span class="special">{</span> | |
234 | <span class="keyword">return</span> <span class="identifier">is_close_to</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">got</span><span class="special">,</span> <span class="identifier">tolerance</span><span class="special">);</span> | |
235 | <span class="special">}</span> | |
236 | <span class="special">}</span> | |
237 | </pre> | |
238 | <p> | |
239 | In practical applications, we might want to know how many iterations, and | |
240 | maybe to limit iterations and perhaps to trade some loss of precision for | |
241 | speed, for example: | |
242 | </p> | |
243 | <pre class="programlisting"><span class="keyword">const</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">maxit</span> <span class="special">=</span> <span class="number">20</span><span class="special">;</span> | |
244 | <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span> | |
245 | <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> | |
246 | <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> | |
247 | <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> | |
248 | </pre> | |
249 | <p> | |
250 | limits to a maximum of 20 iterations (a reasonable estimate for this application, | |
251 | even for higher precision shown later). | |
252 | </p> | |
253 | <p> | |
254 | The parameter <code class="computeroutput"><span class="identifier">it</span></code> is updated | |
255 | to return the actual number of iterations (so it may be useful to also keep | |
256 | a record of the limit in <code class="computeroutput"><span class="identifier">maxit</span></code>). | |
257 | </p> | |
258 | <p> | |
259 | It is neat to avoid showing insignificant digits by computing the number | |
260 | of decimal digits to display. | |
261 | </p> | |
262 | <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="identifier">bits</span><span class="special">));</span> <span class="comment">// Number of significant decimal digits.</span> | |
263 | <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 precision with "</span> <span class="special"><<</span> <span class="identifier">prec</span> | |
264 | <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> | |
265 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> | |
266 | <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.</span> | |
267 | ||
268 | <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> | |
269 | <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> | |
270 | </pre> | |
271 | <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> | |
272 | <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> | |
273 | </pre> | |
274 | <p> | |
275 | We can also half the number of precision bits from 52 to 26. | |
276 | </p> | |
277 | <pre class="programlisting"><span class="identifier">bits</span> <span class="special">/=</span> <span class="number">2</span><span class="special">;</span> <span class="comment">// Half digits precision (effective maximum).</span> | |
278 | <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> | |
279 | ||
280 | <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 precision with "</span> <span class="special"><<</span> <span class="identifier">prec</span> | |
281 | <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> | |
282 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> | |
283 | <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.</span> | |
284 | ||
285 | <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span> | |
286 | <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> | |
287 | <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> | |
288 | <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</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> | |
289 | </pre> | |
290 | <p> | |
291 | showing no change in the result and no change in the number of iterations, | |
292 | as expected. | |
293 | </p> | |
294 | <p> | |
295 | It is only if we reduce the precision to a quarter, specifying only 13 precision | |
296 | bits | |
297 | </p> | |
298 | <pre class="programlisting"><span class="identifier">bits</span> <span class="special">/=</span> <span class="number">2</span><span class="special">;</span> <span class="comment">// Quarter precision.</span> | |
299 | <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> | |
300 | ||
301 | <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 precision with "</span> <span class="special"><<</span> <span class="identifier">prec</span> | |
302 | <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> | |
303 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> | |
304 | <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.</span> | |
305 | ||
306 | <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span> | |
307 | <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> | |
308 | <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> | |
309 | <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> | |
310 | </pre> | |
311 | <p> | |
312 | that we reduce the number of iterations from 10 to 7 and the result significantly | |
313 | differing from <span class="emphasis"><em>one</em></span> and <span class="emphasis"><em>zero</em></span>. | |
314 | </p> | |
315 | <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> | |
316 | <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> | |
317 | </pre> | |
318 | <h6> | |
319 | <a name="math_toolkit.roots.brent_minima.h3"></a> | |
320 | <span class="phrase"><a name="math_toolkit.roots.brent_minima.template"></a></span><a class="link" href="brent_minima.html#math_toolkit.roots.brent_minima.template">Templating | |
321 | on floating-point type</a> | |
322 | </h6> | |
323 | <p> | |
324 | If we want to switch the floating-point type, then the functor must be revised. | |
325 | 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> | |
326 | a template member function: | |
327 | </p> | |
328 | <pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">func</span> | |
329 | <span class="special">{</span> | |
330 | <span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> | |
331 | <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> | |
332 | <span class="special">{</span> <span class="comment">//</span> | |
333 | <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">//</span> | |
334 | <span class="special">}</span> | |
335 | <span class="special">};</span> | |
336 | </pre> | |
337 | <p> | |
338 | The <code class="computeroutput"><span class="identifier">brent_find_minima</span></code> function | |
339 | can now be used in template form. | |
340 | </p> | |
341 | <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="keyword">long</span> <span class="keyword">double</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span> | |
342 | <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> | |
343 | <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> | |
344 | <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> | |
345 | <span class="keyword">const</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">maxit</span> <span class="special">=</span> <span class="number">20</span><span class="special">;</span> | |
346 | <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span> | |
347 | ||
348 | <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> | |
349 | <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> | |
350 | <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> | |
351 | </pre> | |
352 | <p> | |
353 | The form shown uses the floating-point type <code class="computeroutput"><span class="keyword">long</span> | |
354 | <span class="keyword">double</span></code> by deduction, but it is also | |
355 | possible to be more explicit, for example: | |
356 | </p> | |
357 | <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> | |
358 | <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> | |
359 | </pre> | |
360 | <p> | |
361 | In order to show the use of multiprecision below, it may be convenient to | |
362 | write a templated function to use this. | |
363 | </p> | |
364 | <pre class="programlisting"><span class="keyword">template</span> <span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> | |
365 | <span class="keyword">void</span> <span class="identifier">show_minima</span><span class="special">()</span> | |
366 | <span class="special">{</span> | |
367 | <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> | |
368 | <span class="keyword">try</span> | |
369 | <span class="special">{</span> <span class="comment">// Always use try'n'catch blocks with Boost.Math to get any error messages.</span> | |
370 | ||
371 | <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> | |
372 | <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="identifier">bits</span><span class="special">));</span> <span class="comment">// Number of significant decimal digits.</span> | |
373 | <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.</span> | |
374 | ||
375 | <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> | |
376 | <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> | |
377 | <span class="comment">// << ", precision of " << bits << " bits"</span> | |
378 | <span class="special"><<</span> <span class="string">",\n the maximum theoretical precision from Brent 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> | |
379 | <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> | |
380 | <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> | |
381 | ||
382 | <span class="keyword">const</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">maxit</span> <span class="special">=</span> <span class="number">20</span><span class="special">;</span> | |
383 | <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span> | |
384 | <span class="comment">// Construct using string, not double, avoids loss of precision.</span> | |
385 | <span class="comment">//T bracket_min = static_cast<T>("-4");</span> | |
386 | <span class="comment">//T bracket_max = static_cast<T>("1.3333333333333333333333333333333333333333333333333");</span> | |
387 | ||
388 | <span class="comment">// Construction from double may cause loss of precision for multiprecision types like cpp_bin_float.</span> | |
389 | <span class="comment">// but brackets values are good enough for using Brent minimization.</span> | |
390 | <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> | |
391 | <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> | |
392 | ||
393 | <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> | |
394 | ||
395 | <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> | |
396 | <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> | |
397 | <span class="special">{</span> | |
398 | <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> | |
399 | <span class="special">}</span> | |
400 | <span class="keyword">else</span> | |
401 | <span class="special">{</span> | |
402 | <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> | |
403 | <span class="special">}</span> | |
404 | <span class="comment">// Check that result is that expected (compared to theoretical uncertainty).</span> | |
405 | <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> | |
406 | <span class="comment">//std::cout << std::boolalpha << "x == 1 (compared to uncertainty " << uncertainty << ") is " << close(static_cast<T>(1), r.first, uncertainty) << std::endl;</span> | |
407 | <span class="comment">//std::cout << std::boolalpha << "f(x) == (0 compared to uncertainty " << uncertainty << ") is " << close(static_cast<T>(0), r.second, uncertainty) << std::endl;</span> | |
408 | <span class="comment">// Problems with this using multiprecision with expression template on?</span> | |
409 | <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> | |
410 | <span class="special">}</span> | |
411 | <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> | |
412 | <span class="special">{</span> <span class="comment">// Always useful to include try & catch blocks because default policies</span> | |
413 | <span class="comment">// are to throw exceptions on arguments that cause errors like underflow, overflow.</span> | |
414 | <span class="comment">// Lacking try & catch blocks, the program will abort without a message below,</span> | |
415 | <span class="comment">// which may give some helpful clues as to the cause of the exception.</span> | |
416 | <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> | |
417 | <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> | |
418 | <span class="special">}</span> | |
419 | <span class="special">}</span> <span class="comment">// void show_minima()</span> | |
420 | </pre> | |
421 | <p> | |
422 | We can use this with all built-in floating-point types, for example | |
423 | </p> | |
424 | <pre class="programlisting"><span class="identifier">show_minima</span><span class="special"><</span><span class="keyword">float</span><span class="special">>();</span> | |
425 | <span class="identifier">show_minima</span><span class="special"><</span><span class="keyword">double</span><span class="special">>();</span> | |
426 | <span class="identifier">show_minima</span><span class="special"><</span><span class="keyword">long</span> <span class="keyword">double</span><span class="special">>();</span> | |
427 | </pre> | |
428 | <p> | |
429 | and, on platforms that provide it, a <a href="http://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format" target="_top">128-bit | |
430 | quad</a> type. (See <a href="../../../../../../libs/multiprecision/doc/html/boost_multiprecision/tut/floats/float128.html" target="_top">float128</a>). | |
431 | </p> | |
432 | <p> | |
433 | For this optional include, the build should define the macro BOOST_HAVE_QUADMATH: | |
434 | </p> | |
435 | <pre class="programlisting"><span class="preprocessor">#ifdef</span> <span class="identifier">BOOST_HAVE_QUADMATH</span> <span class="comment">// Define only if GCC or Intel and have quadmath.lib or .dll library available.</span> | |
436 | <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> | |
437 | <span class="preprocessor">#endif</span> | |
438 | </pre> | |
439 | <p> | |
440 | or | |
441 | </p> | |
442 | <pre class="programlisting"><span class="comment">// #ifndef _MSC_VER</span> | |
443 | <span class="preprocessor">#ifdef</span> <span class="identifier">BOOST_HAVE_QUADMATH</span> <span class="comment">// Define only if GCC or Intel and have quadmath.lib or .dll library available.</span> | |
444 | <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> | |
445 | <span class="preprocessor">#endif</span> | |
446 | </pre> | |
447 | <h6> | |
448 | <a name="math_toolkit.roots.brent_minima.h4"></a> | |
449 | <span class="phrase"><a name="math_toolkit.roots.brent_minima.multiprecision"></a></span><a class="link" href="brent_minima.html#math_toolkit.roots.brent_minima.multiprecision">Multiprecision</a> | |
450 | </h6> | |
451 | <p> | |
452 | If a higher precision than <code class="computeroutput"><span class="keyword">double</span></code> | |
453 | (or <code class="computeroutput"><span class="keyword">long</span> <span class="keyword">double</span></code> | |
454 | if that is more precise) is required, then this is easily achieved using | |
455 | <a href="../../../../../../libs/multiprecision/doc/html/index.html" target="_top">Boost.Multiprecision</a> | |
456 | with some includes from | |
457 | </p> | |
458 | <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> | |
459 | <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> | |
460 | </pre> | |
461 | <p> | |
462 | and some <code class="computeroutput"><span class="identifier">typdef</span></code>s. | |
463 | </p> | |
464 | <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.</span> | |
465 | ||
466 | <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> | |
467 | <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> | |
468 | <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> | |
469 | ||
470 | <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> | |
471 | <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> | |
472 | <span class="identifier">cpp_bin_float_50_et_off</span><span class="special">;</span> | |
473 | ||
474 | <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.</span> | |
475 | ||
476 | <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> | |
477 | <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> | |
478 | <span class="identifier">cpp_dec_float_50_et_on</span><span class="special">;</span> | |
479 | ||
480 | <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> | |
481 | <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> | |
482 | <span class="identifier">cpp_dec_float_50_et_off</span><span class="special">;</span> | |
483 | </pre> | |
484 | <p> | |
485 | Using thus | |
486 | </p> | |
487 | <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> | |
488 | ||
489 | <span class="identifier">cpp_bin_float_50</span> <span class="identifier">fpv</span><span class="special">(</span><span class="string">"-1.2345"</span><span class="special">);</span> | |
490 | <span class="identifier">cpp_bin_float_50</span> <span class="identifier">absv</span><span class="special">;</span> | |
491 | ||
492 | <span class="identifier">absv</span> <span class="special">=</span> <span class="identifier">fpv</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="number">0</span><span class="special">)</span> <span class="special">?</span> <span class="special">-</span><span class="identifier">fpv</span> <span class="special">:</span> <span class="identifier">fpv</span><span class="special">;</span> | |
493 | <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">fpv</span> <span class="special"><<</span> <span class="char">' '</span> <span class="special"><<</span> <span class="identifier">absv</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span> | |
494 | ||
495 | ||
496 | <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> | |
497 | ||
498 | <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> | |
499 | <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> | |
500 | ||
501 | <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">bracket_min</span> <span class="special"><<</span> <span class="string">" "</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> | |
502 | <span class="keyword">const</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">maxit</span> <span class="special">=</span> <span class="number">20</span><span class="special">;</span> | |
503 | <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span> <span class="identifier">it</span> <span class="special">=</span> <span class="identifier">maxit</span><span class="special">;</span> | |
504 | <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> | |
505 | ||
506 | <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> | |
507 | <span class="comment">// x at minimum = 1, f(1) = 5.04853e-018</span> | |
508 | <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> | |
509 | ||
510 | <span class="identifier">close</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="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">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> | |
511 | </pre> | |
512 | <p> | |
513 | and with our show function | |
514 | </p> | |
515 | <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> | |
516 | </pre> | |
517 | <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> | |
518 | <span class="identifier">epsilon</span> <span class="special">=</span> <span class="number">5.3455294202e-51</span><span class="special">,</span> | |
519 | <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> | |
520 | <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> | |
521 | <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> | |
522 | <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> | |
523 | </pre> | |
524 | <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> | |
525 | </pre> | |
526 | <div class="tip"><table border="0" summary="Tip"> | |
527 | <tr> | |
528 | <td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../../doc/src/images/tip.png"></td> | |
529 | <th align="left">Tip</th> | |
530 | </tr> | |
531 | <tr><td align="left" valign="top"><p> | |
532 | One can usually rely on template argument deduction to avoid specifying | |
533 | the verbose multiprecision types, but great care in needed with the <span class="emphasis"><em>type | |
534 | of the values</em></span> provided to avoid confusing the compiler. | |
535 | </p></td></tr> | |
536 | </table></div> | |
537 | <div class="tip"><table border="0" summary="Tip"> | |
538 | <tr> | |
539 | <td rowspan="2" align="center" valign="top" width="25"><img alt="[Tip]" src="../../../../../../doc/src/images/tip.png"></td> | |
540 | <th align="left">Tip</th> | |
541 | </tr> | |
542 | <tr><td align="left" valign="top"><p> | |
543 | 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> | |
544 | 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> | |
545 | during debugging may be wise because it gives some warning if construction | |
546 | of multiprecision values involves unintended conversion from <code class="computeroutput"><span class="keyword">double</span></code> by showing trailing zero or random | |
547 | digits after <a href="http://en.cppreference.com/w/cpp/types/numeric_limits/max_digits10" target="_top">max_digits10</a>, | |
548 | that is 17 for <code class="computeroutput"><span class="keyword">double</span></code>, digit | |
549 | 18... may be just noise. | |
550 | </p></td></tr> | |
551 | </table></div> | |
552 | <p> | |
553 | The complete example code is at <a href="../../../../example/brent_minimise_example.cpp" target="_top">brent_minimise_example.cpp</a>. | |
554 | </p> | |
555 | <h5> | |
556 | <a name="math_toolkit.roots.brent_minima.h5"></a> | |
557 | <span class="phrase"><a name="math_toolkit.roots.brent_minima.implementation"></a></span><a class="link" href="brent_minima.html#math_toolkit.roots.brent_minima.implementation">Implementation</a> | |
558 | </h5> | |
559 | <p> | |
560 | This is a reasonably faithful implementation of Brent's algorithm. | |
561 | </p> | |
562 | <h5> | |
563 | <a name="math_toolkit.roots.brent_minima.h6"></a> | |
564 | <span class="phrase"><a name="math_toolkit.roots.brent_minima.references"></a></span><a class="link" href="brent_minima.html#math_toolkit.roots.brent_minima.references">References</a> | |
565 | </h5> | |
566 | <div class="orderedlist"><ol class="orderedlist" type="1"> | |
567 | <li class="listitem"> | |
568 | Brent, R.P. 1973, Algorithms for Minimization without Derivatives, (Englewood | |
569 | Cliffs, NJ: Prentice-Hall), Chapter 5. | |
570 | </li> | |
571 | <li class="listitem"> | |
572 | Numerical Recipes in C, The Art of Scientific Computing, Second Edition, | |
573 | William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian | |
574 | P. Flannery. Cambridge University Press. 1988, 1992. | |
575 | </li> | |
576 | <li class="listitem"> | |
577 | An algorithm with guaranteed convergence for finding a zero of a function, | |
578 | R. P. Brent, The Computer Journal, Vol 44, 1971. | |
579 | </li> | |
580 | <li class="listitem"> | |
581 | <a href="http://en.wikipedia.org/wiki/Brent%27s_method" target="_top">Brent's method | |
582 | in Wikipedia.</a> | |
583 | </li> | |
584 | <li class="listitem"> | |
585 | Z. Zhang, An Improvement to the Brent's Method, IJEA, vol. 2, pp. 2 to | |
586 | 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> | |
587 | </li> | |
588 | <li class="listitem"> | |
589 | Steven A. Stage, Comments on An Improvement to the Brent's Method (and | |
590 | 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> | |
591 | Stage concludes that Brent's algorithm is slow to start, but fast to | |
592 | finish convergence, and has good accuracy. | |
593 | </li> | |
594 | </ol></div> | |
595 | </div> | |
596 | <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> | |
597 | <td align="left"></td> | |
598 | <td align="right"><div class="copyright-footer">Copyright © 2006-2010, 2012-2014 Nikhar Agrawal, | |
599 | Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert | |
600 | Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Johan Råde, Gautam Sewani, | |
601 | Benjamin Sobotta, Thijs van den Berg, Daryle Walker and Xiaogang Zhang<p> | |
602 | Distributed under the Boost Software License, Version 1.0. (See accompanying | |
603 | 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>) | |
604 | </p> | |
605 | </div></td> | |
606 | </tr></table> | |
607 | <hr> | |
608 | <div class="spirit-nav"> | |
609 | <a accesskey="p" href="bad_roots.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../roots.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> | |
610 | </div> | |
611 | </body> | |
612 | </html> |