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