]>
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 | ||
92f5a8d4 TL |
305 | // |
306 | // Sanity check - are we allowed to iterate at all? | |
307 | // | |
308 | if (max_iter == 0) | |
309 | return std::make_pair(ax, bx); | |
310 | ||
7c673cae FG |
311 | boost::uintmax_t count = max_iter; |
312 | T a, b, fa, fb, c, u, fu, a0, b0, d, fd, e, fe; | |
313 | static const T mu = 0.5f; | |
314 | ||
315 | // initialise a, b and fa, fb: | |
316 | a = ax; | |
317 | b = bx; | |
318 | if(a >= b) | |
319 | return boost::math::detail::pair_from_single(policies::raise_domain_error( | |
320 | function, | |
321 | "Parameters a and b out of order: a=%1%", a, pol)); | |
322 | fa = fax; | |
323 | fb = fbx; | |
324 | ||
325 | if(tol(a, b) || (fa == 0) || (fb == 0)) | |
326 | { | |
327 | max_iter = 0; | |
328 | if(fa == 0) | |
329 | b = a; | |
330 | else if(fb == 0) | |
331 | a = b; | |
332 | return std::make_pair(a, b); | |
333 | } | |
334 | ||
335 | if(boost::math::sign(fa) * boost::math::sign(fb) > 0) | |
336 | return boost::math::detail::pair_from_single(policies::raise_domain_error( | |
337 | function, | |
338 | "Parameters a and b do not bracket the root: a=%1%", a, pol)); | |
339 | // dummy value for fd, e and fe: | |
340 | fe = e = fd = 1e5F; | |
341 | ||
342 | if(fa != 0) | |
343 | { | |
344 | // | |
345 | // On the first step we take a secant step: | |
346 | // | |
347 | c = detail::secant_interpolate(a, b, fa, fb); | |
348 | detail::bracket(f, a, b, c, fa, fb, d, fd); | |
349 | --count; | |
350 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
351 | ||
352 | if(count && (fa != 0) && !tol(a, b)) | |
353 | { | |
354 | // | |
355 | // On the second step we take a quadratic interpolation: | |
356 | // | |
357 | c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2); | |
358 | e = d; | |
359 | fe = fd; | |
360 | detail::bracket(f, a, b, c, fa, fb, d, fd); | |
361 | --count; | |
362 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
363 | } | |
364 | } | |
365 | ||
366 | while(count && (fa != 0) && !tol(a, b)) | |
367 | { | |
368 | // save our brackets: | |
369 | a0 = a; | |
370 | b0 = b; | |
371 | // | |
372 | // Starting with the third step taken | |
373 | // we can use either quadratic or cubic interpolation. | |
374 | // Cubic interpolation requires that all four function values | |
375 | // fa, fb, fd, and fe are distinct, should that not be the case | |
376 | // then variable prof will get set to true, and we'll end up | |
377 | // taking a quadratic step instead. | |
378 | // | |
379 | T min_diff = tools::min_value<T>() * 32; | |
380 | 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); | |
381 | if(prof) | |
382 | { | |
383 | c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2); | |
384 | BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!"); | |
385 | } | |
386 | else | |
387 | { | |
388 | c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe); | |
389 | } | |
390 | // | |
391 | // re-bracket, and check for termination: | |
392 | // | |
393 | e = d; | |
394 | fe = fd; | |
395 | detail::bracket(f, a, b, c, fa, fb, d, fd); | |
396 | if((0 == --count) || (fa == 0) || tol(a, b)) | |
397 | break; | |
398 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
399 | // | |
400 | // Now another interpolated step: | |
401 | // | |
402 | 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); | |
403 | if(prof) | |
404 | { | |
405 | c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 3); | |
406 | BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!"); | |
407 | } | |
408 | else | |
409 | { | |
410 | c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe); | |
411 | } | |
412 | // | |
413 | // Bracket again, and check termination condition, update e: | |
414 | // | |
415 | detail::bracket(f, a, b, c, fa, fb, d, fd); | |
416 | if((0 == --count) || (fa == 0) || tol(a, b)) | |
417 | break; | |
418 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
419 | // | |
420 | // Now we take a double-length secant step: | |
421 | // | |
422 | if(fabs(fa) < fabs(fb)) | |
423 | { | |
424 | u = a; | |
425 | fu = fa; | |
426 | } | |
427 | else | |
428 | { | |
429 | u = b; | |
430 | fu = fb; | |
431 | } | |
432 | c = u - 2 * (fu / (fb - fa)) * (b - a); | |
433 | if(fabs(c - u) > (b - a) / 2) | |
434 | { | |
435 | c = a + (b - a) / 2; | |
436 | } | |
437 | // | |
438 | // Bracket again, and check termination condition: | |
439 | // | |
440 | e = d; | |
441 | fe = fd; | |
442 | detail::bracket(f, a, b, c, fa, fb, d, fd); | |
443 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
444 | BOOST_MATH_INSTRUMENT_CODE(" tol = " << T((fabs(a) - fabs(b)) / fabs(a))); | |
445 | if((0 == --count) || (fa == 0) || tol(a, b)) | |
446 | break; | |
447 | // | |
448 | // And finally... check to see if an additional bisection step is | |
449 | // to be taken, we do this if we're not converging fast enough: | |
450 | // | |
451 | if((b - a) < mu * (b0 - a0)) | |
452 | continue; | |
453 | // | |
454 | // bracket again on a bisection: | |
455 | // | |
456 | e = d; | |
457 | fe = fd; | |
458 | detail::bracket(f, a, b, T(a + (b - a) / 2), fa, fb, d, fd); | |
459 | --count; | |
460 | BOOST_MATH_INSTRUMENT_CODE("Not converging: Taking a bisection!!!!"); | |
461 | BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b); | |
462 | } // while loop | |
463 | ||
464 | max_iter -= count; | |
465 | if(fa == 0) | |
466 | { | |
467 | b = a; | |
468 | } | |
469 | else if(fb == 0) | |
470 | { | |
471 | a = b; | |
472 | } | |
473 | BOOST_MATH_LOG_COUNT(max_iter) | |
474 | return std::make_pair(a, b); | |
475 | } | |
476 | ||
477 | template <class F, class T, class Tol> | |
478 | 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) | |
479 | { | |
480 | return toms748_solve(f, ax, bx, fax, fbx, tol, max_iter, policies::policy<>()); | |
481 | } | |
482 | ||
483 | template <class F, class T, class Tol, class Policy> | |
484 | 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) | |
485 | { | |
92f5a8d4 TL |
486 | if (max_iter <= 2) |
487 | return std::make_pair(ax, bx); | |
7c673cae FG |
488 | max_iter -= 2; |
489 | std::pair<T, T> r = toms748_solve(f, ax, bx, f(ax), f(bx), tol, max_iter, pol); | |
490 | max_iter += 2; | |
491 | return r; | |
492 | } | |
493 | ||
494 | template <class F, class T, class Tol> | |
495 | inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, boost::uintmax_t& max_iter) | |
496 | { | |
497 | return toms748_solve(f, ax, bx, tol, max_iter, policies::policy<>()); | |
498 | } | |
499 | ||
500 | template <class F, class T, class Tol, class Policy> | |
501 | 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) | |
502 | { | |
503 | BOOST_MATH_STD_USING | |
504 | static const char* function = "boost::math::tools::bracket_and_solve_root<%1%>"; | |
505 | // | |
506 | // Set up inital brackets: | |
507 | // | |
508 | T a = guess; | |
509 | T b = a; | |
510 | T fa = f(a); | |
511 | T fb = fa; | |
512 | // | |
513 | // Set up invocation count: | |
514 | // | |
515 | boost::uintmax_t count = max_iter - 1; | |
516 | ||
517 | int step = 32; | |
518 | ||
519 | if((fa < 0) == (guess < 0 ? !rising : rising)) | |
520 | { | |
521 | // | |
522 | // Zero is to the right of b, so walk upwards | |
523 | // until we find it: | |
524 | // | |
525 | while((boost::math::sign)(fb) == (boost::math::sign)(fa)) | |
526 | { | |
527 | if(count == 0) | |
528 | return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol)); | |
529 | // | |
530 | // Heuristic: normally it's best not to increase the step sizes as we'll just end up | |
531 | // with a really wide range to search for the root. However, if the initial guess was *really* | |
532 | // bad then we need to speed up the search otherwise we'll take forever if we're orders of | |
533 | // magnitude out. This happens most often if the guess is a small value (say 1) and the result | |
534 | // we're looking for is close to std::numeric_limits<T>::min(). | |
535 | // | |
536 | if((max_iter - count) % step == 0) | |
537 | { | |
538 | factor *= 2; | |
539 | if(step > 1) step /= 2; | |
540 | } | |
541 | // | |
542 | // Now go ahead and move our guess by "factor": | |
543 | // | |
544 | a = b; | |
545 | fa = fb; | |
546 | b *= factor; | |
547 | fb = f(b); | |
548 | --count; | |
549 | BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); | |
550 | } | |
551 | } | |
552 | else | |
553 | { | |
554 | // | |
555 | // Zero is to the left of a, so walk downwards | |
556 | // until we find it: | |
557 | // | |
558 | while((boost::math::sign)(fb) == (boost::math::sign)(fa)) | |
559 | { | |
560 | if(fabs(a) < tools::min_value<T>()) | |
561 | { | |
562 | // Escape route just in case the answer is zero! | |
563 | max_iter -= count; | |
564 | max_iter += 1; | |
565 | return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0)); | |
566 | } | |
567 | if(count == 0) | |
568 | return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol)); | |
569 | // | |
570 | // Heuristic: normally it's best not to increase the step sizes as we'll just end up | |
571 | // with a really wide range to search for the root. However, if the initial guess was *really* | |
572 | // bad then we need to speed up the search otherwise we'll take forever if we're orders of | |
573 | // magnitude out. This happens most often if the guess is a small value (say 1) and the result | |
574 | // we're looking for is close to std::numeric_limits<T>::min(). | |
575 | // | |
576 | if((max_iter - count) % step == 0) | |
577 | { | |
578 | factor *= 2; | |
579 | if(step > 1) step /= 2; | |
580 | } | |
581 | // | |
582 | // Now go ahead and move are guess by "factor": | |
583 | // | |
584 | b = a; | |
585 | fb = fa; | |
586 | a /= factor; | |
587 | fa = f(a); | |
588 | --count; | |
589 | BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count); | |
590 | } | |
591 | } | |
592 | max_iter -= count; | |
593 | max_iter += 1; | |
594 | std::pair<T, T> r = toms748_solve( | |
595 | f, | |
596 | (a < 0 ? b : a), | |
597 | (a < 0 ? a : b), | |
598 | (a < 0 ? fb : fa), | |
599 | (a < 0 ? fa : fb), | |
600 | tol, | |
601 | count, | |
602 | pol); | |
603 | max_iter += count; | |
604 | BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count); | |
605 | BOOST_MATH_LOG_COUNT(max_iter) | |
606 | return r; | |
607 | } | |
608 | ||
609 | template <class F, class T, class Tol> | |
610 | 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) | |
611 | { | |
612 | return bracket_and_solve_root(f, guess, factor, rising, tol, max_iter, policies::policy<>()); | |
613 | } | |
614 | ||
615 | } // namespace tools | |
616 | } // namespace math | |
617 | } // namespace boost | |
618 | ||
619 | ||
620 | #endif // BOOST_MATH_TOOLS_SOLVE_ROOT_HPP | |
621 |