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 <limits>
18#include <utility>
19#include <cstdint>
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
29namespace boost{ namespace math{ namespace tools{
30
31template <class T>
32class eps_tolerance
33{
34public:
35 eps_tolerance() : eps(4 * tools::epsilon<T>())
36 {
37
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 }
49private:
50 T eps;
51};
52
53struct equal_floor
54{
55 equal_floor()= default;
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
64struct equal_ceil
65{
66 equal_ceil()= default;
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
75struct equal_nearest_integer
76{
77 equal_nearest_integer()= default;
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
86namespace detail{
87
88template <class F, class T>
89void 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
152template <class T>
153inline 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
169template <class T>
170inline 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
190template <class T>
191T 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
246template <class T>
247T 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
295template <class F, class T, class Tol, class Policy>
296std::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)
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
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
312 std::uintmax_t count = max_iter;
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
478template <class F, class T, class Tol>
479inline 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)
480{
481 return toms748_solve(f, ax, bx, fax, fbx, tol, max_iter, policies::policy<>());
482}
483
484template <class F, class T, class Tol, class Policy>
485inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, std::uintmax_t& max_iter, const Policy& pol)
486{
487 if (max_iter <= 2)
488 return std::make_pair(ax, bx);
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
495template <class F, class T, class Tol>
496inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, std::uintmax_t& max_iter)
497{
498 return toms748_solve(f, ax, bx, tol, max_iter, policies::policy<>());
499}
500
501template <class F, class T, class Tol, class Policy>
502std::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)
503{
504 BOOST_MATH_STD_USING
505 static const char* function = "boost::math::tools::bracket_and_solve_root<%1%>";
506 //
507 // Set up initial brackets:
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 //
516 std::uintmax_t count = max_iter - 1;
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
610template <class F, class T, class Tol>
611inline 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)
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
623

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