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

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