projects / ceph.git / blame
| Commit | Line | Data |
|---|---|---|
| 7c673cae FG |
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> |