]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | |
92f5a8d4 | 12 | #include <boost/math/tools/complex.hpp> // test for multiprecision types. |
7c673cae | 13 | |
92f5a8d4 | 14 | #include <iostream> |
7c673cae FG |
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> | |
f67539c2 | 23 | #include <boost/math/tools/cxx03_warn.hpp> |
7c673cae FG |
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> | |
92f5a8d4 | 35 | #include <boost/math/special_functions/next.hpp> |
7c673cae FG |
36 | #include <boost/math/tools/toms748_solve.hpp> |
37 | #include <boost/math/policies/error_handling.hpp> | |
38 | ||
92f5a8d4 TL |
39 | namespace boost { |
40 | namespace math { | |
41 | namespace tools { | |
7c673cae | 42 | |
92f5a8d4 | 43 | namespace detail { |
7c673cae | 44 | |
92f5a8d4 | 45 | namespace dummy { |
7c673cae FG |
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: | |
92f5a8d4 | 74 | val = get<0>(t); |
7c673cae FG |
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, | |
92f5a8d4 TL |
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>()))) | |
7c673cae | 98 | { |
92f5a8d4 | 99 | if (last_f0 == 0) |
7c673cae FG |
100 | { |
101 | // this must be the first iteration, pretend that we had a | |
102 | // previous one at either min or max: | |
92f5a8d4 | 103 | if (result == min) |
7c673cae FG |
104 | { |
105 | guess = max; | |
106 | } | |
107 | else | |
108 | { | |
109 | guess = min; | |
110 | } | |
111 | unpack_0(f(guess), last_f0); | |
112 | delta = guess - result; | |
113 | } | |
92f5a8d4 | 114 | if (sign(last_f0) * sign(f0) < 0) |
7c673cae FG |
115 | { |
116 | // we've crossed over so move in opposite direction to last step: | |
92f5a8d4 | 117 | if (delta < 0) |
7c673cae FG |
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: | |
92f5a8d4 | 129 | if (delta < 0) |
7c673cae FG |
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> | |
92f5a8d4 | 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>()))) |
7c673cae FG |
144 | { |
145 | T fmin = f(min); | |
146 | T fmax = f(max); | |
92f5a8d4 | 147 | if (fmin == 0) |
7c673cae FG |
148 | { |
149 | max_iter = 2; | |
150 | return std::make_pair(min, min); | |
151 | } | |
92f5a8d4 | 152 | if (fmax == 0) |
7c673cae FG |
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%>"; | |
92f5a8d4 | 162 | if (min >= max) |
7c673cae FG |
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 | } | |
92f5a8d4 | 167 | if (fmin * fmax >= 0) |
7c673cae FG |
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; | |
92f5a8d4 | 177 | if (count < 3) |
7c673cae FG |
178 | count = 0; |
179 | else | |
180 | count -= 3; | |
181 | ||
92f5a8d4 | 182 | while (count && (0 == tol(min, max))) |
7c673cae FG |
183 | { |
184 | T mid = (min + max) / 2; | |
185 | T fmid = f(mid); | |
92f5a8d4 | 186 | if ((mid == max) || (mid == min)) |
7c673cae | 187 | break; |
92f5a8d4 | 188 | if (fmid == 0) |
7c673cae FG |
189 | { |
190 | min = max = mid; | |
191 | break; | |
192 | } | |
92f5a8d4 | 193 | else if (sign(fmid) * sign(fmin) < 0) |
7c673cae FG |
194 | { |
195 | max = mid; | |
7c673cae FG |
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; | |
92f5a8d4 | 211 | if (max_iter > max_count) |
7c673cae FG |
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> | |
92f5a8d4 | 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>()))) |
7c673cae FG |
223 | { |
224 | return bisect(f, min, max, tol, max_iter, policies::policy<>()); | |
225 | } | |
226 | ||
227 | template <class F, class T, class Tol> | |
92f5a8d4 | 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>()))) |
7c673cae FG |
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> | |
92f5a8d4 | 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>()))) |
7c673cae FG |
237 | { |
238 | BOOST_MATH_STD_USING | |
239 | ||
92f5a8d4 TL |
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 | ||
7c673cae FG |
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 | ||
92f5a8d4 TL |
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 | ||
7c673cae FG |
267 | boost::uintmax_t count(max_iter); |
268 | ||
92f5a8d4 TL |
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 { | |
7c673cae FG |
275 | last_f0 = f0; |
276 | delta2 = delta1; | |
277 | delta1 = delta; | |
278 | detail::unpack_tuple(f(result), f0, f1); | |
279 | --count; | |
92f5a8d4 | 280 | if (0 == f0) |
7c673cae | 281 | break; |
92f5a8d4 | 282 | if (f1 == 0) |
7c673cae FG |
283 | { |
284 | // Oops zero derivative!!! | |
285 | #ifdef BOOST_MATH_INSTRUMENT | |
92f5a8d4 | 286 | std::cout << "Newton iteration, zero derivative found!" << std::endl; |
7c673cae FG |
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 | |
92f5a8d4 | 295 | std::cout << "Newton iteration " << max_iter - count << ", delta = " << delta << std::endl; |
7c673cae | 296 | #endif |
92f5a8d4 | 297 | if (fabs(delta * 2) > fabs(delta2)) |
7c673cae | 298 | { |
92f5a8d4 TL |
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; | |
7c673cae FG |
311 | } |
312 | guess = result; | |
313 | result -= delta; | |
92f5a8d4 | 314 | if (result <= min) |
7c673cae FG |
315 | { |
316 | delta = 0.5F * (guess - min); | |
317 | result = guess - delta; | |
92f5a8d4 | 318 | if ((result == min) || (result == max)) |
7c673cae FG |
319 | break; |
320 | } | |
92f5a8d4 | 321 | else if (result >= max) |
7c673cae FG |
322 | { |
323 | delta = 0.5F * (guess - max); | |
324 | result = guess - delta; | |
92f5a8d4 | 325 | if ((result == min) || (result == max)) |
7c673cae FG |
326 | break; |
327 | } | |
92f5a8d4 TL |
328 | // Update brackets: |
329 | if (delta > 0) | |
330 | { | |
7c673cae | 331 | max = guess; |
92f5a8d4 TL |
332 | max_range_f = f0; |
333 | } | |
7c673cae | 334 | else |
92f5a8d4 | 335 | { |
7c673cae | 336 | min = guess; |
92f5a8d4 TL |
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 | } | |
7c673cae FG |
346 | }while(count && (fabs(result * factor) < fabs(delta))); |
347 | ||
348 | max_iter -= count; | |
349 | ||
350 | #ifdef BOOST_MATH_INSTRUMENT | |
92f5a8d4 | 351 | std::cout << "Newton Raphson final iteration count = " << max_iter << std::endl; |
7c673cae FG |
352 | |
353 | static boost::uintmax_t max_count = 0; | |
92f5a8d4 | 354 | if (max_iter > max_count) |
7c673cae FG |
355 | { |
356 | max_count = max_iter; | |
92f5a8d4 TL |
357 | // std::cout << "Maximum iterations: " << max_iter << std::endl; |
358 | // Puzzled what this tells us, so commented out for now? | |
7c673cae FG |
359 | } |
360 | #endif | |
361 | ||
362 | return result; | |
363 | } | |
364 | ||
365 | template <class F, class T> | |
92f5a8d4 | 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>()))) |
7c673cae FG |
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 | ||
92f5a8d4 | 372 | namespace detail { |
7c673cae FG |
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 | ||
92f5a8d4 | 387 | if ((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>())) |
7c673cae FG |
388 | { |
389 | // possible overflow, use Newton step: | |
390 | delta = f0 / f1; | |
391 | } | |
392 | else | |
393 | delta = denom / num; | |
394 | return delta; | |
395 | } | |
396 | }; | |
397 | ||
92f5a8d4 TL |
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 | ||
7c673cae | 512 | template <class Stepper, class F, class T> |
92f5a8d4 | 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>()))) |
7c673cae FG |
514 | { |
515 | BOOST_MATH_STD_USING | |
516 | ||
92f5a8d4 TL |
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; | |
7c673cae FG |
528 | T result = guess; |
529 | ||
92f5a8d4 | 530 | T factor = ldexp(static_cast<T>(1.0), 1 - digits); |
f67539c2 | 531 | T delta = (std::max)(T(10000000 * guess), T(10000000)); // arbitrarily large delta |
7c673cae FG |
532 | T last_f0 = 0; |
533 | T delta1 = delta; | |
534 | T delta2 = delta; | |
7c673cae FG |
535 | bool out_of_bounds_sentry = false; |
536 | ||
92f5a8d4 | 537 | #ifdef BOOST_MATH_INSTRUMENT |
7c673cae | 538 | std::cout << "Second order root iteration, limit = " << factor << std::endl; |
92f5a8d4 TL |
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; | |
7c673cae FG |
553 | |
554 | boost::uintmax_t count(max_iter); | |
555 | ||
92f5a8d4 | 556 | do { |
7c673cae FG |
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 | ||
92f5a8d4 | 567 | if (0 == f0) |
7c673cae | 568 | break; |
92f5a8d4 | 569 | if (f1 == 0) |
7c673cae FG |
570 | { |
571 | // Oops zero derivative!!! | |
92f5a8d4 TL |
572 | #ifdef BOOST_MATH_INSTRUMENT |
573 | std::cout << "Second order root iteration, zero derivative found!" << std::endl; | |
574 | #endif | |
7c673cae FG |
575 | detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max); |
576 | } | |
577 | else | |
578 | { | |
92f5a8d4 | 579 | if (f2 != 0) |
7c673cae FG |
580 | { |
581 | delta = Stepper::step(result, f0, f1, f2); | |
92f5a8d4 | 582 | if (delta * f1 / f0 < 0) |
7c673cae FG |
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; | |
92f5a8d4 | 594 | if (fabs(delta) > 2 * fabs(guess)) |
7c673cae FG |
595 | delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess); |
596 | } | |
597 | } | |
598 | else | |
599 | delta = f0 / f1; | |
600 | } | |
92f5a8d4 | 601 | #ifdef BOOST_MATH_INSTRUMENT |
7c673cae | 602 | std::cout << "Second order root iteration, delta = " << delta << std::endl; |
92f5a8d4 | 603 | #endif |
7c673cae | 604 | T convergence = fabs(delta / delta2); |
92f5a8d4 | 605 | if ((convergence > 0.8) && (convergence < 2)) |
7c673cae | 606 | { |
92f5a8d4 | 607 | // last two steps haven't converged. |
7c673cae | 608 | delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; |
92f5a8d4 TL |
609 | if ((result != 0) && (fabs(delta) > result)) |
610 | delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps! | |
7c673cae FG |
611 | // reset delta2 so that this branch will *not* be taken on the |
612 | // next iteration: | |
613 | delta2 = delta * 3; | |
92f5a8d4 | 614 | delta1 = delta * 3; |
7c673cae FG |
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: | |
92f5a8d4 | 622 | if (result < min) |
7c673cae | 623 | { |
92f5a8d4 TL |
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) | |
7c673cae | 629 | diff = 1 / diff; |
92f5a8d4 | 630 | if (!out_of_bounds_sentry && (diff > 0) && (diff < 3)) |
7c673cae FG |
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 | { | |
92f5a8d4 TL |
640 | if (fabs(float_distance(min, max)) < 2) |
641 | { | |
642 | result = guess = (min + max) / 2; | |
7c673cae | 643 | break; |
92f5a8d4 TL |
644 | } |
645 | delta = bracket_root_towards_min(f, guess, f0, min, max, count); | |
646 | result = guess - delta; | |
647 | guess = min; | |
648 | continue; | |
7c673cae FG |
649 | } |
650 | } | |
92f5a8d4 | 651 | else if (result > max) |
7c673cae FG |
652 | { |
653 | T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max); | |
92f5a8d4 | 654 | if (fabs(diff) < 1) |
7c673cae | 655 | diff = 1 / diff; |
92f5a8d4 | 656 | if (!out_of_bounds_sentry && (diff > 0) && (diff < 3)) |
7c673cae FG |
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 | { | |
92f5a8d4 TL |
666 | if (fabs(float_distance(min, max)) < 2) |
667 | { | |
668 | result = guess = (min + max) / 2; | |
7c673cae | 669 | break; |
92f5a8d4 TL |
670 | } |
671 | delta = bracket_root_towards_max(f, guess, f0, min, max, count); | |
672 | result = guess - delta; | |
673 | guess = min; | |
674 | continue; | |
7c673cae FG |
675 | } |
676 | } | |
677 | // update brackets: | |
92f5a8d4 TL |
678 | if (delta > 0) |
679 | { | |
7c673cae | 680 | max = guess; |
92f5a8d4 TL |
681 | max_range_f = f0; |
682 | } | |
7c673cae | 683 | else |
92f5a8d4 | 684 | { |
7c673cae | 685 | min = guess; |
92f5a8d4 TL |
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 | } | |
7c673cae FG |
695 | } while(count && (fabs(result * factor) < fabs(delta))); |
696 | ||
697 | max_iter -= count; | |
698 | ||
92f5a8d4 TL |
699 | #ifdef BOOST_MATH_INSTRUMENT |
700 | std::cout << "Second order root finder, final iteration count = " << max_iter << std::endl; | |
701 | #endif | |
7c673cae FG |
702 | |
703 | return result; | |
704 | } | |
92f5a8d4 | 705 | } // T second_order_root_finder |
7c673cae FG |
706 | |
707 | template <class F, class T> | |
92f5a8d4 | 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>()))) |
7c673cae FG |
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> | |
92f5a8d4 | 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>()))) |
7c673cae FG |
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 | ||
92f5a8d4 | 720 | namespace detail { |
7c673cae FG |
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 | { | |
92f5a8d4 | 727 | using std::fabs; |
7c673cae FG |
728 | T ratio = f0 / f1; |
729 | T delta; | |
92f5a8d4 | 730 | if ((x != 0) && (fabs(ratio / x) < 0.1)) |
7c673cae FG |
731 | { |
732 | delta = ratio + (f2 / (2 * f1)) * ratio * ratio; | |
733 | // check second derivative doesn't over compensate: | |
92f5a8d4 | 734 | if (delta * ratio < 0) |
7c673cae FG |
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> | |
92f5a8d4 | 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>()))) |
7c673cae FG |
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> | |
92f5a8d4 | 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>()))) |
7c673cae FG |
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 | // | |
f67539c2 | 758 | // These two are the old spelling of this function, retained for backwards compatibility just in case: |
7c673cae FG |
759 | // |
760 | template <class F, class T> | |
92f5a8d4 | 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>()))) |
7c673cae FG |
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> | |
92f5a8d4 | 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>()))) |
7c673cae FG |
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 | ||
92f5a8d4 TL |
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) | |
f67539c2 TL |
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); } | |
92f5a8d4 | 867 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
f67539c2 | 868 | inline long double fma_workaround(long double x, long double y, long double z) { return ::fmal(x, y, z); } |
92f5a8d4 TL |
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 | { | |
f67539c2 TL |
888 | #if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1) |
889 | using boost::math::copysign; | |
890 | #else | |
92f5a8d4 | 891 | using std::copysign; |
f67539c2 | 892 | #endif |
92f5a8d4 TL |
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 | |
7c673cae FG |
1010 | |
1011 | } // namespace tools | |
1012 | } // namespace math | |
1013 | } // namespace boost | |
1014 | ||
1015 | #endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP |