RandLib 2.0.0
Probability distributions, fitting, and heavy-tail tools for C++20
Loading...
Searching...
No Matches
roots.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <randlib/math/constants.hpp>
4
5#include <algorithm>
6#include <cmath>
7#include <concepts>
8#include <functional>
9#include <limits>
10#include <stdexcept>
11#include <tuple>
12#include <type_traits>
13#include <utility>
14
15namespace randlib::math
16{
17
18namespace detail
19{
20
22template <std::floating_point Real>
23bool bracketResolved(Real a, Real b)
24{
25 const Real scale = std::min(std::fabs(a), std::fabs(b));
26 return b - a <= 4 * std::numeric_limits<Real>::epsilon() * scale ||
27 b - a <= 4 * std::numeric_limits<Real>::denorm_min();
28}
29
30template <std::floating_point Real>
31Real safeDivide(Real numerator, Real denominator, Real fallback)
32{
33 if (std::fabs(denominator) < 1 && std::fabs(denominator * std::numeric_limits<Real>::max()) <= std::fabs(numerator))
34 return fallback;
35 return numerator / denominator;
36}
37
38template <std::floating_point Real>
39Real secantStep(Real a, Real b, Real fa, Real fb)
40{
41 const Real tolerance = 5 * std::numeric_limits<Real>::epsilon();
42 Real c = a - (fa / (fb - fa)) * (b - a);
43 if (c <= a + std::fabs(a) * tolerance || c >= b - std::fabs(b) * tolerance)
44 return a + (b - a) / 2;
45 return c;
46}
47
48template <std::floating_point Real>
49Real newtonQuadraticStep(Real a, Real b, Real d, Real fa, Real fb, Real fd, int steps)
50{
51 Real B = safeDivide(fb - fa, b - a, std::numeric_limits<Real>::max());
52 Real A = safeDivide(fd - fb, d - b, std::numeric_limits<Real>::max());
53 A = safeDivide(A - B, d - a, Real(0));
54 if (A == 0)
55 return secantStep(a, b, fa, fb);
56 Real c = (std::signbit(A) == std::signbit(fa)) ? a : b;
57 for (int i = 0; i < steps; ++i)
58 c -= safeDivide(fa + (B + A * (c - b)) * (c - a), B + A * (2 * c - a - b), 1 + c - a);
59 if (!(c > a && c < b))
60 c = secantStep(a, b, fa, fb);
61 return c;
62}
63
64template <std::floating_point Real>
65Real inverseCubicStep(Real a, Real b, Real d, Real e, Real fa, Real fb, Real fd, Real fe)
66{
67 Real q11 = (d - e) * fd / (fe - fd);
68 Real q21 = (b - d) * fb / (fd - fb);
69 Real q31 = (a - b) * fa / (fb - fa);
70 Real d21 = (b - d) * fd / (fd - fb);
71 Real d31 = (a - b) * fb / (fb - fa);
72 Real q22 = (d21 - q11) * fb / (fe - fb);
73 Real q32 = (d31 - q21) * fa / (fd - fa);
74 Real d32 = (d31 - q21) * fd / (fd - fa);
75 Real q33 = (d32 - q22) * fa / (fe - fa);
76 Real c = q31 + q32 + q33 + a;
77 if (!(c > a && c < b))
78 c = newtonQuadraticStep(a, b, d, fa, fb, fd, 3);
79 return c;
80}
81
82} // namespace detail
83
97template <std::floating_point Real, class F>
98std::pair<Real, Real> toms748(const F &f, Real a, Real b, Real fa, Real fb, int maxIterations = 200)
99{
100 if (a > b) {
101 std::swap(a, b);
102 std::swap(fa, fb);
103 }
104 if (fa == 0)
105 return {a, a};
106 if (fb == 0)
107 return {b, b};
108 if (std::signbit(fa) == std::signbit(fb))
109 throw std::invalid_argument("toms748: the root is not bracketed");
110
111 Real d = 0, fd = 0, e = 0, fe = 0;
112 int count = maxIterations;
114 auto bracket = [&f, &a, &b, &fa, &fb, &d, &fd](Real c) {
115 const Real tolerance = 2 * std::numeric_limits<Real>::epsilon();
116 if (b - a < 2 * tolerance * std::fabs(a))
117 c = a + (b - a) / 2;
118 else if (c <= a + std::fabs(a) * tolerance)
119 c = a + std::fabs(a) * tolerance;
120 else if (c >= b - std::fabs(b) * tolerance)
121 c = b - std::fabs(b) * tolerance;
122 if (!(c > a && c < b))
123 c = a + (b - a) / 2;
124 Real fc = f(c);
125 if (fc == 0) {
126 a = b = c;
127 fa = fb = 0;
128 return;
129 }
130 if (std::signbit(fa) != std::signbit(fc)) {
131 d = b;
132 fd = fb;
133 b = c;
134 fb = fc;
135 } else {
136 d = a;
137 fd = fa;
138 a = c;
139 fa = fc;
140 }
141 };
142 auto done = [&] { return count <= 0 || fa == 0 || detail::bracketResolved(a, b); };
143
144 bracket(detail::secantStep(a, b, fa, fb));
145 --count;
146 if (!done()) {
147 e = d;
148 fe = fd;
149 bracket(detail::newtonQuadraticStep(a, b, d, fa, fb, fd, 2));
150 --count;
151 }
152 const Real minDifference = 32 * std::numeric_limits<Real>::min();
153 auto distinct = [&] {
154 return std::fabs(fa - fb) >= minDifference && std::fabs(fa - fd) >= minDifference &&
155 std::fabs(fa - fe) >= minDifference && std::fabs(fb - fd) >= minDifference &&
156 std::fabs(fb - fe) >= minDifference && std::fabs(fd - fe) >= minDifference;
157 };
158 while (!done()) {
159 const Real a0 = a, b0 = b;
160 Real c = distinct() ? detail::inverseCubicStep(a, b, d, e, fa, fb, fd, fe)
161 : detail::newtonQuadraticStep(a, b, d, fa, fb, fd, 2);
162 e = d;
163 fe = fd;
164 bracket(c);
165 --count;
166 if (done())
167 break;
168 c = distinct() ? detail::inverseCubicStep(a, b, d, e, fa, fb, fd, fe)
169 : detail::newtonQuadraticStep(a, b, d, fa, fb, fd, 3);
170 bracket(c);
171 --count;
172 if (done())
173 break;
175 const bool aIsCloser = std::fabs(fa) < std::fabs(fb);
176 const Real u = aIsCloser ? a : b, fu = aIsCloser ? fa : fb;
177 c = u - 2 * (fu / (fb - fa)) * (b - a);
178 if (std::fabs(c - u) > (b - a) / 2)
179 c = a + (b - a) / 2;
180 e = d;
181 fe = fd;
182 bracket(c);
183 --count;
184 if (done())
185 break;
187 if (b - a < (b0 - a0) / 2)
188 continue;
189 e = d;
190 fe = fd;
191 bracket(a + (b - a) / 2);
192 --count;
193 }
194 return {a, b};
195}
196
204template <std::floating_point Real = double, class F>
205Real findRoot(const F &f, std::type_identity_t<Real> a, std::type_identity_t<Real> b, int maxIterations = 200)
206{
207 auto [lo, hi] = toms748<Real>(f, a, b, f(a), f(b), maxIterations);
208 return lo + (hi - lo) / 2;
209}
210
220template <typename RealType>
221bool findRootNewtonFirstOrder(const std::function<DoublePair(RealType)> &funPtr, RealType &root,
222 long double funTol = 1e-10, long double stepTol = 1e-6)
223{
225 funTol = funTol > math::detail::minimumTolerance ? funTol : math::detail::minimumTolerance;
226 stepTol = stepTol > math::detail::minimumTolerance ? stepTol : math::detail::minimumTolerance;
227 static constexpr int MAX_ITER = 1e5;
228 static constexpr double MAX_STEP = 10;
229 int iter = 0;
230 double step = stepTol + 1;
231 DoublePair y = funPtr(root);
232 double fun = y.first;
233 double grad = y.second;
234 if (!std::isfinite(root) || !std::isfinite(fun) || !std::isfinite(grad))
235 return false;
236 if (std::fabs(fun) < math::detail::minimumTolerance)
237 return true;
238 do {
239 double alpha = 1.0;
240 double oldRoot = root;
241 double oldFun = fun;
242 double newtonStep = fun / grad;
243 if (!std::isfinite(newtonStep))
244 return false;
245 step = std::min(MAX_STEP, std::max(-MAX_STEP, newtonStep));
246 bool accepted = false;
247 for (int backtrack = 0; backtrack < 64; ++backtrack) {
248 root = oldRoot - alpha * step;
249 if (!std::isfinite(root)) {
250 alpha *= 0.5;
251 continue;
252 }
253 y = funPtr(root);
254 fun = y.first;
255 grad = y.second;
256 if (std::fabs(fun) < math::detail::minimumTolerance)
257 return true;
258 if (std::isfinite(fun) && std::isfinite(grad) && std::fabs(grad) > math::detail::minimumTolerance &&
259 std::fabs(fun) <= std::fabs(oldFun)) {
260 accepted = true;
261 break;
262 }
263 alpha *= 0.5;
264 }
265 if (!accepted) {
266 root = oldRoot;
267 return false;
268 }
270 double diffX = std::fabs(root - oldRoot);
271 double relDiffX = std::fabs(diffX / oldRoot);
272 if (std::min(diffX, relDiffX) < stepTol) {
273 double diffY = fun - oldFun;
274 double relDiffY = std::fabs(diffY / oldFun);
275 if (std::min(std::fabs(fun), relDiffY) < funTol)
276 return true;
277 }
278 } while (++iter < MAX_ITER);
279 return false;
280}
281
291template <typename RealType>
292bool findRootNewtonSecondOrder(const std::function<DoubleTriplet(RealType)> &funPtr, RealType &root,
293 long double funTol = 1e-10, long double stepTol = 1e-6)
294{
296 funTol = funTol > math::detail::minimumTolerance ? funTol : math::detail::minimumTolerance;
297 stepTol = stepTol > math::detail::minimumTolerance ? stepTol : math::detail::minimumTolerance;
298 static constexpr int MAX_ITER = 1e5;
299 static constexpr double MAX_STEP = 10;
300 int iter = 0;
301 double step = stepTol + 1;
302 auto [f, fx, fxx] = funPtr(root);
303 if (std::fabs(f) < math::detail::minimumTolerance)
304 return true;
305 do {
306 double alpha = 1.0;
307 double oldRoot = root;
308 double oldFun = f;
309 double numerator = 2 * f * fx;
310 double denominator = 2 * fx * fx - f * fxx;
311 step = std::min(MAX_STEP, std::max(-MAX_STEP, numerator / denominator));
312 do {
313 root = oldRoot - alpha * step;
314 std::tie(f, fx, fxx) = funPtr(root);
315 if (std::fabs(f) < math::detail::minimumTolerance)
316 return true;
317 alpha *= 0.5;
318 } while ((std::fabs(fx) <= math::detail::minimumTolerance || std::fabs(oldFun) < std::fabs(f)) && alpha > 0);
320 double diffX = std::fabs(root - oldRoot);
321 double relDiffX = std::fabs(diffX / oldRoot);
322 if (std::min(diffX, relDiffX) < stepTol) {
323 double diffY = f - oldFun;
324 double relDiffY = std::fabs(diffY / oldFun);
325 if (std::min(std::fabs(f), relDiffY) < funTol)
326 return true;
327 }
328 } while (++iter < MAX_ITER);
329 return false;
330}
331
343bool findRootNewtonFirstOrder2d(const std::function<DoublePair(DoublePair)> &funPtr,
344 const std::function<std::tuple<DoublePair, DoublePair>(DoublePair)> &gradPtr,
345 DoublePair &root, long double funTol = 1e-10, long double stepTol = 1e-6);
346
347inline bool findRootNewtonFirstOrder2d(const std::function<DoublePair(DoublePair)> &funPtr,
348 const std::function<std::tuple<DoublePair, DoublePair>(DoublePair)> &gradPtr,
349 DoublePair &root, long double funTol, long double stepTol)
350{
352 funTol = funTol > math::detail::minimumTolerance ? funTol : math::detail::minimumTolerance;
353 stepTol = stepTol > math::detail::minimumTolerance ? stepTol : math::detail::minimumTolerance;
354 static constexpr int MAX_ITER = 1e5;
355 static constexpr double MAX_STEP = 10;
356 int iter = 0;
357 double step1;
358 double step2;
359 DoublePair fun = funPtr(root);
360 double fun1 = fun.first;
361 double fun2 = fun.second;
362 double error = std::max(std::fabs(fun1), std::fabs(fun2));
363 if (error < math::detail::minimumTolerance)
364 return true;
365
366 auto [grad1, grad2] = gradPtr(root);
367 do {
368 double alpha = 1.0;
369 DoublePair oldRoot = root;
370 DoublePair oldFun = fun;
371 double oldError = error;
372 double det = grad1.first * grad2.second - grad1.second * grad2.first;
373 step1 = std::min(MAX_STEP, std::max(-MAX_STEP, (grad2.second * fun1 - grad1.second * fun2) / det));
374 step2 = std::min(MAX_STEP, std::max(-MAX_STEP, (-grad2.first * fun1 + grad1.first * fun2) / det));
375 do {
376 root.first = oldRoot.first - alpha * step1;
377 root.second = oldRoot.second - alpha * step2;
378 fun = funPtr(root);
379 fun1 = fun.first;
380 fun2 = fun.second;
381 error = std::max(std::fabs(fun1), std::fabs(fun2));
382 if (error < math::detail::minimumTolerance)
383 return true;
384 std::tie(grad1, grad2) = gradPtr(root);
385 det = grad1.first * grad2.second - grad1.second * grad2.first;
386 alpha *= 0.5;
387 } while ((std::fabs(det) <= math::detail::minimumTolerance || oldError < error) && alpha > 0);
388
390 double diffX1 = std::fabs(root.first - oldRoot.first);
391 double diffX2 = std::fabs(root.second - oldRoot.second);
392 double diffX = std::max(diffX1, diffX2);
393 double relDiffX1 = std::fabs(diffX1 / oldRoot.first);
394 double relDiffX2 = std::fabs(diffX2 / oldRoot.second);
395 double relDiffX = std::max(relDiffX1, relDiffX2);
396 if (std::min(diffX, relDiffX) < stepTol) {
397 double diffY1 = fun1 - oldFun.first;
398 double diffY2 = fun2 - oldFun.second;
399 double relDiffY1 = std::fabs(diffY1 / oldFun.first);
400 double relDiffY2 = std::fabs(diffY2 / oldFun.second);
401 double relDiffY = std::max(relDiffY1, relDiffY2);
402 if (std::min(std::max(std::fabs(fun1), std::fabs(fun2)), relDiffY) < funTol)
403 return true;
404 }
405 } while (++iter < MAX_ITER);
406 return false;
407}
408
409} // namespace randlib::math
randlib::math::detail::safeDivide
Real safeDivide(Real numerator, Real denominator, Real fallback)
Definition roots.hpp:31
randlib::math::detail::newtonQuadraticStep
Real newtonQuadraticStep(Real a, Real b, Real d, Real fa, Real fb, Real fd, int steps)
Definition roots.hpp:49
randlib::math::detail::secantStep
Real secantStep(Real a, Real b, Real fa, Real fb)
Definition roots.hpp:39
randlib::math::detail::minimumTolerance
constexpr long double minimumTolerance
Definition constants.hpp:44
randlib::math::detail::bracketResolved
bool bracketResolved(Real a, Real b)
|b - a| no longer resolvable at the relative precision of Real
Definition roots.hpp:23
randlib::math::detail::inverseCubicStep
Real inverseCubicStep(Real a, Real b, Real d, Real e, Real fa, Real fb, Real fd, Real fe)
Definition roots.hpp:65
randlib::math
Beta-related functions.
Definition math/beta.hpp:14
randlib::math::findRootNewtonFirstOrder
bool findRootNewtonFirstOrder(const std::function< DoublePair(RealType)> &funPtr, RealType &root, long double funTol=1e-10, long double stepTol=1e-6)
Newton's root-finding procedure, using first derivative.
Definition roots.hpp:221
randlib::math::findRootNewtonSecondOrder
bool findRootNewtonSecondOrder(const std::function< DoubleTriplet(RealType)> &funPtr, RealType &root, long double funTol=1e-10, long double stepTol=1e-6)
Newton's root-finding procedure, using first and second derivatives.
Definition roots.hpp:292
randlib::math::findRoot
Real findRoot(const F &f, std::type_identity_t< Real > a, std::type_identity_t< Real > b, int maxIterations=200)
Definition roots.hpp:205
randlib::math::findRootNewtonFirstOrder2d
bool findRootNewtonFirstOrder2d(const std::function< DoublePair(DoublePair)> &funPtr, const std::function< std::tuple< DoublePair, DoublePair >(DoublePair)> &gradPtr, DoublePair &root, long double funTol=1e-10, long double stepTol=1e-6)
Newton's root-finding procedure for 2 functions of 2 parameters, using Jacobian matrix.
Definition roots.hpp:347
randlib::math::toms748
std::pair< Real, Real > toms748(const F &f, Real a, Real b, Real fa, Real fb, int maxIterations=200)
Bracketing root finder of Alefeld, Potra and Shi (ACM TOMS algorithm 748): inverse cubic and Newton-q...
Definition roots.hpp:98
randlib::DoublePair
Pair< double > DoublePair
Definition fwd.hpp:39
randlib::DoubleTriplet
Triplet< double > DoubleTriplet
Definition fwd.hpp:43