Eigen-unsupported  5.0.1-dev
 
Loading...
Searching...
No Matches
Polynomials
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4//
5// This Source Code Form is subject to the terms of the Mozilla
6// Public License v. 2.0. If a copy of the MPL was not distributed
7// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
8
9#ifndef EIGEN_POLYNOMIALS_MODULE_H
10#define EIGEN_POLYNOMIALS_MODULE_H
11
12#include "../../Eigen/Core"
13
14#include "../../Eigen/Eigenvalues"
15
16#include "../../Eigen/src/Core/util/DisableStupidWarnings.h"
17
18// Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module
19#if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS >= 2)
20#ifndef EIGEN_HIDE_HEAVY_CODE
21#define EIGEN_HIDE_HEAVY_CODE
22#endif
23#elif defined EIGEN_HIDE_HEAVY_CODE
24#undef EIGEN_HIDE_HEAVY_CODE
25#endif
26
27/**
28 * \defgroup Polynomials_Module Polynomials module
29 * \brief This module provides a QR based polynomial solver.
30 *
31 * To use this module, add
32 * \code
33 * #include <unsupported/Eigen/Polynomials>
34 * \endcode
35 * at the start of your source file.
36 */
37
38// IWYU pragma: begin_exports
39#include "src/Polynomials/PolynomialUtils.h"
40#include "src/Polynomials/Companion.h"
41#include "src/Polynomials/PolynomialSolver.h"
42// IWYU pragma: end_exports
43
44/**
45 \page polynomials Polynomials defines functions for dealing with polynomials
46 and a QR based polynomial solver.
47 \ingroup Polynomials_Module
48
49 The remainder of the page documents first the functions for evaluating, computing
50 polynomials, computing estimates about polynomials and next the QR based polynomial
51 solver.
52
53 \section polynomialUtils convenient functions to deal with polynomials
54 \subsection roots_to_monicPolynomial
55 The function
56 \code
57 void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
58 \endcode
59 computes the coefficients \f$ a_i \f$ of
60
61 \f$ p(x) = a_0 + a_{1}x + ... + a_{n-1}x^{n-1} + x^n \f$
62
63 where \f$ p \f$ is known through its roots i.e. \f$ p(x) = (x-r_1)(x-r_2)...(x-r_n) \f$.
64
65 \subsection poly_eval
66 The function
67 \code
68 T poly_eval( const Polynomials& poly, const T& x )
69 \endcode
70 evaluates a polynomial at a given point using stabilized H&ouml;rner method.
71
72 The following code: first computes the coefficients in the monomial basis of the monic polynomial that has the
73 provided roots; then, it evaluates the computed polynomial, using a stabilized H&ouml;rner method.
74
75 \include PolynomialUtils1.cpp
76 Output: \verbinclude PolynomialUtils1.out
77
78 \subsection Cauchy bounds
79 The function
80 \code
81 Real cauchy_max_bound( const Polynomial& poly )
82 \endcode
83 provides a maximum bound (the Cauchy one: \f$C(p)\f$) for the absolute value of a root of the given polynomial
84 i.e. \f$ \forall r_i \f$ root of \f$ p(x) = \sum_{k=0}^d a_k x^k \f$, \f$ |r_i| \le C(p) = \sum_{k=0}^{d} \left |
85 \frac{a_k}{a_d} \right | \f$ The leading coefficient \f$ p \f$: should be non zero \f$a_d \neq 0\f$.
86
87
88 The function
89 \code
90 Real cauchy_min_bound( const Polynomial& poly )
91 \endcode
92 provides a minimum bound (the Cauchy one: \f$c(p)\f$) for the absolute value of a non zero root of the given
93 polynomial i.e. \f$ \forall r_i \neq 0 \f$ root of \f$ p(x) = \sum_{k=0}^d a_k x^k \f$, \f$ |r_i| \ge c(p) = \left(
94 \sum_{k=0}^{d} \left | \frac{a_k}{a_0} \right | \right)^{-1} \f$
95
96
97
98
99 \section QR polynomial solver class
100 Computes the complex roots of a polynomial by computing the eigenvalues of the associated companion matrix with
101 the QR algorithm.
102
103 The roots of \f$ p(x) = a_0 + a_1 x + a_2 x^2 + a_{3} x^3 + x^4 \f$ are the eigenvalues of
104 \f$
105 \left [
106 \begin{array}{cccc}
107 0 & 0 & 0 & a_0 \\
108 1 & 0 & 0 & a_1 \\
109 0 & 1 & 0 & a_2 \\
110 0 & 0 & 1 & a_3
111 \end{array} \right ]
112 \f$
113
114 However, the QR algorithm is not guaranteed to converge when there are several eigenvalues with same modulus.
115
116 Therefore the current polynomial solver is guaranteed to provide a correct result only when the complex roots
117 \f$r_1,r_2,...,r_d\f$ have distinct moduli i.e.
118
119 \f$ \forall i,j \in [1;d],~ \| r_i \| \neq \| r_j \| \f$.
120
121 With 32bit (float) floating types this problem shows up frequently.
122 However, almost always, correct accuracy is reached even in these cases for 64bit
123 (double) floating types and small polynomial degree (<20).
124
125 \include PolynomialSolver1.cpp
126
127 In the above example:
128
129 -# a simple use of the polynomial solver is shown;
130 -# the accuracy problem with the QR algorithm is presented: a polynomial with almost conjugate roots is provided
131 to the solver. Those roots have almost same module therefore the QR algorithm failed to converge: the accuracy of the
132 last root is bad;
133 -# a simple way to circumvent the problem is shown: use doubles instead of floats.
134
135 Output: \verbinclude PolynomialSolver1.out
136*/
137
138#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
139
140#endif // EIGEN_POLYNOMIALS_MODULE_H