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