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 | |
39 | namespace boost { |
40 | namespace math { |
41 | namespace tools { |
42 | |
43 | namespace detail { |
44 | |
45 | namespace dummy { |
46 | |
47 | template<int n, class T> |
48 | typename T::value_type get(const T&) BOOST_MATH_NOEXCEPT(T); |
49 | } |
50 | |
51 | template <class Tuple, class T> |
52 | void 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 | } |
59 | template <class Tuple, class T> |
60 | void 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 | |
69 | template <class Tuple, class T> |
70 | inline 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 | |
77 | template <class T, class U, class V> |
78 | inline 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 | } |
83 | template <class T, class U, class V> |
84 | inline void unpack_0(const std::pair<T, U>& p, V& a) BOOST_MATH_NOEXCEPT(T) |
85 | { |
86 | a = p.first; |
87 | } |
88 | |
89 | template <class F, class T> |
90 | void 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 | |
142 | template <class F, class T, class Tol, class Policy> |
143 | std::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 | |
221 | template <class F, class T, class Tol> |
222 | inline 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 | |
227 | template <class F, class T, class Tol> |
228 | inline 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 | |
235 | template <class F, class T> |
236 | T 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 | |
365 | template <class F, class T> |
366 | inline 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 | |
372 | namespace 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 | |
707 | template <class F, class T> |
708 | T 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 | |
713 | template <class F, class T> |
714 | inline 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 | |
720 | namespace 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 | |
745 | template <class F, class T> |
746 | T 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 | |
751 | template <class F, class T> |
752 | inline 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 | // |
760 | template <class F, class T> |
761 | T 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 | |
766 | template <class F, class T> |
767 | inline 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 | */ |
780 | template<class Complex, class F> |
781 | Complex 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 |
862 | namespace detail |
863 | { |
864 | #if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1) |
865 | inline float fma_workaround(float x, float y, float z) { return ::fmaf(x, y, z); } |
866 | inline double fma_workaround(double x, double y, double z) { return ::fma(x, y, z); } |
867 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
868 | inline long double fma_workaround(long double x, long double y, long double z) { return ::fmal(x, y, z); } |
869 | #endif |
870 | #endif |
871 | template<class T> |
872 | inline 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 | |
885 | template<class T> |
886 | std::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 | |
1002 | template<class T1, class T2 = T1, class T3 = T1> |
1003 | inline 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 | |