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
28namespace boost{ namespace math{ namespace tools{
29
30template <class T>
31class eps_tolerance
32{
33public:
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 }
48private:
49 T eps;
50};
51
52struct 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
63struct 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
74struct 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
85namespace detail{
86
87template <class F, class T>
88void 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
151template <class T>
152inline 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
168template <class T>
169inline 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
189template <class T>
190T 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
245template <class T>
246T 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
294template <class F, class T, class Tol, class Policy>
295std::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 //
306 // Sanity check - are we allowed to iterate at all?
307 //
308 if (max_iter == 0)
309 return std::make_pair(ax, bx);
310
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
477template <class F, class T, class Tol>
478inline 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
483template <class F, class T, class Tol, class Policy>
484inline 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{
486 if (max_iter <= 2)
487 return std::make_pair(ax, bx);
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
494template <class F, class T, class Tol>
495inline 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
500template <class F, class T, class Tol, class Policy>
501std::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 initial 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
609template <class F, class T, class Tol>
610inline 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
622

source code of include/boost/math/tools/toms748_solve.hpp