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