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_NEWTON_SOLVER_HPP
7#define BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP
8
9#ifdef _MSC_VER
10#pragma once
11#endif
12#include <boost/math/tools/complex.hpp> // test for multiprecision types in complex Newton
13
14#include <utility>
15#include <cmath>
16#include <tuple>
17#include <cstdint>
18
19#include <boost/math/tools/config.hpp>
20#include <boost/math/tools/cxx03_warn.hpp>
21
22#include <boost/math/special_functions/sign.hpp>
23#include <boost/math/special_functions/next.hpp>
24#include <boost/math/tools/toms748_solve.hpp>
25#include <boost/math/policies/error_handling.hpp>
26
27namespace boost {
28namespace math {
29namespace tools {
30
31namespace detail {
32
33namespace dummy {
34
35 template<int n, class T>
36 typename T::value_type get(const T&) BOOST_MATH_NOEXCEPT(T);
37}
38
39template <class Tuple, class T>
40void unpack_tuple(const Tuple& t, T& a, T& b) BOOST_MATH_NOEXCEPT(T)
41{
42 using dummy::get;
43 // Use ADL to find the right overload for get:
44 a = get<0>(t);
45 b = get<1>(t);
46}
47template <class Tuple, class T>
48void unpack_tuple(const Tuple& t, T& a, T& b, T& c) BOOST_MATH_NOEXCEPT(T)
49{
50 using dummy::get;
51 // Use ADL to find the right overload for get:
52 a = get<0>(t);
53 b = get<1>(t);
54 c = get<2>(t);
55}
56
57template <class Tuple, class T>
58inline void unpack_0(const Tuple& t, T& val) BOOST_MATH_NOEXCEPT(T)
59{
60 using dummy::get;
61 // Rely on ADL to find the correct overload of get:
62 val = get<0>(t);
63}
64
65template <class T, class U, class V>
66inline void unpack_tuple(const std::pair<T, U>& p, V& a, V& b) BOOST_MATH_NOEXCEPT(T)
67{
68 a = p.first;
69 b = p.second;
70}
71template <class T, class U, class V>
72inline void unpack_0(const std::pair<T, U>& p, V& a) BOOST_MATH_NOEXCEPT(T)
73{
74 a = p.first;
75}
76
77template <class F, class T>
78void handle_zero_derivative(F f,
79 T& last_f0,
80 const T& f0,
81 T& delta,
82 T& result,
83 T& guess,
84 const T& min,
85 const T& max) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
86{
87 if (last_f0 == 0)
88 {
89 // this must be the first iteration, pretend that we had a
90 // previous one at either min or max:
91 if (result == min)
92 {
93 guess = max;
94 }
95 else
96 {
97 guess = min;
98 }
99 unpack_0(f(guess), last_f0);
100 delta = guess - result;
101 }
102 if (sign(last_f0) * sign(f0) < 0)
103 {
104 // we've crossed over so move in opposite direction to last step:
105 if (delta < 0)
106 {
107 delta = (result - min) / 2;
108 }
109 else
110 {
111 delta = (result - max) / 2;
112 }
113 }
114 else
115 {
116 // move in same direction as last step:
117 if (delta < 0)
118 {
119 delta = (result - max) / 2;
120 }
121 else
122 {
123 delta = (result - min) / 2;
124 }
125 }
126}
127
128} // namespace
129
130template <class F, class T, class Tol, class Policy>
131std::pair<T, T> bisect(F f, T min, T max, Tol tol, std::uintmax_t& max_iter, const Policy& pol) noexcept(policies::is_noexcept_error_policy<Policy>::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
132{
133 T fmin = f(min);
134 T fmax = f(max);
135 if (fmin == 0)
136 {
137 max_iter = 2;
138 return std::make_pair(min, min);
139 }
140 if (fmax == 0)
141 {
142 max_iter = 2;
143 return std::make_pair(max, max);
144 }
145
146 //
147 // Error checking:
148 //
149 static const char* function = "boost::math::tools::bisect<%1%>";
150 if (min >= max)
151 {
152 return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function,
153 "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol));
154 }
155 if (fmin * fmax >= 0)
156 {
157 return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function,
158 "No change of sign in boost::math::tools::bisect, either there is no root to find, or there are multiple roots in the interval (f(min) = %1%).", fmin, pol));
159 }
160
161 //
162 // Three function invocations so far:
163 //
164 std::uintmax_t count = max_iter;
165 if (count < 3)
166 count = 0;
167 else
168 count -= 3;
169
170 while (count && (0 == tol(min, max)))
171 {
172 T mid = (min + max) / 2;
173 T fmid = f(mid);
174 if ((mid == max) || (mid == min))
175 break;
176 if (fmid == 0)
177 {
178 min = max = mid;
179 break;
180 }
181 else if (sign(fmid) * sign(fmin) < 0)
182 {
183 max = mid;
184 }
185 else
186 {
187 min = mid;
188 fmin = fmid;
189 }
190 --count;
191 }
192
193 max_iter -= count;
194
195#ifdef BOOST_MATH_INSTRUMENT
196 std::cout << "Bisection required " << max_iter << " iterations.\n";
197#endif
198
199 return std::make_pair(min, max);
200}
201
202template <class F, class T, class Tol>
203inline std::pair<T, T> bisect(F f, T min, T max, Tol tol, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
204{
205 return bisect(f, min, max, tol, max_iter, policies::policy<>());
206}
207
208template <class F, class T, class Tol>
209inline std::pair<T, T> bisect(F f, T min, T max, Tol tol) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
210{
211 std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)();
212 return bisect(f, min, max, tol, m, policies::policy<>());
213}
214
215
216template <class F, class T>
217T newton_raphson_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
218{
219 BOOST_MATH_STD_USING
220
221 static const char* function = "boost::math::tools::newton_raphson_iterate<%1%>";
222 if (min > max)
223 {
224 return policies::raise_evaluation_error(function, "Range arguments in wrong order in boost::math::tools::newton_raphson_iterate(first arg=%1%)", min, boost::math::policies::policy<>());
225 }
226
227 T f0(0), f1, last_f0(0);
228 T result = guess;
229
230 T factor = static_cast<T>(ldexp(1.0, 1 - digits));
231 T delta = tools::max_value<T>();
232 T delta1 = tools::max_value<T>();
233 T delta2 = tools::max_value<T>();
234
235 //
236 // We use these to sanity check that we do actually bracket a root,
237 // we update these to the function value when we update the endpoints
238 // of the range. Then, provided at some point we update both endpoints
239 // checking that max_range_f * min_range_f <= 0 verifies there is a root
240 // to be found somewhere. Note that if there is no root, and we approach
241 // a local minima, then the derivative will go to zero, and hence the next
242 // step will jump out of bounds (or at least past the minima), so this
243 // check *should* happen in pathological cases.
244 //
245 T max_range_f = 0;
246 T min_range_f = 0;
247
248 std::uintmax_t count(max_iter);
249
250#ifdef BOOST_MATH_INSTRUMENT
251 std::cout << "Newton_raphson_iterate, guess = " << guess << ", min = " << min << ", max = " << max
252 << ", digits = " << digits << ", max_iter = " << max_iter << "\n";
253#endif
254
255 do {
256 last_f0 = f0;
257 delta2 = delta1;
258 delta1 = delta;
259 detail::unpack_tuple(f(result), f0, f1);
260 --count;
261 if (0 == f0)
262 break;
263 if (f1 == 0)
264 {
265 // Oops zero derivative!!!
266 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
267 }
268 else
269 {
270 delta = f0 / f1;
271 }
272#ifdef BOOST_MATH_INSTRUMENT
273 std::cout << "Newton iteration " << max_iter - count << ", delta = " << delta << ", residual = " << f0 << "\n";
274#endif
275 if (fabs(delta * 2) > fabs(delta2))
276 {
277 // Last two steps haven't converged.
278 T shift = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
279 if ((result != 0) && (fabs(shift) > fabs(result)))
280 {
281 delta = sign(delta) * fabs(result) * 1.1f; // Protect against huge jumps!
282 //delta = sign(delta) * result; // Protect against huge jumps! Failed for negative result. https://github.com/boostorg/math/issues/216
283 }
284 else
285 delta = shift;
286 // reset delta1/2 so we don't take this branch next time round:
287 delta1 = 3 * delta;
288 delta2 = 3 * delta;
289 }
290 guess = result;
291 result -= delta;
292 if (result <= min)
293 {
294 delta = 0.5F * (guess - min);
295 result = guess - delta;
296 if ((result == min) || (result == max))
297 break;
298 }
299 else if (result >= max)
300 {
301 delta = 0.5F * (guess - max);
302 result = guess - delta;
303 if ((result == min) || (result == max))
304 break;
305 }
306 // Update brackets:
307 if (delta > 0)
308 {
309 max = guess;
310 max_range_f = f0;
311 }
312 else
313 {
314 min = guess;
315 min_range_f = f0;
316 }
317 //
318 // Sanity check that we bracket the root:
319 //
320 if (max_range_f * min_range_f > 0)
321 {
322 return policies::raise_evaluation_error(function, "There appears to be no root to be found in boost::math::tools::newton_raphson_iterate, perhaps we have a local minima near current best guess of %1%", guess, boost::math::policies::policy<>());
323 }
324 }while(count && (fabs(result * factor) < fabs(delta)));
325
326 max_iter -= count;
327
328#ifdef BOOST_MATH_INSTRUMENT
329 std::cout << "Newton Raphson required " << max_iter << " iterations\n";
330#endif
331
332 return result;
333}
334
335template <class F, class T>
336inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
337{
338 std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)();
339 return newton_raphson_iterate(f, guess, min, max, digits, m);
340}
341
342namespace detail {
343
344 struct halley_step
345 {
346 template <class T>
347 static T step(const T& /*x*/, const T& f0, const T& f1, const T& f2) noexcept(BOOST_MATH_IS_FLOAT(T))
348 {
349 using std::fabs;
350 T denom = 2 * f0;
351 T num = 2 * f1 - f0 * (f2 / f1);
352 T delta;
353
354 BOOST_MATH_INSTRUMENT_VARIABLE(denom);
355 BOOST_MATH_INSTRUMENT_VARIABLE(num);
356
357 if ((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>()))
358 {
359 // possible overflow, use Newton step:
360 delta = f0 / f1;
361 }
362 else
363 delta = denom / num;
364 return delta;
365 }
366 };
367
368 template <class F, class T>
369 T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())));
370
371 template <class F, class T>
372 T bracket_root_towards_max(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
373 {
374 using std::fabs;
375 using std::ldexp;
376 using std::abs;
377 using std::frexp;
378 if(count < 2)
379 return guess - (max + min) / 2; // Not enough counts left to do anything!!
380 //
381 // Move guess towards max until we bracket the root, updating min and max as we go:
382 //
383 int e;
384 frexp(max / guess, &e);
385 e = abs(e);
386 T guess0 = guess;
387 T multiplier = e < 64 ? static_cast<T>(2) : static_cast<T>(ldexp(T(1), e / 32));
388 T f_current = f0;
389 if (fabs(min) < fabs(max))
390 {
391 while (--count && ((f_current < 0) == (f0 < 0)))
392 {
393 min = guess;
394 guess *= multiplier;
395 if (guess > max)
396 {
397 guess = max;
398 f_current = -f_current; // There must be a change of sign!
399 break;
400 }
401 multiplier *= e > 1024 ? 8 : 2;
402 unpack_0(f(guess), f_current);
403 }
404 }
405 else
406 {
407 //
408 // If min and max are negative we have to divide to head towards max:
409 //
410 while (--count && ((f_current < 0) == (f0 < 0)))
411 {
412 min = guess;
413 guess /= multiplier;
414 if (guess > max)
415 {
416 guess = max;
417 f_current = -f_current; // There must be a change of sign!
418 break;
419 }
420 multiplier *= e > 1024 ? 8 : 2;
421 unpack_0(f(guess), f_current);
422 }
423 }
424
425 if (count)
426 {
427 max = guess;
428 if (multiplier > 16)
429 return (guess0 - guess) + bracket_root_towards_min(f, guess, f_current, min, max, count);
430 }
431 return guess0 - (max + min) / 2;
432 }
433
434 template <class F, class T>
435 T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
436 {
437 using std::fabs;
438 using std::ldexp;
439 using std::abs;
440 using std::frexp;
441 if (count < 2)
442 return guess - (max + min) / 2; // Not enough counts left to do anything!!
443 //
444 // Move guess towards min until we bracket the root, updating min and max as we go:
445 //
446 int e;
447 frexp(guess / min, &e);
448 e = abs(e);
449 T guess0 = guess;
450 T multiplier = e < 64 ? static_cast<T>(2) : static_cast<T>(ldexp(T(1), e / 32));
451 T f_current = f0;
452
453 if (fabs(min) < fabs(max))
454 {
455 while (--count && ((f_current < 0) == (f0 < 0)))
456 {
457 max = guess;
458 guess /= multiplier;
459 if (guess < min)
460 {
461 guess = min;
462 f_current = -f_current; // There must be a change of sign!
463 break;
464 }
465 multiplier *= e > 1024 ? 8 : 2;
466 unpack_0(f(guess), f_current);
467 }
468 }
469 else
470 {
471 //
472 // If min and max are negative we have to multiply to head towards min:
473 //
474 while (--count && ((f_current < 0) == (f0 < 0)))
475 {
476 max = guess;
477 guess *= multiplier;
478 if (guess < min)
479 {
480 guess = min;
481 f_current = -f_current; // There must be a change of sign!
482 break;
483 }
484 multiplier *= e > 1024 ? 8 : 2;
485 unpack_0(f(guess), f_current);
486 }
487 }
488
489 if (count)
490 {
491 min = guess;
492 if (multiplier > 16)
493 return (guess0 - guess) + bracket_root_towards_max(f, guess, f_current, min, max, count);
494 }
495 return guess0 - (max + min) / 2;
496 }
497
498 template <class Stepper, class F, class T>
499 T second_order_root_finder(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
500 {
501 BOOST_MATH_STD_USING
502
503#ifdef BOOST_MATH_INSTRUMENT
504 std::cout << "Second order root iteration, guess = " << guess << ", min = " << min << ", max = " << max
505 << ", digits = " << digits << ", max_iter = " << max_iter << "\n";
506#endif
507 static const char* function = "boost::math::tools::halley_iterate<%1%>";
508 if (min >= max)
509 {
510 return policies::raise_evaluation_error(function, "Range arguments in wrong order in boost::math::tools::halley_iterate(first arg=%1%)", min, boost::math::policies::policy<>());
511 }
512
513 T f0(0), f1, f2;
514 T result = guess;
515
516 T factor = ldexp(static_cast<T>(1.0), 1 - digits);
517 T delta = (std::max)(T(10000000 * guess), T(10000000)); // arbitrarily large delta
518 T last_f0 = 0;
519 T delta1 = delta;
520 T delta2 = delta;
521 bool out_of_bounds_sentry = false;
522
523 #ifdef BOOST_MATH_INSTRUMENT
524 std::cout << "Second order root iteration, limit = " << factor << "\n";
525 #endif
526
527 //
528 // We use these to sanity check that we do actually bracket a root,
529 // we update these to the function value when we update the endpoints
530 // of the range. Then, provided at some point we update both endpoints
531 // checking that max_range_f * min_range_f <= 0 verifies there is a root
532 // to be found somewhere. Note that if there is no root, and we approach
533 // a local minima, then the derivative will go to zero, and hence the next
534 // step will jump out of bounds (or at least past the minima), so this
535 // check *should* happen in pathological cases.
536 //
537 T max_range_f = 0;
538 T min_range_f = 0;
539
540 std::uintmax_t count(max_iter);
541
542 do {
543 last_f0 = f0;
544 delta2 = delta1;
545 delta1 = delta;
546#ifndef BOOST_NO_EXCEPTIONS
547 try
548#endif
549 {
550 detail::unpack_tuple(f(result), f0, f1, f2);
551 }
552#ifndef BOOST_NO_EXCEPTIONS
553 catch (const std::overflow_error&)
554 {
555 f0 = max > 0 ? tools::max_value<T>() : -tools::min_value<T>();
556 f1 = f2 = 0;
557 }
558#endif
559 --count;
560
561 BOOST_MATH_INSTRUMENT_VARIABLE(f0);
562 BOOST_MATH_INSTRUMENT_VARIABLE(f1);
563 BOOST_MATH_INSTRUMENT_VARIABLE(f2);
564
565 if (0 == f0)
566 break;
567 if (f1 == 0)
568 {
569 // Oops zero derivative!!!
570 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
571 }
572 else
573 {
574 if (f2 != 0)
575 {
576 delta = Stepper::step(result, f0, f1, f2);
577 if (delta * f1 / f0 < 0)
578 {
579 // Oh dear, we have a problem as Newton and Halley steps
580 // disagree about which way we should move. Probably
581 // there is cancelation error in the calculation of the
582 // Halley step, or else the derivatives are so small
583 // that their values are basically trash. We will move
584 // in the direction indicated by a Newton step, but
585 // by no more than twice the current guess value, otherwise
586 // we can jump way out of bounds if we're not careful.
587 // See https://svn.boost.org/trac/boost/ticket/8314.
588 delta = f0 / f1;
589 if (fabs(delta) > 2 * fabs(result))
590 delta = (delta < 0 ? -1 : 1) * 2 * fabs(result);
591 }
592 }
593 else
594 delta = f0 / f1;
595 }
596 #ifdef BOOST_MATH_INSTRUMENT
597 std::cout << "Second order root iteration, delta = " << delta << ", residual = " << f0 << "\n";
598 #endif
599 T convergence = fabs(delta / delta2);
600 if ((convergence > 0.8) && (convergence < 2))
601 {
602 // last two steps haven't converged.
603 if (fabs(min) < 1 ? fabs(1000 * min) < fabs(max) : fabs(max / min) > 1000)
604 {
605 if(delta > 0)
606 delta = bracket_root_towards_min(f, result, f0, min, max, count);
607 else
608 delta = bracket_root_towards_max(f, result, f0, min, max, count);
609 }
610 else
611 {
612 delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
613 if ((result != 0) && (fabs(delta) > result))
614 delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps!
615 }
616 // reset delta2 so that this branch will *not* be taken on the
617 // next iteration:
618 delta2 = delta * 3;
619 delta1 = delta * 3;
620 BOOST_MATH_INSTRUMENT_VARIABLE(delta);
621 }
622 guess = result;
623 result -= delta;
624 BOOST_MATH_INSTRUMENT_VARIABLE(result);
625
626 // check for out of bounds step:
627 if (result < min)
628 {
629 T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min)))
630 ? T(1000)
631 : (fabs(min) < 1) && (fabs(tools::max_value<T>() * min) < fabs(result))
632 ? ((min < 0) != (result < 0)) ? -tools::max_value<T>() : tools::max_value<T>() : T(result / min);
633 if (fabs(diff) < 1)
634 diff = 1 / diff;
635 if (!out_of_bounds_sentry && (diff > 0) && (diff < 3))
636 {
637 // Only a small out of bounds step, lets assume that the result
638 // is probably approximately at min:
639 delta = 0.99f * (guess - min);
640 result = guess - delta;
641 out_of_bounds_sentry = true; // only take this branch once!
642 }
643 else
644 {
645 if (fabs(float_distance(min, max)) < 2)
646 {
647 result = guess = (min + max) / 2;
648 break;
649 }
650 delta = bracket_root_towards_min(f, guess, f0, min, max, count);
651 result = guess - delta;
652 if (result <= min)
653 result = float_next(min);
654 if (result >= max)
655 result = float_prior(max);
656 guess = min;
657 continue;
658 }
659 }
660 else if (result > max)
661 {
662 T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max);
663 if (fabs(diff) < 1)
664 diff = 1 / diff;
665 if (!out_of_bounds_sentry && (diff > 0) && (diff < 3))
666 {
667 // Only a small out of bounds step, lets assume that the result
668 // is probably approximately at min:
669 delta = 0.99f * (guess - max);
670 result = guess - delta;
671 out_of_bounds_sentry = true; // only take this branch once!
672 }
673 else
674 {
675 if (fabs(float_distance(min, max)) < 2)
676 {
677 result = guess = (min + max) / 2;
678 break;
679 }
680 delta = bracket_root_towards_max(f, guess, f0, min, max, count);
681 result = guess - delta;
682 if (result >= max)
683 result = float_prior(max);
684 if (result <= min)
685 result = float_next(min);
686 guess = min;
687 continue;
688 }
689 }
690 // update brackets:
691 if (delta > 0)
692 {
693 max = guess;
694 max_range_f = f0;
695 }
696 else
697 {
698 min = guess;
699 min_range_f = f0;
700 }
701 //
702 // Sanity check that we bracket the root:
703 //
704 if (max_range_f * min_range_f > 0)
705 {
706 return policies::raise_evaluation_error(function, "There appears to be no root to be found in boost::math::tools::newton_raphson_iterate, perhaps we have a local minima near current best guess of %1%", guess, boost::math::policies::policy<>());
707 }
708 } while(count && (fabs(result * factor) < fabs(delta)));
709
710 max_iter -= count;
711
712 #ifdef BOOST_MATH_INSTRUMENT
713 std::cout << "Second order root finder required " << max_iter << " iterations.\n";
714 #endif
715
716 return result;
717 }
718} // T second_order_root_finder
719
720template <class F, class T>
721T halley_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
722{
723 return detail::second_order_root_finder<detail::halley_step>(f, guess, min, max, digits, max_iter);
724}
725
726template <class F, class T>
727inline T halley_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
728{
729 std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)();
730 return halley_iterate(f, guess, min, max, digits, m);
731}
732
733namespace detail {
734
735 struct schroder_stepper
736 {
737 template <class T>
738 static T step(const T& x, const T& f0, const T& f1, const T& f2) noexcept(BOOST_MATH_IS_FLOAT(T))
739 {
740 using std::fabs;
741 T ratio = f0 / f1;
742 T delta;
743 if ((x != 0) && (fabs(ratio / x) < 0.1))
744 {
745 delta = ratio + (f2 / (2 * f1)) * ratio * ratio;
746 // check second derivative doesn't over compensate:
747 if (delta * ratio < 0)
748 delta = ratio;
749 }
750 else
751 delta = ratio; // fall back to Newton iteration.
752 return delta;
753 }
754 };
755
756}
757
758template <class F, class T>
759T schroder_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
760{
761 return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);
762}
763
764template <class F, class T>
765inline T schroder_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
766{
767 std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)();
768 return schroder_iterate(f, guess, min, max, digits, m);
769}
770//
771// These two are the old spelling of this function, retained for backwards compatibility just in case:
772//
773template <class F, class T>
774T schroeder_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
775{
776 return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);
777}
778
779template <class F, class T>
780inline T schroeder_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
781{
782 std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)();
783 return schroder_iterate(f, guess, min, max, digits, m);
784}
785
786#ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS
787/*
788 * Why do we set the default maximum number of iterations to the number of digits in the type?
789 * Because for double roots, the number of digits increases linearly with the number of iterations,
790 * so this default should recover full precision even in this somewhat pathological case.
791 * For isolated roots, the problem is so rapidly convergent that this doesn't matter at all.
792 */
793template<class Complex, class F>
794Complex complex_newton(F g, Complex guess, int max_iterations = std::numeric_limits<typename Complex::value_type>::digits)
795{
796 typedef typename Complex::value_type Real;
797 using std::norm;
798 using std::abs;
799 using std::max;
800 // z0, z1, and z2 cannot be the same, in case we immediately need to resort to Muller's Method:
801 Complex z0 = guess + Complex(1, 0);
802 Complex z1 = guess + Complex(0, 1);
803 Complex z2 = guess;
804
805 do {
806 auto pair = g(z2);
807 if (norm(pair.second) == 0)
808 {
809 // Muller's method. Notation follows Numerical Recipes, 9.5.2:
810 Complex q = (z2 - z1) / (z1 - z0);
811 auto P0 = g(z0);
812 auto P1 = g(z1);
813 Complex qp1 = static_cast<Complex>(1) + q;
814 Complex A = q * (pair.first - qp1 * P1.first + q * P0.first);
815
816 Complex B = (static_cast<Complex>(2) * q + static_cast<Complex>(1)) * pair.first - qp1 * qp1 * P1.first + q * q * P0.first;
817 Complex C = qp1 * pair.first;
818 Complex rad = sqrt(B * B - static_cast<Complex>(4) * A * C);
819 Complex denom1 = B + rad;
820 Complex denom2 = B - rad;
821 Complex correction = (z1 - z2) * static_cast<Complex>(2) * C;
822 if (norm(denom1) > norm(denom2))
823 {
824 correction /= denom1;
825 }
826 else
827 {
828 correction /= denom2;
829 }
830
831 z0 = z1;
832 z1 = z2;
833 z2 = z2 + correction;
834 }
835 else
836 {
837 z0 = z1;
838 z1 = z2;
839 z2 = z2 - (pair.first / pair.second);
840 }
841
842 // See: https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root
843 // If f' is continuous, then convergence of x_n -> x* implies f(x*) = 0.
844 // This condition approximates this convergence condition by requiring three consecutive iterates to be clustered.
845 Real tol = (max)(abs(z2) * std::numeric_limits<Real>::epsilon(), std::numeric_limits<Real>::epsilon());
846 bool real_close = abs(z0.real() - z1.real()) < tol && abs(z0.real() - z2.real()) < tol && abs(z1.real() - z2.real()) < tol;
847 bool imag_close = abs(z0.imag() - z1.imag()) < tol && abs(z0.imag() - z2.imag()) < tol && abs(z1.imag() - z2.imag()) < tol;
848 if (real_close && imag_close)
849 {
850 return z2;
851 }
852
853 } while (max_iterations--);
854
855 // The idea is that if we can get abs(f) < eps, we should, but if we go through all these iterations
856 // and abs(f) < sqrt(eps), then roundoff error simply does not allow that we can evaluate f to < eps
857 // This is somewhat awkward as it isn't scale invariant, but using the Daubechies coefficient example code,
858 // I found this condition generates correct roots, whereas the scale invariant condition discussed here:
859 // https://scicomp.stackexchange.com/questions/30597/defining-a-condition-number-and-termination-criteria-for-newtons-method
860 // allows nonroots to be passed off as roots.
861 auto pair = g(z2);
862 if (abs(pair.first) < sqrt(std::numeric_limits<Real>::epsilon()))
863 {
864 return z2;
865 }
866
867 return { std::numeric_limits<Real>::quiet_NaN(),
868 std::numeric_limits<Real>::quiet_NaN() };
869}
870#endif
871
872
873#if !defined(BOOST_NO_CXX17_IF_CONSTEXPR)
874// https://stackoverflow.com/questions/48979861/numerically-stable-method-for-solving-quadratic-equations/50065711
875namespace detail
876{
877#if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1)
878inline float fma_workaround(float x, float y, float z) { return ::fmaf(x, y, z); }
879inline double fma_workaround(double x, double y, double z) { return ::fma(x, y, z); }
880#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
881inline long double fma_workaround(long double x, long double y, long double z) { return ::fmal(x, y, z); }
882#endif
883#endif
884template<class T>
885inline T discriminant(T const& a, T const& b, T const& c)
886{
887 T w = 4 * a * c;
888#if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1)
889 T e = fma_workaround(-c, 4 * a, w);
890 T f = fma_workaround(b, b, -w);
891#else
892 T e = std::fma(-c, 4 * a, w);
893 T f = std::fma(b, b, -w);
894#endif
895 return f + e;
896}
897
898template<class T>
899std::pair<T, T> quadratic_roots_imp(T const& a, T const& b, T const& c)
900{
901#if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1)
902 using boost::math::copysign;
903#else
904 using std::copysign;
905#endif
906 using std::sqrt;
907 if constexpr (std::is_floating_point<T>::value)
908 {
909 T nan = std::numeric_limits<T>::quiet_NaN();
910 if (a == 0)
911 {
912 if (b == 0 && c != 0)
913 {
914 return std::pair<T, T>(nan, nan);
915 }
916 else if (b == 0 && c == 0)
917 {
918 return std::pair<T, T>(0, 0);
919 }
920 return std::pair<T, T>(-c / b, -c / b);
921 }
922 if (b == 0)
923 {
924 T x0_sq = -c / a;
925 if (x0_sq < 0) {
926 return std::pair<T, T>(nan, nan);
927 }
928 T x0 = sqrt(x0_sq);
929 return std::pair<T, T>(-x0, x0);
930 }
931 T discriminant = detail::discriminant(a, b, c);
932 // Is there a sane way to flush very small negative values to zero?
933 // If there is I don't know of it.
934 if (discriminant < 0)
935 {
936 return std::pair<T, T>(nan, nan);
937 }
938 T q = -(b + copysign(sqrt(discriminant), b)) / T(2);
939 T x0 = q / a;
940 T x1 = c / q;
941 if (x0 < x1)
942 {
943 return std::pair<T, T>(x0, x1);
944 }
945 return std::pair<T, T>(x1, x0);
946 }
947 else if constexpr (boost::math::tools::is_complex_type<T>::value)
948 {
949 typename T::value_type nan = std::numeric_limits<typename T::value_type>::quiet_NaN();
950 if (a.real() == 0 && a.imag() == 0)
951 {
952 using std::norm;
953 if (b.real() == 0 && b.imag() && norm(c) != 0)
954 {
955 return std::pair<T, T>({ nan, nan }, { nan, nan });
956 }
957 else if (b.real() == 0 && b.imag() && c.real() == 0 && c.imag() == 0)
958 {
959 return std::pair<T, T>({ 0,0 }, { 0,0 });
960 }
961 return std::pair<T, T>(-c / b, -c / b);
962 }
963 if (b.real() == 0 && b.imag() == 0)
964 {
965 T x0_sq = -c / a;
966 T x0 = sqrt(x0_sq);
967 return std::pair<T, T>(-x0, x0);
968 }
969 // There's no fma for complex types:
970 T discriminant = b * b - T(4) * a * c;
971 T q = -(b + sqrt(discriminant)) / T(2);
972 return std::pair<T, T>(q / a, c / q);
973 }
974 else // Most likely the type is a boost.multiprecision.
975 { //There is no fma for multiprecision, and in addition it doesn't seem to be useful, so revert to the naive computation.
976 T nan = std::numeric_limits<T>::quiet_NaN();
977 if (a == 0)
978 {
979 if (b == 0 && c != 0)
980 {
981 return std::pair<T, T>(nan, nan);
982 }
983 else if (b == 0 && c == 0)
984 {
985 return std::pair<T, T>(0, 0);
986 }
987 return std::pair<T, T>(-c / b, -c / b);
988 }
989 if (b == 0)
990 {
991 T x0_sq = -c / a;
992 if (x0_sq < 0) {
993 return std::pair<T, T>(nan, nan);
994 }
995 T x0 = sqrt(x0_sq);
996 return std::pair<T, T>(-x0, x0);
997 }
998 T discriminant = b * b - 4 * a * c;
999 if (discriminant < 0)
1000 {
1001 return std::pair<T, T>(nan, nan);
1002 }
1003 T q = -(b + copysign(sqrt(discriminant), b)) / T(2);
1004 T x0 = q / a;
1005 T x1 = c / q;
1006 if (x0 < x1)
1007 {
1008 return std::pair<T, T>(x0, x1);
1009 }
1010 return std::pair<T, T>(x1, x0);
1011 }
1012}
1013} // namespace detail
1014
1015template<class T1, class T2 = T1, class T3 = T1>
1016inline std::pair<typename tools::promote_args<T1, T2, T3>::type, typename tools::promote_args<T1, T2, T3>::type> quadratic_roots(T1 const& a, T2 const& b, T3 const& c)
1017{
1018 typedef typename tools::promote_args<T1, T2, T3>::type value_type;
1019 return detail::quadratic_roots_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(c));
1020}
1021
1022#endif
1023
1024} // namespace tools
1025} // namespace math
1026} // namespace boost
1027
1028#endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP
1029

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