Theoretica
A C++ numerical and automatic mathematical library
roots.h
Go to the documentation of this file.
1 
5 
6 #ifndef THEORETICA_ROOTS_H
7 #define THEORETICA_ROOTS_H
8 
9 #include "../core/function.h"
10 #include "../calculus/deriv.h"
11 #include "../autodiff/dual.h"
12 #include "../autodiff/dual2.h"
13 #include "../algebra/vec.h"
14 #include "../complex/complex.h"
15 
16 
17 namespace theoretica {
18 
19 
26  template<typename RealFunction>
27  inline std::vector<vec2> find_root_intervals(
28  RealFunction f, real a, real b, unsigned int steps = 10) {
29 
30  std::vector<vec2> res;
31  const real dx = (b - a) / (real) steps;
32 
33  for (unsigned int i = 0; i < steps; ++i) {
34 
35  const real x1 = a + i * dx;
36  const real x2 = a + (i + 1) * dx;
37 
38  if(f(x1) * f(x2) <= 0)
39  res.push_back({x1, x2});
40  }
41 
42  return res;
43  }
44 
45 
55  template<typename RealFunction>
56  inline real root_bisection(
57  RealFunction f, real a, real b, real tolerance = OPTIMIZATION_TOL) {
58 
59  if(f(a) * f(b) >= 0) {
60  TH_MATH_ERROR("root_bisection", f(a) * f(b), INVALID_ARGUMENT);
61  return nan();
62  }
63 
64  real x_avg = a;
65 
66  real x_min = a;
67  real x_max = b;
68 
69  unsigned int iter = 0;
70 
71  while((x_max - x_min) > tolerance && iter <= OPTIMIZATION_BISECTION_ITER) {
72 
73  x_avg = (x_max + x_min) / 2.0;
74 
75  if(f(x_avg) * f(x_min) > 0)
76  x_min = x_avg;
77  else
78  x_max = x_avg;
79 
80  ++iter;
81  }
82 
83  if(iter > OPTIMIZATION_BISECTION_ITER) {
84  TH_MATH_ERROR("root_bisection", x_avg, NO_ALGO_CONVERGENCE);
85  return nan();
86  }
87 
88  return x_avg;
89  }
90 
91 
99  template<typename RealFunction>
100  inline real root_newton(RealFunction f, RealFunction Df, real guess = 0) {
101 
102 
103  real x = guess;
104  unsigned int iter = 0;
105 
106  while(abs(f(x)) > OPTIMIZATION_TOL && iter <= OPTIMIZATION_NEWTON_ITER) {
107  x = x - (f(x) / Df(x));
108  iter++;
109  }
110 
111  if(iter > OPTIMIZATION_NEWTON_ITER) {
112  TH_MATH_ERROR("root_newton", x, NO_ALGO_CONVERGENCE);
113  return nan();
114  }
115 
116  return x;
117  }
118 
119 
128  inline real root_newton(dual(*f)(dual), real guess = 0) {
129 
130  real x = guess;
131  unsigned int iter = 0;
132 
133  dual s = f(dual(x, 1));
134 
135  while(abs(s.Re()) > OPTIMIZATION_TOL && iter <= OPTIMIZATION_NEWTON_ITER) {
136 
137  s = f(dual(x, 1));
138 
139  x = x - (s.Re() / s.Dual());
140  iter++;
141  }
142 
143  if(iter > OPTIMIZATION_NEWTON_ITER) {
144  TH_MATH_ERROR("root_newton", x, NO_ALGO_CONVERGENCE);
145  return nan();
146  }
147 
148  return x;
149  }
150 
151 
158  inline real root_newton(const polynomial<real>& p, real guess = 0) {
159 
160  real x = guess;
161  polynomial<real> Dp = deriv(p);
162  unsigned int iter = 0;
163 
164  while(abs(p(x)) > OPTIMIZATION_TOL && iter <= OPTIMIZATION_NEWTON_ITER) {
165  x = x - (p(x) / Dp(x));
166  iter++;
167  }
168 
169  if(iter > OPTIMIZATION_NEWTON_ITER) {
170  TH_MATH_ERROR("root_newton", x, NO_ALGO_CONVERGENCE);
171  return nan();
172  }
173 
174  return x;
175  }
176 
177 
190  inline complex<> root_newton(
191  complex<>(*f)(complex<>),
192  complex<>(*df)(complex<>),
193  complex<> guess = complex<>(0, 0),
194  real tolerance = OPTIMIZATION_TOL,
195  unsigned int max_iter = OPTIMIZATION_TOL) {
196 
197  complex<> z = guess;
198  unsigned int iter = 0;
199 
200  while(abs(f(z)) > tolerance && iter <= max_iter) {
201  z = z - (f(z) / df(z));
202  iter++;
203  }
204 
205  if(iter > max_iter) {
206  TH_MATH_ERROR("root_newton", z.Re(), NO_ALGO_CONVERGENCE);
207  return nan();
208  }
209 
210  return z;
211  }
212 
213 
222  template<typename RealFunction>
223  inline real root_halley(RealFunction f, RealFunction Df,
224  RealFunction D2f, real guess = 0) {
225 
226  real x = guess;
227  unsigned int iter = 0;
228 
229  while(abs(f(x)) > OPTIMIZATION_TOL && iter <= OPTIMIZATION_HALLEY_ITER) {
230  x = x - (2 * f(x) * Df(x)) / (2 * Df(x) - f(x) * D2f(x));
231  iter++;
232  }
233 
234  if(iter > OPTIMIZATION_HALLEY_ITER) {
235  TH_MATH_ERROR("root_halley", x, NO_ALGO_CONVERGENCE);
236  return nan();
237  }
238 
239  return x;
240  }
241 
242 
252  inline real root_halley(dual2(*f)(dual2), real guess = 0) {
253 
254  real x = guess;
255  unsigned int iter = 0;
256 
257  dual2 s = f(dual2(x, 1.0, 0.0));
258 
259  while(abs(s.Re()) > OPTIMIZATION_TOL && iter <= OPTIMIZATION_HALLEY_ITER) {
260 
261  s = f(dual2(x, 1, 0));
262 
263  x = x - (2 * s.Re() * s.Dual1())
264  / (2 * s.Dual1() - s.Re() * s.Dual2());
265  iter++;
266  }
267 
268  if(iter > OPTIMIZATION_HALLEY_ITER) {
269  TH_MATH_ERROR("root_halley", x, NO_ALGO_CONVERGENCE);
270  return nan();
271  }
272 
273  return x;
274  }
275 
276 
283  inline real root_halley(const polynomial<real>& p, real guess = 0) {
284 
285  polynomial<real> Dp = deriv(p);
286  polynomial<real> D2p = deriv(Dp);
287 
288  real x = guess;
289  unsigned int iter = 0;
290 
291  while(abs(p(x)) > OPTIMIZATION_TOL && iter <= OPTIMIZATION_HALLEY_ITER) {
292  x = x - (2 * p(x) * Dp(x)) / (2 * Dp(x) - p(x) * D2p(x));
293  iter++;
294  }
295 
296  if(iter > OPTIMIZATION_HALLEY_ITER) {
297  TH_MATH_ERROR("root_halley", x, NO_ALGO_CONVERGENCE);
298  return nan();
299  }
300 
301  return x;
302  }
303 
304 
311  template<typename RealFunction>
312  inline real root_steffensen(RealFunction f, real guess = 0) {
313 
314  real x = guess;
315  unsigned int iter = 0;
316 
317  while(abs(f(x)) > OPTIMIZATION_TOL && iter <= OPTIMIZATION_STEFFENSEN_ITER) {
318 
319  const real f_x = f(x);
320  x = x - (f_x / ((f(x + f_x) / f_x) - 1));
321  iter++;
322  }
323 
324  if(iter > OPTIMIZATION_STEFFENSEN_ITER) {
325  TH_MATH_ERROR("root_steffensen", x, NO_ALGO_CONVERGENCE);
326  return nan();
327  }
328 
329  return x;
330  }
331 
332 
339  inline real root_steffensen(const polynomial<real>& p, real guess = 0) {
340 
341  real x = guess;
342  unsigned int iter = 0;
343 
344  while(abs(p(x)) > OPTIMIZATION_TOL && iter <= OPTIMIZATION_STEFFENSEN_ITER) {
345  x = x - (p(x) / ((p(x + p(x)) / p(x)) - 1));
346  iter++;
347  }
348 
349  if(iter > OPTIMIZATION_STEFFENSEN_ITER) {
350  TH_MATH_ERROR("root_steffensen", x, NO_ALGO_CONVERGENCE);
351  return nan();
352  }
353 
354  return x;
355  }
356 
357 
366  template<typename RealFunction>
367  inline real root_chebyshev(
368  RealFunction f, RealFunction Df,
369  RealFunction D2f, real guess = 0) {
370 
371  real x = guess;
372  unsigned int iter = 0;
373 
374  while(abs(f(x)) > OPTIMIZATION_TOL && iter <= OPTIMIZATION_CHEBYSHEV_ITER) {
375  x = x - (f(x) / Df(x)) - square((f(x) / Df(x))) * (Df(x) / (D2f(x) * 2));
376  iter++;
377  }
378 
379  if(iter > OPTIMIZATION_CHEBYSHEV_ITER) {
380  TH_MATH_ERROR("root_chebyshev", x, NO_ALGO_CONVERGENCE);
381  return nan();
382  }
383 
384  return x;
385  }
386 
387 
396  inline real root_chebyshev(dual2(*f)(dual2), real guess = 0) {
397 
398  real x = guess;
399  unsigned int iter = 0;
400 
401  dual2 s = f(dual2(x, 1.0, 0.0));
402 
403  while(abs(s.Re()) > OPTIMIZATION_TOL && iter <= OPTIMIZATION_CHEBYSHEV_ITER) {
404 
405  s = f(dual2(x, 1, 0));
406 
407  x = x - (s.Re() / s.Dual1())
408  - square((s.Re() / s.Dual1()))
409  * (s.Dual1() / (s.Dual2() * 2));
410  iter++;
411  }
412 
413  if(iter > OPTIMIZATION_CHEBYSHEV_ITER) {
414  TH_MATH_ERROR("root_chebyshev", x, NO_ALGO_CONVERGENCE);
415  return nan();
416  }
417 
418  return x;
419  }
420 
421 
428  inline real root_chebyshev(const polynomial<real>& p, real guess = 0) {
429 
430  polynomial<real> Dp = deriv(p);
431  polynomial<real> D2p = deriv(p);
432 
433  real x = guess;
434  unsigned int iter = 0;
435 
436  while(abs(p(x)) > OPTIMIZATION_TOL && iter <= OPTIMIZATION_CHEBYSHEV_ITER) {
437  x = x - (p(x) / Dp(x)) - square((p(x) / Dp(x))) * (Dp(x) / (D2p(x) * 2));
438  iter++;
439  }
440 
441  if(iter > OPTIMIZATION_CHEBYSHEV_ITER) {
442  TH_MATH_ERROR("root_chebyshev", x, NO_ALGO_CONVERGENCE);
443  return nan();
444  }
445 
446  return x;
447  }
448 
449 
461  template<typename RealFunction>
462  inline std::vector<real> roots(
463  RealFunction f, real a, real b,
464  real tolerance = OPTIMIZATION_TOL, real steps = 10) {
465 
466  if(steps == 0) {
467  TH_MATH_ERROR("roots", steps, DIV_BY_ZERO);
468  return {nan()};
469  }
470 
471  // Find candidate intervals
472  std::vector<vec2> intervals = find_root_intervals(f, a, b, steps);
473 
474  std::vector<real> res;
475  res.reserve(intervals.size());
476 
477  // Iterate through all candidate intervals and refine the results
478  for (unsigned int i = 0; i < intervals.size(); ++i) {
479 
480  // Check whether the extremes of the candidate intervals
481  // happen to coincide with the roots
482  if(abs(f(intervals[i][0])) <= MACH_EPSILON) {
483  res.push_back(intervals[i][0]);
484  continue;
485  }
486 
487  if(abs(f(intervals[i][1])) <= MACH_EPSILON) {
488  res.push_back(intervals[i][1]);
489  continue;
490  }
491 
492  // Approximate the roots using bisection inside each interval
493  res.push_back(
494  root_bisection(f, intervals[i][0], intervals[i][1], tolerance));
495  }
496 
497  return res;
498  }
499 
500 
508  template<typename Field>
509  inline std::vector<Field> roots(
510  const polynomial<Field>& p,
511  real tolerance = OPTIMIZATION_TOL,
512  unsigned int steps = 0) {
513 
514  // Effective order of the polynomial
515  const unsigned int n = p.find_order();
516  p /= p.coeff[n];
517 
518  // Absolute value of the highest coefficient
519  Field a_hi = abs(p.coeff[n]);
520  Field a_sum = 0;
521 
522  // Sum the absolute values of the lesser coefficients
523  for (unsigned int i = 0; i < n; ++i)
524  a_sum += abs(p.coeff[i]);
525 
526  // The roots are bounded in absolute value by the maximum
527  const Field M = max(a_hi, a_sum);
528 
529  // Back track from the bounds to the first sign inversion ?
530 
531  return roots(
532  [p](real x) { return p(x); },
533  -M, M, tolerance, steps ? steps : (2 * n));
534  }
535 
536 }
537 
538 #endif
theoretica::complex
Complex number in algebraic form .
Definition: complex.h:26
theoretica::complex::Re
Type Re() const
Get the real part of the complex number.
Definition: complex.h:74
theoretica::dual2
Second order dual number class.
Definition: dual2.h:29
theoretica::dual2::Dual2
real Dual2() const
Return second order dual part.
Definition: dual2.h:95
theoretica::dual2::Re
real Re() const
Return real part.
Definition: dual2.h:85
theoretica::dual2::Dual1
real Dual1() const
Return first order dual part.
Definition: dual2.h:90
theoretica::dual
Dual number class.
Definition: dual.h:28
theoretica::dual::Re
const real & Re() const
Return real part.
Definition: dual.h:75
theoretica::dual::Dual
const real & Dual() const
Return dual part.
Definition: dual.h:96
theoretica::polynomial
A polynomial of arbitrary order.
Definition: polynomial.h:25
TH_MATH_ERROR
#define TH_MATH_ERROR(F_NAME, VALUE, EXCEPTION)
TH_MATH_ERROR is a macro which throws exceptions or modifies errno (depending on which compiling opti...
Definition: error.h:219
theoretica::autodiff::deriv
real deriv(DualFunction f, real x)
Compute the derivative of a function at the given point using univariate automatic differentiation.
Definition: autodiff.h:41
theoretica::parallel::square
Vector square(const Vector &v)
Parallel element-wise evaluation of the square function.
Definition: parallel.h:53
theoretica::parallel::abs
Vector abs(const Vector &v)
Parallel element-wise evaluation of the abs function.
Definition: parallel.h:89
theoretica
Main namespace of the library which contains all functions and objects.
Definition: algebra.h:27
theoretica::root_halley
real root_halley(const polynomial< real > &p, real guess=0)
Approximate a root of a polynomial using Halley's method.
Definition: roots.h:283
theoretica::real
double real
A real number, defined as a floating point type.
Definition: constants.h:188
theoretica::root_steffensen
real root_steffensen(const polynomial< real > &p, real guess=0)
Approximate a root of a polynomial using Steffensen's method.
Definition: roots.h:339
theoretica::OPTIMIZATION_TOL
constexpr real OPTIMIZATION_TOL
Approximation tolerance for root finding.
Definition: constants.h:278
theoretica::max
auto max(const Vector &X)
Finds the maximum value inside a dataset.
Definition: dataset.h:330
theoretica::MACH_EPSILON
constexpr real MACH_EPSILON
Machine epsilon for the real type.
Definition: constants.h:197
theoretica::roots
std::vector< Field > roots(const polynomial< Field > &p, real tolerance=OPTIMIZATION_TOL, unsigned int steps=0)
Find all the roots of a polynomial.
Definition: roots.h:509
theoretica::root_newton
complex root_newton(complex<>(*f)(complex<>), complex<>(*df)(complex<>), complex<> guess=complex<>(0, 0), real tolerance=OPTIMIZATION_TOL, unsigned int max_iter=OPTIMIZATION_TOL)
Approximate a root of an arbitrary complex function using Newton's method,.
Definition: roots.h:190
theoretica::find_root_intervals
std::vector< vec2 > find_root_intervals(RealFunction f, real a, real b, unsigned int steps=10)
Find candidate intervals for root finding.
Definition: roots.h:27
theoretica::OPTIMIZATION_CHEBYSHEV_ITER
constexpr unsigned int OPTIMIZATION_CHEBYSHEV_ITER
Maximum number of iterations for the Chebyshev algorithm.
Definition: constants.h:296
theoretica::OPTIMIZATION_BISECTION_ITER
constexpr unsigned int OPTIMIZATION_BISECTION_ITER
Maximum number of iterations for the bisection algorithm.
Definition: constants.h:281
theoretica::OPTIMIZATION_STEFFENSEN_ITER
constexpr unsigned int OPTIMIZATION_STEFFENSEN_ITER
Maximum number of iterations for the Steffensen algorithm.
Definition: constants.h:293
theoretica::OPTIMIZATION_NEWTON_ITER
constexpr unsigned int OPTIMIZATION_NEWTON_ITER
Maximum number of iterations for the Newton-Raphson algorithm.
Definition: constants.h:290
theoretica::OPTIMIZATION_HALLEY_ITER
constexpr unsigned int OPTIMIZATION_HALLEY_ITER
Maximum number of iterations for Halley's method.
Definition: constants.h:287
theoretica::nan
real nan()
Return a quiet NaN number in floating point representation.
Definition: error.h:54
theoretica::root_chebyshev
real root_chebyshev(const polynomial< real > &p, real guess=0)
Approximate a root of a polynomial using Chebyshev's method.
Definition: roots.h:428
theoretica::root_bisection
real root_bisection(RealFunction f, real a, real b, real tolerance=OPTIMIZATION_TOL)
Approximate a root of an arbitrary function using bisection inside a compact interval [a,...
Definition: roots.h:56