1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Root Finding With Derivatives: Newton-Raphson, Halley & Schröder</title> 5<link rel="stylesheet" href="../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../root_finding.html" title="Chapter 10. Root Finding & Minimization Algorithms"> 9<link rel="prev" href="roots_noderiv/implementation.html" title="Implementation"> 10<link rel="next" href="root_finding_examples.html" title="Examples of Root-Finding (with and without derivatives)"> 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="roots_noderiv/implementation.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../root_finding.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="root_finding_examples.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h2 class="title" style="clear: both"> 27<a name="math_toolkit.roots_deriv"></a><a class="link" href="roots_deriv.html" title="Root Finding With Derivatives: Newton-Raphson, Halley & Schröder">Root Finding With Derivatives: 28 Newton-Raphson, Halley & Schröder</a> 29</h2></div></div></div> 30<h5> 31<a name="math_toolkit.roots_deriv.h0"></a> 32 <span class="phrase"><a name="math_toolkit.roots_deriv.synopsis"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.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">roots</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 35</pre> 36<pre class="programlisting"><span class="keyword">namespace</span> <span class="identifier">boost</span> <span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span> <span class="special">{</span> 37<span class="keyword">namespace</span> <span class="identifier">tools</span> <span class="special">{</span> <span class="comment">// Note namespace boost::math::tools.</span> 38<span class="comment">// Newton-Raphson</span> 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">T</span> <span class="identifier">newton_raphson_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">);</span> 41 42<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> 43<span class="identifier">T</span> <span class="identifier">newton_raphson_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&</span> <span class="identifier">max_iter</span><span class="special">);</span> 44 45<span class="comment">// Halley</span> 46<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> 47<span class="identifier">T</span> <span class="identifier">halley_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">);</span> 48 49<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> 50<span class="identifier">T</span> <span class="identifier">halley_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&</span> <span class="identifier">max_iter</span><span class="special">);</span> 51 52<span class="comment">// Schr'''&#xf6;'''der</span> 53<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> 54<span class="identifier">T</span> <span class="identifier">schroder_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">);</span> 55 56<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> 57<span class="identifier">T</span> <span class="identifier">schroder_iterate</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">guess</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">min</span><span class="special">,</span> <span class="identifier">T</span> <span class="identifier">max</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">digits</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">&</span> <span class="identifier">max_iter</span><span class="special">);</span> 58 59<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">Complex</span><span class="special">></span> 60<span class="identifier">Complex</span> <span class="identifier">complex_newton</span><span class="special">(</span><span class="identifier">F</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Complex</span> <span class="identifier">guess</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">max_iterations</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">typename</span> <span class="identifier">Complex</span><span class="special">::</span><span class="identifier">value_type</span><span class="special">>::</span><span class="identifier">digits</span><span class="special">);</span> 61 62<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">T</span><span class="special">></span> 63<span class="keyword">auto</span> <span class="identifier">quadratic_roots</span><span class="special">(</span><span class="identifier">T</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">a</span><span class="special">,</span> <span class="identifier">T</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">b</span><span class="special">,</span> <span class="identifier">T</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">c</span><span class="special">);</span> 64 65<span class="special">}}}</span> <span class="comment">// namespaces boost::math::tools.</span> 66</pre> 67<h5> 68<a name="math_toolkit.roots_deriv.h1"></a> 69 <span class="phrase"><a name="math_toolkit.roots_deriv.description"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.description">Description</a> 70 </h5> 71<p> 72 These functions all perform iterative root-finding <span class="bold"><strong>using 73 derivatives</strong></span>: 74 </p> 75<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 76<li class="listitem"> 77 <code class="computeroutput"><span class="identifier">newton_raphson_iterate</span></code> 78 performs second-order <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.newton">Newton-Raphson 79 iteration</a>. 80 </li> 81<li class="listitem"> 82 <code class="computeroutput"><span class="identifier">halley_iterate</span></code> and <code class="computeroutput"><span class="identifier">schroder_iterate</span></code> perform third-order 83 <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.halley">Halley</a> and <a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.schroder">Schröder</a> iteration. 84 </li> 85<li class="listitem"> 86 <code class="computeroutput"><span class="identifier">complex_newton</span></code> performs 87 Newton's method on complex-analytic functions. 88 </li> 89<li class="listitem"> 90 <code class="computeroutput"><span class="identifier">solve_quadratic</span></code> solves 91 quadratic equations using various tricks to keep catastrophic cancellation 92 from occurring in computation of the discriminant. 93 </li> 94</ul></div> 95<div class="variablelist"> 96<p class="title"><b>Parameters of the real-valued root finding functions</b></p> 97<dl class="variablelist"> 98<dt><span class="term">F f</span></dt> 99<dd> 100<p> 101 Type F must be a callable function object (or C++ lambda) that accepts 102 one parameter and returns a <a class="link" href="internals/tuples.html" title="Tuples">std::pair, 103 std::tuple, boost::tuple or boost::fusion::tuple</a>: 104 </p> 105<p> 106 For second-order iterative method (<a href="http://en.wikipedia.org/wiki/Newton_Raphson" target="_top">Newton 107 Raphson</a>) the <code class="computeroutput"><span class="identifier">tuple</span></code> 108 should have <span class="bold"><strong>two</strong></span> elements containing 109 the evaluation of the function and its first derivative. 110 </p> 111<p> 112 For the third-order methods (<a href="http://en.wikipedia.org/wiki/Halley%27s_method" target="_top">Halley</a> 113 and Schröder) the <code class="computeroutput"><span class="identifier">tuple</span></code> 114 should have <span class="bold"><strong>three</strong></span> elements containing 115 the evaluation of the function and its first and second derivatives. 116 </p> 117</dd> 118<dt><span class="term">T guess</span></dt> 119<dd><p> 120 The initial starting value. A good guess is crucial to quick convergence! 121 </p></dd> 122<dt><span class="term">T min</span></dt> 123<dd><p> 124 The minimum possible value for the result, this is used as an initial 125 lower bracket. 126 </p></dd> 127<dt><span class="term">T max</span></dt> 128<dd><p> 129 The maximum possible value for the result, this is used as an initial 130 upper bracket. 131 </p></dd> 132<dt><span class="term">int digits</span></dt> 133<dd><p> 134 The desired number of binary digits precision. 135 </p></dd> 136<dt><span class="term">uintmax_t& max_iter</span></dt> 137<dd><p> 138 An optional maximum number of iterations to perform. On exit, this is 139 updated to the actual number of iterations performed. 140 </p></dd> 141</dl> 142</div> 143<p> 144 When using these functions you should note that: 145 </p> 146<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 147<li class="listitem"> 148 Default <code class="computeroutput"><span class="identifier">max_iter</span> <span class="special">=</span> 149 <span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">uintmax_t</span><span class="special">>::</span><span class="identifier">max</span><span class="special">)()</span></code> is effectively 'iterate for ever'. 150 </li> 151<li class="listitem"> 152 They may be very sensitive to the initial guess, typically they converge 153 very rapidly if the initial guess has two or three decimal digits correct. 154 However convergence can be no better than <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a>, 155 or in some rare cases, even worse than <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a> 156 if the initial guess is a long way from the correct value and the derivatives 157 are close to zero. 158 </li> 159<li class="listitem"> 160 These functions include special cases to handle zero first (and second 161 where appropriate) derivatives, and fall back to <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a> 162 in this case. However, it is helpful if functor F is defined to return 163 an arbitrarily small value <span class="emphasis"><em>of the correct sign</em></span> rather 164 than zero. 165 </li> 166<li class="listitem"> 167 The functions will raise an <a class="link" href="error_handling.html#math_toolkit.error_handling.evaluation_error">evaluation_error</a> 168 if arguments <code class="computeroutput"><span class="identifier">min</span></code> and <code class="computeroutput"><span class="identifier">max</span></code> are the wrong way around or if they 169 converge to a local minima. 170 </li> 171<li class="listitem"> 172 If the derivative at the current best guess for the result is infinite 173 (or very close to being infinite) then these functions may terminate prematurely. 174 A large first derivative leads to a very small next step, triggering the 175 termination condition. Derivative based iteration may not be appropriate 176 in such cases. 177 </li> 178<li class="listitem"> 179 If the function is 'Really Well Behaved' (is monotonic and has only one 180 root) the bracket bounds <span class="emphasis"><em>min</em></span> and <span class="emphasis"><em>max</em></span> 181 may as well be set to the widest limits like zero and <code class="computeroutput"><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">max</span><span class="special">()</span></code>. 182 </li> 183<li class="listitem"> 184 But if the function more complex and may have more than one root or a pole, 185 the choice of bounds is protection against jumping out to seek the 'wrong' 186 root. 187 </li> 188<li class="listitem"> 189 These functions fall back to <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisect</a> 190 if the next computed step would take the next value out of bounds. The 191 bounds are updated after each step to ensure this leads to convergence. 192 However, a good initial guess backed up by asymptotically-tight bounds 193 will improve performance no end - rather than relying on <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisection</a>. 194 </li> 195<li class="listitem"> 196 The value of <span class="emphasis"><em>digits</em></span> is crucial to good performance 197 of these functions, if it is set too high then at best you will get one 198 extra (unnecessary) iteration, and at worst the last few steps will proceed 199 by <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisection</a>. 200 Remember that the returned value can never be more accurate than <span class="emphasis"><em>f(x)</em></span> 201 can be evaluated, and that if <span class="emphasis"><em>f(x)</em></span> suffers from cancellation 202 errors as it tends to zero then the computed steps will be effectively 203 random. The value of <span class="emphasis"><em>digits</em></span> should be set so that 204 iteration terminates before this point: remember that for second and third 205 order methods the number of correct digits in the result is increasing 206 quite substantially with each iteration, <span class="emphasis"><em>digits</em></span> should 207 be set by experiment so that the final iteration just takes the next value 208 into the zone where <span class="emphasis"><em>f(x)</em></span> becomes inaccurate. A good 209 starting point for <span class="emphasis"><em>digits</em></span> would be 0.6*D for Newton 210 and 0.4*D for Halley or Shröder iteration, where D is <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">T</span><span class="special">>::</span><span class="identifier">digits</span></code>. 211 </li> 212<li class="listitem"> 213 If you need some diagnostic output to see what is going on, you can <code class="computeroutput"><span class="preprocessor">#define</span> <span class="identifier">BOOST_MATH_INSTRUMENT</span></code> 214 before the <code class="computeroutput"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">tools</span><span class="special">/</span><span class="identifier">roots</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span></code>, and also ensure that display of all 215 the significant digits with <code class="computeroutput"> <span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">)</span></code>: or even possibly significant digits with 216 <code class="computeroutput"> <span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">)</span></code>: 217 but be warned, this may produce copious output! 218 </li> 219<li class="listitem"> 220 Finally: you may well be able to do better than these functions by hand-coding 221 the heuristics used so that they are tailored to a specific function. You 222 may also be able to compute the ratio of derivatives used by these methods 223 more efficiently than computing the derivatives themselves. As ever, algebraic 224 simplification can be a big win. 225 </li> 226</ul></div> 227<h5> 228<a name="math_toolkit.roots_deriv.h2"></a> 229 <span class="phrase"><a name="math_toolkit.roots_deriv.newton"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.newton">Newton 230 Raphson Method</a> 231 </h5> 232<p> 233 Given an initial guess <span class="emphasis"><em>x0</em></span> the subsequent values are computed 234 using: 235 </p> 236<div class="blockquote"><blockquote class="blockquote"><p> 237 <span class="inlinemediaobject"><img src="../../equations/roots1.svg"></span> 238 239 </p></blockquote></div> 240<p> 241 Out-of-bounds steps revert to <a class="link" href="roots_noderiv/bisect.html" title="Bisection">bisection</a> 242 of the current bounds. 243 </p> 244<p> 245 Under ideal conditions, the number of correct digits doubles with each iteration. 246 </p> 247<h5> 248<a name="math_toolkit.roots_deriv.h3"></a> 249 <span class="phrase"><a name="math_toolkit.roots_deriv.halley"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.halley">Halley's 250 Method</a> 251 </h5> 252<p> 253 Given an initial guess <span class="emphasis"><em>x0</em></span> the subsequent values are computed 254 using: 255 </p> 256<div class="blockquote"><blockquote class="blockquote"><p> 257 <span class="inlinemediaobject"><img src="../../equations/roots2.svg"></span> 258 259 </p></blockquote></div> 260<p> 261 Over-compensation by the second derivative (one which would proceed in the 262 wrong direction) causes the method to revert to a Newton-Raphson step. 263 </p> 264<p> 265 Out of bounds steps revert to bisection of the current bounds. 266 </p> 267<p> 268 Under ideal conditions, the number of correct digits trebles with each iteration. 269 </p> 270<h5> 271<a name="math_toolkit.roots_deriv.h4"></a> 272 <span class="phrase"><a name="math_toolkit.roots_deriv.schroder"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.schroder">Schröder's 273 Method</a> 274 </h5> 275<p> 276 Given an initial guess x0 the subsequent values are computed using: 277 </p> 278<div class="blockquote"><blockquote class="blockquote"><p> 279 <span class="inlinemediaobject"><img src="../../equations/roots3.svg"></span> 280 281 </p></blockquote></div> 282<p> 283 Over-compensation by the second derivative (one which would proceed in the 284 wrong direction) causes the method to revert to a Newton-Raphson step. Likewise 285 a Newton step is used whenever that Newton step would change the next value 286 by more than 10%. 287 </p> 288<p> 289 Out of bounds steps revert to <a href="https://en.wikipedia.org/wiki/Bisection" target="_top">bisection</a> 290 of the current bounds. 291 </p> 292<p> 293 Under ideal conditions, the number of correct digits trebles with each iteration. 294 </p> 295<p> 296 This is Schröder's general result (equation 18 from <a href="http://drum.lib.umd.edu/handle/1903/577" target="_top">Stewart, 297 G. W. "On Infinitely Many Algorithms for Solving Equations." English 298 translation of Schröder's original paper. College Park, MD: University of Maryland, 299 Institute for Advanced Computer Studies, Department of Computer Science, 1993</a>.) 300 </p> 301<p> 302 This method guarantees at least quadratic convergence (the same as Newton's 303 method), and is known to work well in the presence of multiple roots: something 304 that neither Newton nor Halley can do. 305 </p> 306<p> 307 The complex Newton method works slightly differently than the rest of the methods: 308 Since there is no way to bracket roots in the complex plane, the <code class="computeroutput"><span class="identifier">min</span></code> and <code class="computeroutput"><span class="identifier">max</span></code> 309 arguments are not accepted. Failure to reach a root is communicated by returning 310 <code class="computeroutput"><span class="identifier">nan</span></code>s. Remember that if a function 311 has many roots, then which root the complex Newton's method converges to is 312 essentially impossible to predict a priori; see the Newton's fractal for more 313 information. 314 </p> 315<p> 316 Finally, the derivative of <span class="emphasis"><em>f</em></span> must be continuous at the 317 root or else non-roots can be found; see <a href="https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root" target="_top">here</a> 318 for an example. 319 </p> 320<p> 321 An example usage of <code class="computeroutput"><span class="identifier">complex_newton</span></code> 322 is given in <code class="computeroutput"><span class="identifier">examples</span><span class="special">/</span><span class="identifier">daubechies_coefficients</span><span class="special">.</span><span class="identifier">cpp</span></code>. 323 </p> 324<h5> 325<a name="math_toolkit.roots_deriv.h5"></a> 326 <span class="phrase"><a name="math_toolkit.roots_deriv.quadratics"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.quadratics">Quadratics</a> 327 </h5> 328<p> 329 To solve a quadratic <span class="emphasis"><em>ax</em></span><sup>2</sup> + <span class="emphasis"><em>bx</em></span> + <span class="emphasis"><em>c</em></span> 330 = 0, we may use 331 </p> 332<pre class="programlisting"><span class="keyword">auto</span> <span class="special">[</span><span class="identifier">x0</span><span class="special">,</span> <span class="identifier">x1</span><span class="special">]</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">quadratic_roots</span><span class="special">(</span><span class="identifier">a</span><span class="special">,</span> <span class="identifier">b</span><span class="special">,</span> <span class="identifier">c</span><span class="special">);</span> 333</pre> 334<p> 335 If the roots are real, they are arranged so that <code class="computeroutput"><span class="identifier">x0</span></code> 336 ≤ <code class="computeroutput"><span class="identifier">x1</span></code>. If the roots are 337 complex and the inputs are real, <code class="computeroutput"><span class="identifier">x0</span></code> 338 and <code class="computeroutput"><span class="identifier">x1</span></code> are both <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="identifier">Real</span><span class="special">>::</span><span class="identifier">quiet_NaN</span><span class="special">()</span></code>. In this case we must cast <code class="computeroutput"><span class="identifier">a</span></code>, <code class="computeroutput"><span class="identifier">b</span></code> 339 and <code class="computeroutput"><span class="identifier">c</span></code> to a complex type to 340 extract the complex roots. If <code class="computeroutput"><span class="identifier">a</span></code>, 341 <code class="computeroutput"><span class="identifier">b</span></code> and <code class="computeroutput"><span class="identifier">c</span></code> 342 are integral, then the roots are of type double. The routine is much faster 343 if the fused-multiply-add instruction is available on your architecture. If 344 the fma is not available, the function resorts to slow emulation. Finally, 345 speed is improved if you compile for your particular architecture. For instance, 346 if you compile without any architecture flags, then the <code class="computeroutput"><span class="identifier">std</span><span class="special">::</span><span class="identifier">fma</span></code> call 347 compiles down to <code class="computeroutput"><span class="identifier">call</span> <span class="identifier">_fma</span></code>, 348 which dynamically chooses to emulate or execute the <code class="computeroutput"><span class="identifier">vfmadd132sd</span></code> 349 instruction based on the capabilities of the architecture. If instead, you 350 compile with (say) <code class="computeroutput"><span class="special">-</span><span class="identifier">march</span><span class="special">=</span><span class="identifier">native</span></code> then 351 no dynamic choice is made: The <code class="computeroutput"><span class="identifier">vfmadd132sd</span></code> 352 instruction is always executed if available and emulation is used if not. 353 </p> 354<h5> 355<a name="math_toolkit.roots_deriv.h6"></a> 356 <span class="phrase"><a name="math_toolkit.roots_deriv.examples"></a></span><a class="link" href="roots_deriv.html#math_toolkit.roots_deriv.examples">Examples</a> 357 </h5> 358<p> 359 See <a class="link" href="root_finding_examples.html" title="Examples of Root-Finding (with and without derivatives)">root-finding examples</a>. 360 </p> 361</div> 362<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 363<td align="left"></td> 364<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 365 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 366 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 367 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 368 Daryle Walker and Xiaogang Zhang<p> 369 Distributed under the Boost Software License, Version 1.0. (See accompanying 370 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>) 371 </p> 372</div></td> 373</tr></table> 374<hr> 375<div class="spirit-nav"> 376<a accesskey="p" href="roots_noderiv/implementation.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../root_finding.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="root_finding_examples.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 377</div> 378</body> 379</html> 380
