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