]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright John Maddock 2006. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #ifndef BOOST_MATH_TOOLS_SOLVE_ROOT_HPP | |
7 | #define BOOST_MATH_TOOLS_SOLVE_ROOT_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
13 | #include <boost/math/tools/precision.hpp> | |
14 | #include <boost/math/policies/error_handling.hpp> | |
15 | #include <boost/math/tools/config.hpp> | |
16 | #include <boost/math/special_functions/sign.hpp> | |
17 | #include <boost/cstdint.hpp> | |
18 | #include <limits> | |
19 | ||
20 | #ifdef BOOST_MATH_LOG_ROOT_ITERATIONS | |
21 | # define BOOST_MATH_LOGGER_INCLUDE <boost/math/tools/iteration_logger.hpp> | |
22 | # include BOOST_MATH_LOGGER_INCLUDE | |
23 | # undef BOOST_MATH_LOGGER_INCLUDE | |
24 | #else | |
25 | # define BOOST_MATH_LOG_COUNT(count) | |
26 | #endif | |
27 | ||
28 | namespace boost{ namespace math{ namespace tools{ | |
29 | ||
30 | template <class T> | |
31 | class eps_tolerance | |
32 | { | |
33 | public: | |
34 | eps_tolerance() | |
35 | { | |
36 | eps = 4 * tools::epsilon<T>(); | |
37 | } | |
38 | eps_tolerance(unsigned bits) | |
39 | { | |
40 | BOOST_MATH_STD_USING | |
41 | eps = (std::max)(T(ldexp(1.0F, 1-bits)), T(4 * tools::epsilon<T>())); | |
42 | } | |
43 | bool operator()(const T& a, const T& b) | |
44 | { | |
45 | BOOST_MATH_STD_USING | |
46 | return fabs(a - b) <= (eps * (std::min)(fabs(a), fabs(b))); | |
47 | } | |
48 | private: | |
49 | T eps; | |
50 | }; | |
51 | ||
52 | struct equal_floor | |
53 | { | |
54 | equal_floor(){} | |
55 | template <class T> | |
56 | bool operator()(const T& a, const T& b) | |
57 | { | |
58 | BOOST_MATH_STD_USING | |
59 | return floor(a) == floor(b); | |
60 | } | |
61 | }; | |
62 | ||
63 | struct equal_ceil | |
64 | { | |
65 | equal_ceil(){} | |
66 | template <class T> | |
67 | bool operator()(const T& a, const T& b) | |
68 | { | |
69 | BOOST_MATH_STD_USING | |
70 | return ceil(a) == ceil(b); | |
71 | } | |
72 | }; | |
73 | ||
74 | struct equal_nearest_integer | |
75 | { | |
76 | equal_nearest_integer(){} | |
77 | template <class T> | |
78 | bool operator()(const T& a, const T& b) | |
79 | { | |
80 | BOOST_MATH_STD_USING | |
81 | return floor(a + 0.5f) == floor(b + 0.5f); | |
82 | } | |
83 | }; | |
84 | ||
85 | namespace detail{ | |
86 | ||
87 | template <class F, class T> | |
88 | void bracket(F f, T& a, T& b, T c, T& fa, T& fb, T& d, T& fd) | |
89 | { | |
90 | // | |
91 | // Given a point c inside the existing enclosing interval | |
92 | // [a, b] sets a = c if f(c) == 0, otherwise finds the new | |
93 | // enclosing interval: either [a, c] or [c, b] and sets | |
94 | // d and fd to the point that has just been removed from | |
95 | // the interval. In other words d is the third best guess | |
96 | // to the root. | |
97 | // | |
98 | BOOST_MATH_STD_USING // For ADL of std math functions | |
99 | T tol = tools::epsilon<T>() * 2; | |
100 | // | |
101 | // If the interval [a,b] is very small, or if c is too close | |
102 | // to one end of the interval then we need to adjust the | |
103 | // location of c accordingly: | |
104 | // | |
105 | if((b - a) < 2 * tol * a) | |
106 | { | |
107 | c = a + (b - a) / 2; | |
108 | } | |
109 | else if(c <= a + fabs(a) * tol) | |
110 | { | |
111 | c = a + fabs(a) * tol; | |
112 | } | |
113 | else if(c >= b - fabs(b) * tol) | |
114 | { | |
115 | c = b - fabs(b) * tol; | |
116 | } | |
117 | // | |
118 | // OK, lets invoke f(c): | |
119 | // | |
120 | T fc = f(c); | |
121 | // | |
122 | // if we have a zero then we have an exact solution to the root: | |
123 | // | |
124 | if(fc == 0) | |
125 | { | |
126 | a = c; | |
127 | fa = 0; | |
128 | d = 0; | |
129 | fd = 0; | |
130 | return; | |
131 | } | |
132 | // | |
133 | // Non-zero fc, update the interval: | |
134 | // | |
135 | if(boost::math::sign(fa) * boost::math::sign(fc) < 0) | |
136 | { | |
137 | d = b; | |
138 | fd = fb; | |
139 | b = c; | |
140 | fb = fc; | |
141 | } | |
142 | else | |
143 | { | |
144 | d = a; | |
145 | fd = fa; | |
146 | a = c; | |
147 | fa= fc; | |
148 | } | |
149 | } | |
150 | ||
151 | template <class T> | |
152 | inline T safe_div(T num, T denom, T r) | |
153 | { | |
154 | // | |
155 | // return num / denom without overflow, | |
156 | // return r if overflow would occur. | |
157 | // | |
158 | BOOST_MATH_STD_USING // For ADL of std math functions | |
159 | ||
160 | if(fabs(denom) < 1) | |
161 | { | |
162 | if(fabs(denom * tools::max_value<T>()) <= fabs(num)) | |
163 | return r; | |
164 | } | |
165 | return num / denom; | |
166 | } | |
167 | ||
168 | template <class T> | |
169 | inline T secant_interpolate(const T& a, const T& b, const T& fa, const T& fb) | |
170 | { | |
171 | // | |
172 | // Performs standard secant interpolation of [a,b] given | |
173 | // function evaluations f(a) and f(b). Performs a bisection | |
174 | // if secant interpolation would leave us very close to either | |
175 | // a or b. Rationale: we only call this function when at least | |
176 | // one other form of interpolation has already failed, so we know | |
177 | // that the function is unlikely to be smooth with a root very | |
178 | // close to a or b. | |
179 | // | |
180 | BOOST_MATH_STD_USING // For ADL of std math functions | |
181 | ||
182 | T tol = tools::epsilon<T>() * 5; | |
183 | T c = a - (fa / (fb - fa)) * (b - a); | |
184 | if((c <= a + fabs(a) * tol) || (c >= b - fabs(b) * tol)) | |
185 | return (a + b) / 2; | |
186 | return c; | |
187 | } | |
188 | ||
189 | template <class T> | |
190 | T quadratic_interpolate(const T& a, const T& b, T const& d, | |
191 | const T& fa, const T& fb, T const& fd, | |
192 | unsigned count) | |
193 | { | |
194 | // | |
195 | // Performs quadratic interpolation to determine the next point, | |
196 | // takes count Newton steps to find the location of the | |
197 | // quadratic polynomial. | |
198 | // | |
199 | // Point d must lie outside of the interval [a,b], it is the third | |
200 | // best approximation to the root, after a and b. | |
201 | // | |
202 | // Note: this does not guarantee to find a root | |
203 | // inside [a, b], so we fall back to a secant step should | |
204 | // the result be out of range. | |
205 | // | |
206 | // Start by obtaining the coefficients of the quadratic polynomial: | |
207 | // | |
208 | T B = safe_div(T(fb - fa), T(b - a), tools::max_value<T>()); | |
209 | T A = safe_div(T(fd - fb), T(d - b), tools::max_value<T>()); | |
210 | A = safe_div(T(A - B), T(d - a), T(0)); | |
211 | ||
212 | if(A == 0) | |
213 | { | |
214 | // failure to determine coefficients, try a secant step: | |
215 | return secant_interpolate(a, b, fa, fb); | |
216 | } | |
217 | // | |
218 | // Determine the starting point of the Newton steps: | |
219 | // | |
220 | T c; | |
221 | if(boost::math::sign(A) * boost::math::sign(fa) > 0) | |
222 | { | |
223 | c = a; | |
224 | } | |
225 | else | |
226 | { | |
227 | c = b; | |
228 | } | |
229 | // | |
230 | // Take the Newton steps: | |
231 | // | |
232 | for(unsigned i = 1; i <= count; ++i) | |
233 | { | |
234 | //c -= safe_div(B * c, (B + A * (2 * c - a - b)), 1 + c - a); | |
235 | c -= safe_div(T(fa+(B+A*(c-b))*(c-a)), T(B + A * (2 * c - a - b)), T(1 + c - a)); | |
236 | } | |
237 | if((c <= a) || (c >= b)) | |
238 | { | |
239 | // Oops, failure, try a secant step: | |
240 | c = secant_interpolate(a, b, fa, fb); | |
241 | } | |
242 | return c; | |
243 | } | |
244 | ||
245 | template <class T> | |
246 | T cubic_interpolate(const T& a, const T& b, const T& d, | |
247 | const T& e, const T& fa, const T& fb, | |
248 | const T& fd, const T& fe) | |
249 | { | |
250 | // | |
251 | // Uses inverse cubic interpolation of f(x) at points | |
252 | // [a,b,d,e] to obtain an approximate root of f(x). | |
253 | // Points d and e lie outside the interval [a,b] | |
254 | // and are the third and forth best approximations | |
255 | // to the root that we have found so far. | |
256 | // | |
257 | // Note: this does not guarantee to find a root | |
258 | // inside [a, b], so we fall back to quadratic | |
259 | // interpolation in case of an erroneous result. | |
260 | // | |
261 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b | |
262 | << " d = " << d << " e = " << e << " fa = " << fa << " fb = " << fb | |
263 | << " fd = " << fd << " fe = " << fe); | |
264 | T q11 = (d - e) * fd / (fe - fd); | |
265 | T q21 = (b - d) * fb / (fd - fb); | |
266 | T q31 = (a - b) * fa / (fb - fa); | |
267 | T d21 = (b - d) * fd / (fd - fb); | |
268 | T d31 = (a - b) * fb / (fb - fa); | |
269 | BOOST_MATH_INSTRUMENT_CODE( | |
270 | "q11 = " << q11 << " q21 = " << q21 << " q31 = " << q31 | |
271 | << " d21 = " << d21 << " d31 = " << d31); | |
272 | T q22 = (d21 - q11) * fb / (fe - fb); | |
273 | T q32 = (d31 - q21) * fa / (fd - fa); | |
274 | T d32 = (d31 - q21) * fd / (fd - fa); | |
275 | T q33 = (d32 - q22) * fa / (fe - fa); | |
276 | T c = q31 + q32 + q33 + a; | |
277 | BOOST_MATH_INSTRUMENT_CODE( | |
278 | "q22 = " << q22 << " q32 = " << q32 << " d32 = " << d32 | |
279 | << " q33 = " << q33 << " c = " << c); | |
280 | ||
281 | if((c <= a) || (c >= b)) | |
282 | { | |
283 | // Out of bounds step, fall back to quadratic interpolation: | |
284 | c = quadratic_interpolate(a, b, d, fa, fb, fd, 3); | |
285 | BOOST_MATH_INSTRUMENT_CODE( | |
286 | "Out of bounds interpolation, falling back to quadratic interpolation. c = " << c); | |
287 | } | |
288 | ||
289 | return c; | |
290 | } | |
291 | ||
292 | } // namespace detail | |
293 | ||
294 | template <class F, class T, class Tol, class Policy> | |
295 | std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, boost::uintmax_t& max_iter, const Policy& pol) | |
296 | { | |
297 | // | |
298 | // Main entry point and logic for Toms Algorithm 748 | |
299 | // root finder. | |
300 | // | |
301 | BOOST_MATH_STD_USING // For ADL of std math functions | |
302 | ||
303 | static const char* function = "boost::math::tools::toms748_solve<%1%>"; | |
304 | ||
305 | boost::uintmax_t count = max_iter; | |
306 | T a, b, fa, fb, c, u, fu, a0, b0, d, fd, e, fe; | |
307 | static const T mu = 0.5f; | |
308 | ||
309 | // initialise a, b and fa, fb: | |
310 | a = ax; | |
311 | b = bx; | |
312 | if(a >= b) | |
313 | return boost::math::detail::pair_from_single(policies::raise_domain_error( | |
314 | function, | |
315 | "Parameters a and b out of order: a=%1%", a, pol)); | |
316 | fa = fax; | |
317 | fb = fbx; | |
318 | ||
319 | if(tol(a, b) || (fa == 0) || (fb == 0)) | |
320 | { | |
321 | max_iter = 0; | |
322 | if(fa == 0) | |
323 | b = a; | |
324 | else if(fb == 0) | |
325 | a = b; | |
326 | return std::make_pair(a, b); | |
327 | } | |
328 | ||
329 | if(boost::math::sign(fa) * boost::math::sign(fb) > 0) | |
330 | return boost::math::detail::pair_from_single(policies::raise_domain_error( | |
331 | function, | |
332 | "Parameters a and b do not bracket the root: a=%1%", a, pol)); | |
333 | // dummy value for fd, e and fe: | |
334 | fe = e = fd = 1e5F; | |
335 | ||
336 | if(fa != 0) | |
337 | { | |
338 | // | |
339 | // On the first step we take a secant step: | |
340 | // | |
341 | c = detail::secant_interpolate(a, b, fa, fb); | |
342 | detail::bracket(f, a, b, c, fa, fb, d, fd); | |
343 | --count; | |
344 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
345 | ||
346 | if(count && (fa != 0) && !tol(a, b)) | |
347 | { | |
348 | // | |
349 | // On the second step we take a quadratic interpolation: | |
350 | // | |
351 | c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2); | |
352 | e = d; | |
353 | fe = fd; | |
354 | detail::bracket(f, a, b, c, fa, fb, d, fd); | |
355 | --count; | |
356 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
357 | } | |
358 | } | |
359 | ||
360 | while(count && (fa != 0) && !tol(a, b)) | |
361 | { | |
362 | // save our brackets: | |
363 | a0 = a; | |
364 | b0 = b; | |
365 | // | |
366 | // Starting with the third step taken | |
367 | // we can use either quadratic or cubic interpolation. | |
368 | // Cubic interpolation requires that all four function values | |
369 | // fa, fb, fd, and fe are distinct, should that not be the case | |
370 | // then variable prof will get set to true, and we'll end up | |
371 | // taking a quadratic step instead. | |
372 | // | |
373 | T min_diff = tools::min_value<T>() * 32; | |
374 | bool prof = (fabs(fa - fb) < min_diff) || (fabs(fa - fd) < min_diff) || (fabs(fa - fe) < min_diff) || (fabs(fb - fd) < min_diff) || (fabs(fb - fe) < min_diff) || (fabs(fd - fe) < min_diff); | |
375 | if(prof) | |
376 | { | |
377 | c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2); | |
378 | BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!"); | |
379 | } | |
380 | else | |
381 | { | |
382 | c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe); | |
383 | } | |
384 | // | |
385 | // re-bracket, and check for termination: | |
386 | // | |
387 | e = d; | |
388 | fe = fd; | |
389 | detail::bracket(f, a, b, c, fa, fb, d, fd); | |
390 | if((0 == --count) || (fa == 0) || tol(a, b)) | |
391 | break; | |
392 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
393 | // | |
394 | // Now another interpolated step: | |
395 | // | |
396 | prof = (fabs(fa - fb) < min_diff) || (fabs(fa - fd) < min_diff) || (fabs(fa - fe) < min_diff) || (fabs(fb - fd) < min_diff) || (fabs(fb - fe) < min_diff) || (fabs(fd - fe) < min_diff); | |
397 | if(prof) | |
398 | { | |
399 | c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 3); | |
400 | BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!"); | |
401 | } | |
402 | else | |
403 | { | |
404 | c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe); | |
405 | } | |
406 | // | |
407 | // Bracket again, and check termination condition, update e: | |
408 | // | |
409 | detail::bracket(f, a, b, c, fa, fb, d, fd); | |
410 | if((0 == --count) || (fa == 0) || tol(a, b)) | |
411 | break; | |
412 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
413 | // | |
414 | // Now we take a double-length secant step: | |
415 | // | |
416 | if(fabs(fa) < fabs(fb)) | |
417 | { | |
418 | u = a; | |
419 | fu = fa; | |
420 | } | |
421 | else | |
422 | { | |
423 | u = b; | |
424 | fu = fb; | |
425 | } | |
426 | c = u - 2 * (fu / (fb - fa)) * (b - a); | |
427 | if(fabs(c - u) > (b - a) / 2) | |
428 | { | |
429 | c = a + (b - a) / 2; | |
430 | } | |
431 | // | |
432 | // Bracket again, and check termination condition: | |
433 | // | |
434 | e = d; | |
435 | fe = fd; | |
436 | detail::bracket(f, a, b, c, fa, fb, d, fd); | |
437 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
438 | BOOST_MATH_INSTRUMENT_CODE(" tol = " << T((fabs(a) - fabs(b)) / fabs(a))); | |
439 | if((0 == --count) || (fa == 0) || tol(a, b)) | |
440 | break; | |
441 | // | |
442 | // And finally... check to see if an additional bisection step is | |
443 | // to be taken, we do this if we're not converging fast enough: | |
444 | // | |
445 | if((b - a) < mu * (b0 - a0)) | |
446 | continue; | |
447 | // | |
448 | // bracket again on a bisection: | |
449 | // | |
450 | e = d; | |
451 | fe = fd; | |
452 | detail::bracket(f, a, b, T(a + (b - a) / 2), fa, fb, d, fd); | |
453 | --count; | |
454 | BOOST_MATH_INSTRUMENT_CODE("Not converging: Taking a bisection!!!!"); | |
455 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
456 | } // while loop | |
457 | ||
458 | max_iter -= count; | |
459 | if(fa == 0) | |
460 | { | |
461 | b = a; | |
462 | } | |
463 | else if(fb == 0) | |
464 | { | |
465 | a = b; | |
466 | } | |
467 | BOOST_MATH_LOG_COUNT(max_iter) | |
468 | return std::make_pair(a, b); | |
469 | } | |
470 | ||
471 | template <class F, class T, class Tol> | |
472 | inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, boost::uintmax_t& max_iter) | |
473 | { | |
474 | return toms748_solve(f, ax, bx, fax, fbx, tol, max_iter, policies::policy<>()); | |
475 | } | |
476 | ||
477 | template <class F, class T, class Tol, class Policy> | |
478 | inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, boost::uintmax_t& max_iter, const Policy& pol) | |
479 | { | |
480 | max_iter -= 2; | |
481 | std::pair<T, T> r = toms748_solve(f, ax, bx, f(ax), f(bx), tol, max_iter, pol); | |
482 | max_iter += 2; | |
483 | return r; | |
484 | } | |
485 | ||
486 | template <class F, class T, class Tol> | |
487 | inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, boost::uintmax_t& max_iter) | |
488 | { | |
489 | return toms748_solve(f, ax, bx, tol, max_iter, policies::policy<>()); | |
490 | } | |
491 | ||
492 | template <class F, class T, class Tol, class Policy> | |
493 | std::pair<T, T> bracket_and_solve_root(F f, const T& guess, T factor, bool rising, Tol tol, boost::uintmax_t& max_iter, const Policy& pol) | |
494 | { | |
495 | BOOST_MATH_STD_USING | |
496 | static const char* function = "boost::math::tools::bracket_and_solve_root<%1%>"; | |
497 | // | |
498 | // Set up inital brackets: | |
499 | // | |
500 | T a = guess; | |
501 | T b = a; | |
502 | T fa = f(a); | |
503 | T fb = fa; | |
504 | // | |
505 | // Set up invocation count: | |
506 | // | |
507 | boost::uintmax_t count = max_iter - 1; | |
508 | ||
509 | int step = 32; | |
510 | ||
511 | if((fa < 0) == (guess < 0 ? !rising : rising)) | |
512 | { | |
513 | // | |
514 | // Zero is to the right of b, so walk upwards | |
515 | // until we find it: | |
516 | // | |
517 | while((boost::math::sign)(fb) == (boost::math::sign)(fa)) | |
518 | { | |
519 | if(count == 0) | |
520 | return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol)); | |
521 | // | |
522 | // Heuristic: normally it's best not to increase the step sizes as we'll just end up | |
523 | // with a really wide range to search for the root. However, if the initial guess was *really* | |
524 | // bad then we need to speed up the search otherwise we'll take forever if we're orders of | |
525 | // magnitude out. This happens most often if the guess is a small value (say 1) and the result | |
526 | // we're looking for is close to std::numeric_limits<T>::min(). | |
527 | // | |
528 | if((max_iter - count) % step == 0) | |
529 | { | |
530 | factor *= 2; | |
531 | if(step > 1) step /= 2; | |
532 | } | |
533 | // | |
534 | // Now go ahead and move our guess by "factor": | |
535 | // | |
536 | a = b; | |
537 | fa = fb; | |
538 | b *= factor; | |
539 | fb = f(b); | |
540 | --count; | |
541 | BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); | |
542 | } | |
543 | } | |
544 | else | |
545 | { | |
546 | // | |
547 | // Zero is to the left of a, so walk downwards | |
548 | // until we find it: | |
549 | // | |
550 | while((boost::math::sign)(fb) == (boost::math::sign)(fa)) | |
551 | { | |
552 | if(fabs(a) < tools::min_value<T>()) | |
553 | { | |
554 | // Escape route just in case the answer is zero! | |
555 | max_iter -= count; | |
556 | max_iter += 1; | |
557 | return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0)); | |
558 | } | |
559 | if(count == 0) | |
560 | return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol)); | |
561 | // | |
562 | // Heuristic: normally it's best not to increase the step sizes as we'll just end up | |
563 | // with a really wide range to search for the root. However, if the initial guess was *really* | |
564 | // bad then we need to speed up the search otherwise we'll take forever if we're orders of | |
565 | // magnitude out. This happens most often if the guess is a small value (say 1) and the result | |
566 | // we're looking for is close to std::numeric_limits<T>::min(). | |
567 | // | |
568 | if((max_iter - count) % step == 0) | |
569 | { | |
570 | factor *= 2; | |
571 | if(step > 1) step /= 2; | |
572 | } | |
573 | // | |
574 | // Now go ahead and move are guess by "factor": | |
575 | // | |
576 | b = a; | |
577 | fb = fa; | |
578 | a /= factor; | |
579 | fa = f(a); | |
580 | --count; | |
581 | BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); | |
582 | } | |
583 | } | |
584 | max_iter -= count; | |
585 | max_iter += 1; | |
586 | std::pair<T, T> r = toms748_solve( | |
587 | f, | |
588 | (a < 0 ? b : a), | |
589 | (a < 0 ? a : b), | |
590 | (a < 0 ? fb : fa), | |
591 | (a < 0 ? fa : fb), | |
592 | tol, | |
593 | count, | |
594 | pol); | |
595 | max_iter += count; | |
596 | BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count); | |
597 | BOOST_MATH_LOG_COUNT(max_iter) | |
598 | return r; | |
599 | } | |
600 | ||
601 | template <class F, class T, class Tol> | |
602 | inline std::pair<T, T> bracket_and_solve_root(F f, const T& guess, const T& factor, bool rising, Tol tol, boost::uintmax_t& max_iter) | |
603 | { | |
604 | return bracket_and_solve_root(f, guess, factor, rising, tol, max_iter, policies::policy<>()); | |
605 | } | |
606 | ||
607 | } // namespace tools | |
608 | } // namespace math | |
609 | } // namespace boost | |
610 | ||
611 | ||
612 | #endif // BOOST_MATH_TOOLS_SOLVE_ROOT_HPP | |
613 |