diff options
author | Benoit Jacob <jacob.benoit.1@gmail.com> | 2010-04-16 11:25:50 -0400 |
---|---|---|
committer | Benoit Jacob <jacob.benoit.1@gmail.com> | 2010-04-16 11:25:50 -0400 |
commit | 0ab431d7b860afc6766c7c20f7bb39a1d71bff62 (patch) | |
tree | f8da6ce3cc7738735f315f7954bbbabf48e0c621 /unsupported/Eigen | |
parent | ff6a46105d86e92753858c1b2aea8bcaf4575819 (diff) | |
parent | ea1a2df37092f88f5594dfea1f7e4996dd8e612d (diff) |
* merge with mainline
* adapt Eigenvalues module to the new rule that the RowMajorBit must have the proper value for vectors
* Fix RowMajorBit in ei_traits<ProductBase>
* Fix vectorizability logic in CoeffBasedProduct
Diffstat (limited to 'unsupported/Eigen')
-rw-r--r-- | unsupported/Eigen/CMakeLists.txt | 2 | ||||
-rw-r--r-- | unsupported/Eigen/FFT | 2 | ||||
-rw-r--r-- | unsupported/Eigen/MatrixFunctions | 213 | ||||
-rw-r--r-- | unsupported/Eigen/Polynomials | 137 | ||||
-rw-r--r-- | unsupported/Eigen/src/CMakeLists.txt | 1 | ||||
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h | 50 | ||||
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h | 96 | ||||
-rw-r--r-- | unsupported/Eigen/src/Polynomials/CMakeLists.txt | 6 | ||||
-rw-r--r-- | unsupported/Eigen/src/Polynomials/Companion.h | 281 | ||||
-rw-r--r-- | unsupported/Eigen/src/Polynomials/PolynomialSolver.h | 395 | ||||
-rw-r--r-- | unsupported/Eigen/src/Polynomials/PolynomialUtils.h | 153 |
11 files changed, 1184 insertions, 152 deletions
diff --git a/unsupported/Eigen/CMakeLists.txt b/unsupported/Eigen/CMakeLists.txt index f9d79c51a..87cc4be1e 100644 --- a/unsupported/Eigen/CMakeLists.txt +++ b/unsupported/Eigen/CMakeLists.txt @@ -1,4 +1,4 @@ -set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3) +set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3 Polynomials) install(FILES ${Eigen_HEADERS} diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT index 045c44611..2a1e21007 100644 --- a/unsupported/Eigen/FFT +++ b/unsupported/Eigen/FFT @@ -320,7 +320,7 @@ class FFT // if the vector is strided, then we need to copy it to a packed temporary Matrix<src_type,1,Dynamic> tmp; if ( resize_input ) { - size_t ncopy = min(src.size(),src.size() + resize_input); + size_t ncopy = std::min(src.size(),src.size() + resize_input); tmp.setZero(src.size() + resize_input); if ( realfft && HasFlag(HalfSpectrum) ) { // pad at the Nyquist bin diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index 292357fda..e0bc4732c 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -40,6 +40,22 @@ namespace Eigen { * \brief This module aims to provide various methods for the computation of * matrix functions. * + * To use this module, add + * \code + * #include <unsupported/Eigen/MatrixFunctions> + * \endcode + * at the start of your source file. + * + * This module defines the following MatrixBase methods. + * - \ref matrixbase_cos "MatrixBase::cos()", for computing the matrix cosine + * - \ref matrixbase_cosh "MatrixBase::cosh()", for computing the matrix hyperbolic cosine + * - \ref matrixbase_exp "MatrixBase::exp()", for computing the matrix exponential + * - \ref matrixbase_matrixfunction "MatrixBase::matrixFunction()", for computing general matrix functions + * - \ref matrixbase_sin "MatrixBase::sin()", for computing the matrix sine + * - \ref matrixbase_sinh "MatrixBase::sinh()", for computing the matrix hyperbolic sine + * + * These methods are the main entry points to this module. + * * %Matrix functions are defined as follows. Suppose that \f$ f \f$ * is an entire function (that is, a function on the complex plane * that is everywhere complex differentiable). Then its Taylor @@ -49,16 +65,205 @@ namespace Eigen { * function by the same series: * \f[ f(M) = f(0) + f'(0) M + \frac{f''(0)}{2} M^2 + \frac{f'''(0)}{3!} M^3 + \cdots \f] * - * \code - * #include <unsupported/Eigen/MatrixFunctions> - * \endcode */ #include "src/MatrixFunctions/MatrixExponential.h" #include "src/MatrixFunctions/MatrixFunction.h" -} +/** +\page matrixbaseextra MatrixBase methods defined in the MatrixFunctions module +\ingroup MatrixFunctions_Module + +The remainder of the page documents the following MatrixBase methods +which are defined in the MatrixFunctions module. + + + +\section matrixbase_cos MatrixBase::cos() + +Compute the matrix cosine. + +\code +const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const +\endcode + +\param[in] M a square matrix. +\returns expression representing \f$ \cos(M) \f$. + +This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cos(). + +\sa \ref matrixbase_sin "sin()" for an example. + + + +\section matrixbase_cosh MatrixBase::cosh() + +Compute the matrix hyberbolic cosine. + +\code +const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const +\endcode + +\param[in] M a square matrix. +\returns expression representing \f$ \cosh(M) \f$ + +This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::cosh(). + +\sa \ref matrixbase_sinh "sinh()" for an example. + + + +\section matrixbase_exp MatrixBase::exp() + +Compute the matrix exponential. + +\code +const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const +\endcode + +\param[in] M matrix whose exponential is to be computed. +\returns expression representing the matrix exponential of \p M. + +The matrix exponential of \f$ M \f$ is defined by +\f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f] +The matrix exponential can be used to solve linear ordinary +differential equations: the solution of \f$ y' = My \f$ with the +initial condition \f$ y(0) = y_0 \f$ is given by +\f$ y(t) = \exp(M) y_0 \f$. + +The cost of the computation is approximately \f$ 20 n^3 \f$ for +matrices of size \f$ n \f$. The number 20 depends weakly on the +norm of the matrix. + +The matrix exponential is computed using the scaling-and-squaring +method combined with Padé approximation. The matrix is first +rescaled, then the exponential of the reduced matrix is computed +approximant, and then the rescaling is undone by repeated +squaring. The degree of the Padé approximant is chosen such +that the approximation error is less than the round-off +error. However, errors may accumulate during the squaring phase. + +Details of the algorithm can be found in: Nicholas J. Higham, "The +scaling and squaring method for the matrix exponential revisited," +<em>SIAM J. %Matrix Anal. Applic.</em>, <b>26</b>:1179–1193, +2005. + +Example: The following program checks that +\f[ \exp \left[ \begin{array}{ccc} + 0 & \frac14\pi & 0 \\ + -\frac14\pi & 0 & 0 \\ + 0 & 0 & 0 + \end{array} \right] = \left[ \begin{array}{ccc} + \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ + \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ + 0 & 0 & 1 + \end{array} \right]. \f] +This corresponds to a rotation of \f$ \frac14\pi \f$ radians around +the z-axis. + +\include MatrixExponential.cpp +Output: \verbinclude MatrixExponential.out + +\note \p M has to be a matrix of \c float, \c double, +\c complex<float> or \c complex<double> . + + + +\section matrixbase_matrixfunction MatrixBase::matrixFunction() + +Compute a matrix function. + +\code +const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f) const +\endcode + +\param[in] M argument of matrix function, should be a square matrix. +\param[in] f an entire function; \c f(x,n) should compute the n-th +derivative of f at x. +\returns expression representing \p f applied to \p M. + +Suppose that \p M is a matrix whose entries have type \c Scalar. +Then, the second argument, \p f, should be a function with prototype +\code +ComplexScalar f(ComplexScalar, int) +\endcode +where \c ComplexScalar = \c std::complex<Scalar> if \c Scalar is +real (e.g., \c float or \c double) and \c ComplexScalar = +\c Scalar if \c Scalar is complex. The return value of \c f(x,n) +should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x. + +This routine uses the algorithm described in: +Philip Davies and Nicholas J. Higham, +"A Schur-Parlett algorithm for computing matrix functions", +<em>SIAM J. %Matrix Anal. Applic.</em>, <b>25</b>:464–485, 2003. + +The actual work is done by the MatrixFunction class. + +Example: The following program checks that +\f[ \exp \left[ \begin{array}{ccc} + 0 & \frac14\pi & 0 \\ + -\frac14\pi & 0 & 0 \\ + 0 & 0 & 0 + \end{array} \right] = \left[ \begin{array}{ccc} + \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ + \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ + 0 & 0 & 1 + \end{array} \right]. \f] +This corresponds to a rotation of \f$ \frac14\pi \f$ radians around +the z-axis. This is the same example as used in the documentation +of \ref matrixbase_exp "exp()". + +\include MatrixFunction.cpp +Output: \verbinclude MatrixFunction.out + +Note that the function \c expfn is defined for complex numbers +\c x, even though the matrix \c A is over the reals. Instead of +\c expfn, we could also have used StdStemFunctions::exp: +\code +A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B); +\endcode + + + +\section matrixbase_sin MatrixBase::sin() + +Compute the matrix sine. + +\code +const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const +\endcode + +\param[in] M a square matrix. +\returns expression representing \f$ \sin(M) \f$. + +This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sin(). + +Example: \include MatrixSine.cpp +Output: \verbinclude MatrixSine.out + + + +\section matrixbase_sinh const MatrixBase::sinh() + +Compute the matrix hyperbolic sine. + +\code +MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const +\endcode + +\param[in] M a square matrix. +\returns expression representing \f$ \sinh(M) \f$ + +This function calls \ref matrixbase_matrixfunction "matrixFunction()" with StdStemFunctions::sinh(). + +Example: \include MatrixSinh.cpp +Output: \verbinclude MatrixSinh.out + +*/ + +} + #endif // EIGEN_MATRIX_FUNCTIONS diff --git a/unsupported/Eigen/Polynomials b/unsupported/Eigen/Polynomials new file mode 100644 index 000000000..d69abb1be --- /dev/null +++ b/unsupported/Eigen/Polynomials @@ -0,0 +1,137 @@ +#ifndef EIGEN_POLYNOMIALS_MODULE_H +#define EIGEN_POLYNOMIALS_MODULE_H + +#include <Eigen/Core> + +#include <Eigen/src/Core/util/DisableMSVCWarnings.h> + +#include <Eigen/QR> + +// Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module +#if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2) + #ifndef EIGEN_HIDE_HEAVY_CODE + #define EIGEN_HIDE_HEAVY_CODE + #endif +#elif defined EIGEN_HIDE_HEAVY_CODE + #undef EIGEN_HIDE_HEAVY_CODE +#endif + +namespace Eigen { + +/** \ingroup Unsupported_modules + * \defgroup Polynomials_Module Polynomials module + * + * \nonstableyet + * + * \brief This module provides a QR based polynomial solver. + * + * To use this module, add + * \code + * #include <unsupported/Eigen/Polynomials> + * \endcode + * at the start of your source file. + */ + +#include "src/Polynomials/PolynomialUtils.h" +#include "src/Polynomials/Companion.h" +#include "src/Polynomials/PolynomialSolver.h" + +/** + \page polynomials Polynomials defines functions for dealing with polynomials + and a QR based polynomial solver. + \ingroup Polynomials_Module + + The remainder of the page documents first the functions for evaluating, computing + polynomials, computing estimates about polynomials and next the QR based polynomial + solver. + + \section polynomialUtils convenient functions to deal with polynomials + \subsection roots_to_monicPolynomial + The function + \code + void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly ) + \endcode + computes the coefficients \f$ a_i \f$ of + + \f$ p(x) = a_0 + a_{1}x + ... + a_{n-1}x^{n-1} + x^n \f$ + + where \f$ p \f$ is known through its roots i.e. \f$ p(x) = (x-r_1)(x-r_2)...(x-r_n) \f$. + + \subsection poly_eval + The function + \code + T poly_eval( const Polynomials& poly, const T& x ) + \endcode + evaluates a polynomial at a given point using stabilized Hörner method. + + The following code: first computes the coefficients in the monomial basis of the monic polynomial that has the provided roots; + then, it evaluates the computed polynomial, using a stabilized Hörner method. + + \include PolynomialUtils1.cpp + Output: \verbinclude PolynomialUtils1.out + + \subsection Cauchy bounds + The function + \code + Real cauchy_max_bound( const Polynomial& poly ) + \endcode + provides a maximum bound (the Cauchy one: \f$C(p)\f$) for the absolute value of a root of the given polynomial 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 | \frac{a_k}{a_d} \right | \f$ + The leading coefficient \f$ p \f$: should be non zero \f$a_d \neq 0\f$. + + + The function + \code + Real cauchy_min_bound( const Polynomial& poly ) + \endcode + 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. + \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( \sum_{k=0}^{d} \left | \frac{a_k}{a_0} \right | \right)^{-1} \f$ + + + + + \section QR polynomial solver class + Computes the complex roots of a polynomial by computing the eigenvalues of the associated companion matrix with the QR algorithm. + + 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 + \f$ + \left [ + \begin{array}{cccc} + 0 & 0 & 0 & a_0 \\ + 1 & 0 & 0 & a_1 \\ + 0 & 1 & 0 & a_2 \\ + 0 & 0 & 1 & a_3 + \end{array} \right ] + \f$ + + However, the QR algorithm is not guaranteed to converge when there are several eigenvalues with same modulus. + + 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. + + \f$ \forall i,j \in [1;d],~ \| r_i \| \neq \| r_j \| \f$. + + With 32bit (float) floating types this problem shows up frequently. + However, almost always, correct accuracy is reached even in these cases for 64bit + (double) floating types and small polynomial degree (<20). + + \include PolynomialSolver1.cpp + + In the above example: + + -# a simple use of the polynomial solver is shown; + -# the accuracy problem with the QR algorithm is presented: a polynomial with almost conjugate roots is provided to the solver. + Those roots have almost same module therefore the QR algorithm failed to converge: the accuracy + of the last root is bad; + -# a simple way to circumvent the problem is shown: use doubles instead of floats. + + Output: \verbinclude PolynomialSolver1.out +*/ + +} // namespace Eigen + +#include <Eigen/src/Core/util/EnableMSVCWarnings.h> + +#endif // EIGEN_POLYNOMIALS_MODULE_H +/* vim: set filetype=cpp et sw=2 ts=2 ai: */ diff --git a/unsupported/Eigen/src/CMakeLists.txt b/unsupported/Eigen/src/CMakeLists.txt index a0870d3af..035a84ca8 100644 --- a/unsupported/Eigen/src/CMakeLists.txt +++ b/unsupported/Eigen/src/CMakeLists.txt @@ -5,3 +5,4 @@ ADD_SUBDIRECTORY(MoreVectorization) # ADD_SUBDIRECTORY(FFT) # ADD_SUBDIRECTORY(Skyline) ADD_SUBDIRECTORY(MatrixFunctions) +ADD_SUBDIRECTORY(Polynomials) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index 3ef91a7b2..d7eb1b2fe 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -330,56 +330,6 @@ struct ei_traits<MatrixExponentialReturnValue<Derived> > typedef typename Derived::PlainObject ReturnType; }; -/** \ingroup MatrixFunctions_Module - * - * \brief Compute the matrix exponential. - * - * \param[in] M matrix whose exponential is to be computed. - * \returns expression representing the matrix exponential of \p M. - * - * The matrix exponential of \f$ M \f$ is defined by - * \f[ \exp(M) = \sum_{k=0}^\infty \frac{M^k}{k!}. \f] - * The matrix exponential can be used to solve linear ordinary - * differential equations: the solution of \f$ y' = My \f$ with the - * initial condition \f$ y(0) = y_0 \f$ is given by - * \f$ y(t) = \exp(M) y_0 \f$. - * - * The cost of the computation is approximately \f$ 20 n^3 \f$ for - * matrices of size \f$ n \f$. The number 20 depends weakly on the - * norm of the matrix. - * - * The matrix exponential is computed using the scaling-and-squaring - * method combined with Padé approximation. The matrix is first - * rescaled, then the exponential of the reduced matrix is computed - * approximant, and then the rescaling is undone by repeated - * squaring. The degree of the Padé approximant is chosen such - * that the approximation error is less than the round-off - * error. However, errors may accumulate during the squaring phase. - * - * Details of the algorithm can be found in: Nicholas J. Higham, "The - * scaling and squaring method for the matrix exponential revisited," - * <em>SIAM J. %Matrix Anal. Applic.</em>, <b>26</b>:1179–1193, - * 2005. - * - * Example: The following program checks that - * \f[ \exp \left[ \begin{array}{ccc} - * 0 & \frac14\pi & 0 \\ - * -\frac14\pi & 0 & 0 \\ - * 0 & 0 & 0 - * \end{array} \right] = \left[ \begin{array}{ccc} - * \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ - * \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ - * 0 & 0 & 1 - * \end{array} \right]. \f] - * This corresponds to a rotation of \f$ \frac14\pi \f$ radians around - * the z-axis. - * - * \include MatrixExponential.cpp - * Output: \verbinclude MatrixExponential.out - * - * \note \p M has to be a matrix of \c float, \c double, - * \c complex<float> or \c complex<double> . - */ template <typename Derived> const MatrixExponentialReturnValue<Derived> MatrixBase<Derived>::exp() const { diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h index 72e42aa7c..270d23b8b 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h @@ -536,56 +536,6 @@ struct ei_traits<MatrixFunctionReturnValue<Derived> > /********** MatrixBase methods **********/ -/** \ingroup MatrixFunctions_Module - * - * \brief Compute a matrix function. - * - * \param[in] M argument of matrix function, should be a square matrix. - * \param[in] f an entire function; \c f(x,n) should compute the n-th - * derivative of f at x. - * \returns expression representing \p f applied to \p M. - * - * Suppose that \p M is a matrix whose entries have type \c Scalar. - * Then, the second argument, \p f, should be a function with prototype - * \code - * ComplexScalar f(ComplexScalar, int) - * \endcode - * where \c ComplexScalar = \c std::complex<Scalar> if \c Scalar is - * real (e.g., \c float or \c double) and \c ComplexScalar = - * \c Scalar if \c Scalar is complex. The return value of \c f(x,n) - * should be \f$ f^{(n)}(x) \f$, the n-th derivative of f at x. - * - * This routine uses the algorithm described in: - * Philip Davies and Nicholas J. Higham, - * "A Schur-Parlett algorithm for computing matrix functions", - * <em>SIAM J. %Matrix Anal. Applic.</em>, <b>25</b>:464–485, 2003. - * - * The actual work is done by the MatrixFunction class. - * - * Example: The following program checks that - * \f[ \exp \left[ \begin{array}{ccc} - * 0 & \frac14\pi & 0 \\ - * -\frac14\pi & 0 & 0 \\ - * 0 & 0 & 0 - * \end{array} \right] = \left[ \begin{array}{ccc} - * \frac12\sqrt2 & -\frac12\sqrt2 & 0 \\ - * \frac12\sqrt2 & \frac12\sqrt2 & 0 \\ - * 0 & 0 & 1 - * \end{array} \right]. \f] - * This corresponds to a rotation of \f$ \frac14\pi \f$ radians around - * the z-axis. This is the same example as used in the documentation - * of MatrixBase::exp(). - * - * \include MatrixFunction.cpp - * Output: \verbinclude MatrixFunction.out - * - * Note that the function \c expfn is defined for complex numbers - * \c x, even though the matrix \c A is over the reals. Instead of - * \c expfn, we could also have used StdStemFunctions::exp: - * \code - * A.matrixFunction(StdStemFunctions<std::complex<double> >::exp, &B); - * \endcode - */ template <typename Derived> const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f) const { @@ -593,18 +543,6 @@ const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::matrixFunction(typ return MatrixFunctionReturnValue<Derived>(derived(), f); } -/** \ingroup MatrixFunctions_Module - * - * \brief Compute the matrix sine. - * - * \param[in] M a square matrix. - * \returns expression representing \f$ \sin(M) \f$. - * - * This function calls matrixFunction() with StdStemFunctions::sin(). - * - * \include MatrixSine.cpp - * Output: \verbinclude MatrixSine.out - */ template <typename Derived> const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const { @@ -613,17 +551,6 @@ const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sin() const return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sin); } -/** \ingroup MatrixFunctions_Module - * - * \brief Compute the matrix cosine. - * - * \param[in] M a square matrix. - * \returns expression representing \f$ \cos(M) \f$. - * - * This function calls matrixFunction() with StdStemFunctions::cos(). - * - * \sa ei_matrix_sin() for an example. - */ template <typename Derived> const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const { @@ -632,18 +559,6 @@ const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cos() const return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::cos); } -/** \ingroup MatrixFunctions_Module - * - * \brief Compute the matrix hyperbolic sine. - * - * \param[in] M a square matrix. - * \returns expression representing \f$ \sinh(M) \f$ - * - * This function calls matrixFunction() with StdStemFunctions::sinh(). - * - * \include MatrixSinh.cpp - * Output: \verbinclude MatrixSinh.out - */ template <typename Derived> const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const { @@ -652,17 +567,6 @@ const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::sinh() const return MatrixFunctionReturnValue<Derived>(derived(), StdStemFunctions<ComplexScalar>::sinh); } -/** \ingroup MatrixFunctions_Module - * - * \brief Compute the matrix hyberbolic cosine. - * - * \param[in] M a square matrix. - * \returns expression representing \f$ \cosh(M) \f$ - * - * This function calls matrixFunction() with StdStemFunctions::cosh(). - * - * \sa ei_matrix_sinh() for an example. - */ template <typename Derived> const MatrixFunctionReturnValue<Derived> MatrixBase<Derived>::cosh() const { diff --git a/unsupported/Eigen/src/Polynomials/CMakeLists.txt b/unsupported/Eigen/src/Polynomials/CMakeLists.txt new file mode 100644 index 000000000..51f13f3cb --- /dev/null +++ b/unsupported/Eigen/src/Polynomials/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_Polynomials_SRCS "*.h") + +INSTALL(FILES + ${Eigen_Polynomials_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/Polynomials COMPONENT Devel + ) diff --git a/unsupported/Eigen/src/Polynomials/Companion.h b/unsupported/Eigen/src/Polynomials/Companion.h new file mode 100644 index 000000000..7c9e9c518 --- /dev/null +++ b/unsupported/Eigen/src/Polynomials/Companion.h @@ -0,0 +1,281 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_COMPANION_H +#define EIGEN_COMPANION_H + +// This file requires the user to include +// * Eigen/Core +// * Eigen/src/PolynomialSolver.h + +#ifndef EIGEN_PARSED_BY_DOXYGEN + +template <typename T> +T ei_radix(){ return 2; } + +template <typename T> +T ei_radix2(){ return ei_radix<T>()*ei_radix<T>(); } + +template<int Size> +struct ei_decrement_if_fixed_size +{ + enum { + ret = (Size == Dynamic) ? Dynamic : Size-1 }; +}; + +#endif + +template< typename _Scalar, int _Deg > +class ei_companion +{ + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) + + enum { + Deg = _Deg, + Deg_1=ei_decrement_if_fixed_size<Deg>::ret + }; + + typedef _Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, Deg, 1> RightColumn; + //typedef DiagonalMatrix< Scalar, Deg_1, Deg_1 > BottomLeftDiagonal; + typedef Matrix<Scalar, Deg_1, 1> BottomLeftDiagonal; + + typedef Matrix<Scalar, Deg, Deg> DenseCompanionMatrixType; + typedef Matrix< Scalar, _Deg, Deg_1 > LeftBlock; + typedef Matrix< Scalar, Deg_1, Deg_1 > BottomLeftBlock; + typedef Matrix< Scalar, 1, Deg_1 > LeftBlockFirstRow; + + public: + EIGEN_STRONG_INLINE const _Scalar operator()( int row, int col ) const + { + if( m_bl_diag.rows() > col ) + { + if( 0 < row ){ return m_bl_diag[col]; } + else{ return 0; } + } + else{ return m_monic[row]; } + } + + public: + template<typename VectorType> + void setPolynomial( const VectorType& poly ) + { + const int deg = poly.size()-1; + m_monic = -1/poly[deg] * poly.head(deg); + //m_bl_diag.setIdentity( deg-1 ); + m_bl_diag.setOnes(deg-1); + } + + template<typename VectorType> + ei_companion( const VectorType& poly ){ + setPolynomial( poly ); } + + public: + DenseCompanionMatrixType denseMatrix() const + { + const int deg = m_monic.size(); + const int deg_1 = deg-1; + DenseCompanionMatrixType companion(deg,deg); + companion << + ( LeftBlock(deg,deg_1) + << LeftBlockFirstRow::Zero(1,deg_1), + BottomLeftBlock::Identity(deg-1,deg-1)*m_bl_diag.asDiagonal() ).finished() + , m_monic; + return companion; + } + + + + protected: + /** Helper function for the balancing algorithm. + * \returns true if the row and the column, having colNorm and rowNorm + * as norms, are balanced, false otherwise. + * colB and rowB are repectively the multipliers for + * the column and the row in order to balance them. + * */ + bool balanced( Scalar colNorm, Scalar rowNorm, + bool& isBalanced, Scalar& colB, Scalar& rowB ); + + /** Helper function for the balancing algorithm. + * \returns true if the row and the column, having colNorm and rowNorm + * as norms, are balanced, false otherwise. + * colB and rowB are repectively the multipliers for + * the column and the row in order to balance them. + * */ + bool balancedR( Scalar colNorm, Scalar rowNorm, + bool& isBalanced, Scalar& colB, Scalar& rowB ); + + public: + /** + * Balancing algorithm from B. N. PARLETT and C. REINSCH (1969) + * "Balancing a matrix for calculation of eigenvalues and eigenvectors" + * adapted to the case of companion matrices. + * A matrix with non zero row and non zero column is balanced + * for a certain norm if the i-th row and the i-th column + * have same norm for all i. + */ + void balance(); + + protected: + RightColumn m_monic; + BottomLeftDiagonal m_bl_diag; +}; + + + +template< typename _Scalar, int _Deg > +inline +bool ei_companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm, + bool& isBalanced, Scalar& colB, Scalar& rowB ) +{ + if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; } + else + { + //To find the balancing coefficients, if the radix is 2, + //one finds \f$ \sigma \f$ such that + // \f$ 2^{2\sigma-1} < rowNorm / colNorm \le 2^{2\sigma+1} \f$ + // then the balancing coefficient for the row is \f$ 1/2^{\sigma} \f$ + // and the balancing coefficient for the column is \f$ 2^{\sigma} \f$ + rowB = rowNorm / ei_radix<Scalar>(); + colB = Scalar(1); + const Scalar s = colNorm + rowNorm; + + while (colNorm < rowB) + { + colB *= ei_radix<Scalar>(); + colNorm *= ei_radix2<Scalar>(); + } + + rowB = rowNorm * ei_radix<Scalar>(); + + while (colNorm >= rowB) + { + colB /= ei_radix<Scalar>(); + colNorm /= ei_radix2<Scalar>(); + } + + //This line is used to avoid insubstantial balancing + if ((rowNorm + colNorm) < Scalar(0.95) * s * colB) + { + isBalanced = false; + rowB = Scalar(1) / colB; + return false; + } + else{ + return true; } + } +} + +template< typename _Scalar, int _Deg > +inline +bool ei_companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm, + bool& isBalanced, Scalar& colB, Scalar& rowB ) +{ + if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; } + else + { + /** + * Set the norm of the column and the row to the geometric mean + * of the row and column norm + */ + const _Scalar q = colNorm/rowNorm; + if( !ei_isApprox( q, _Scalar(1) ) ) + { + rowB = ei_sqrt( colNorm/rowNorm ); + colB = Scalar(1)/rowB; + + isBalanced = false; + return false; + } + else{ + return true; } + } +} + + +template< typename _Scalar, int _Deg > +void ei_companion<_Scalar,_Deg>::balance() +{ + EIGEN_STATIC_ASSERT( 1 < Deg, YOU_MADE_A_PROGRAMMING_MISTAKE ); + const int deg = m_monic.size(); + const int deg_1 = deg-1; + + bool hasConverged=false; + while( !hasConverged ) + { + hasConverged = true; + Scalar colNorm,rowNorm; + Scalar colB,rowB; + + //First row, first column excluding the diagonal + //============================================== + colNorm = ei_abs(m_bl_diag[0]); + rowNorm = ei_abs(m_monic[0]); + + //Compute balancing of the row and the column + if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) ) + { + m_bl_diag[0] *= colB; + m_monic[0] *= rowB; + } + + //Middle rows and columns excluding the diagonal + //============================================== + for( int i=1; i<deg_1; ++i ) + { + // column norm, excluding the diagonal + colNorm = ei_abs(m_bl_diag[i]); + + // row norm, excluding the diagonal + rowNorm = ei_abs(m_bl_diag[i-1]) + ei_abs(m_monic[i]); + + //Compute balancing of the row and the column + if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) ) + { + m_bl_diag[i] *= colB; + m_bl_diag[i-1] *= rowB; + m_monic[i] *= rowB; + } + } + + //Last row, last column excluding the diagonal + //============================================ + const int ebl = m_bl_diag.size()-1; + VectorBlock<RightColumn,Deg_1> headMonic( m_monic, 0, deg_1 ); + colNorm = headMonic.array().abs().sum(); + rowNorm = ei_abs( m_bl_diag[ebl] ); + + //Compute balancing of the row and the column + if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) ) + { + headMonic *= colB; + m_bl_diag[ebl] *= rowB; + } + } +} + + +#endif // EIGEN_COMPANION_H diff --git a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h new file mode 100644 index 000000000..49a5ffffa --- /dev/null +++ b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h @@ -0,0 +1,395 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_POLYNOMIAL_SOLVER_H +#define EIGEN_POLYNOMIAL_SOLVER_H + +/** \ingroup Polynomials_Module + * \class PolynomialSolverBase. + * + * \brief Defined to be inherited by polynomial solvers: it provides + * convenient methods such as + * - real roots, + * - greatest, smallest complex roots, + * - real roots with greatest, smallest absolute real value, + * - greatest, smallest real roots. + * + * It stores the set of roots as a vector of complexes. + * + */ +template< typename _Scalar, int _Deg > +class PolynomialSolverBase +{ + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) + + typedef _Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef std::complex<RealScalar> RootType; + typedef Matrix<RootType,_Deg,1> RootsType; + + protected: + template< typename OtherPolynomial > + inline void setPolynomial( const OtherPolynomial& poly ){ + m_roots.resize(poly.size()); } + + public: + template< typename OtherPolynomial > + inline PolynomialSolverBase( const OtherPolynomial& poly ){ + setPolynomial( poly() ); } + + inline PolynomialSolverBase(){} + + public: + /** \returns the complex roots of the polynomial */ + inline const RootsType& roots() const { return m_roots; } + + public: + /** Clear and fills the back insertion sequence with the real roots of the polynomial + * i.e. the real part of the complex roots that have an imaginary part which + * absolute value is smaller than absImaginaryThreshold. + * absImaginaryThreshold takes the dummy_precision associated + * with the _Scalar template parameter of the PolynomialSolver class as the default value. + * + * \param[out] bi_seq : the back insertion sequence (stl concept) + * \param[in] absImaginaryThreshold : the maximum bound of the imaginary part of a complex + * number that is considered as real. + * */ + template<typename Stl_back_insertion_sequence> + inline void realRoots( Stl_back_insertion_sequence& bi_seq, + const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const + { + bi_seq.clear(); + for( int i=0; i<m_roots.size(); ++i ) + { + if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold ){ + bi_seq.push_back( m_roots[i].real() ); } + } + } + + protected: + template<typename squaredNormBinaryPredicate> + inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const + { + int res=0; + RealScalar norm2 = ei_abs2( m_roots[0] ); + for( int i=1; i<m_roots.size(); ++i ) + { + const RealScalar currNorm2 = ei_abs2( m_roots[i] ); + if( pred( currNorm2, norm2 ) ){ + res=i; norm2=currNorm2; } + } + return m_roots[res]; + } + + public: + /** + * \returns the complex root with greatest norm. + */ + inline const RootType& greatestRoot() const + { + std::greater<Scalar> greater; + return selectComplexRoot_withRespectToNorm( greater ); + } + + /** + * \returns the complex root with smallest norm. + */ + inline const RootType& smallestRoot() const + { + std::less<Scalar> less; + return selectComplexRoot_withRespectToNorm( less ); + } + + protected: + template<typename squaredRealPartBinaryPredicate> + inline const RealScalar& selectRealRoot_withRespectToAbsRealPart( + squaredRealPartBinaryPredicate& pred, + bool& hasArealRoot, + const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const + { + hasArealRoot = false; + int res=0; + RealScalar abs2; + + for( int i=0; i<m_roots.size(); ++i ) + { + if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold ) + { + if( !hasArealRoot ) + { + hasArealRoot = true; + res = i; + abs2 = m_roots[i].real() * m_roots[i].real(); + } + else + { + const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real(); + if( pred( currAbs2, abs2 ) ) + { + abs2 = currAbs2; + res = i; + } + } + } + else + { + if( ei_abs( m_roots[i].imag() ) < ei_abs( m_roots[res].imag() ) ){ + res = i; } + } + } + return m_roots[res].real(); + } + + + template<typename RealPartBinaryPredicate> + inline const RealScalar& selectRealRoot_withRespectToRealPart( + RealPartBinaryPredicate& pred, + bool& hasArealRoot, + const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const + { + hasArealRoot = false; + int res=0; + RealScalar val; + + for( int i=0; i<m_roots.size(); ++i ) + { + if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold ) + { + if( !hasArealRoot ) + { + hasArealRoot = true; + res = i; + val = m_roots[i].real(); + } + else + { + const RealScalar curr = m_roots[i].real(); + if( pred( curr, val ) ) + { + val = curr; + res = i; + } + } + } + else + { + if( ei_abs( m_roots[i].imag() ) < ei_abs( m_roots[res].imag() ) ){ + res = i; } + } + } + return m_roots[res].real(); + } + + public: + /** + * \returns a real root with greatest absolute magnitude. + * A real root is defined as the real part of a complex root with absolute imaginary + * part smallest than absImaginaryThreshold. + * absImaginaryThreshold takes the dummy_precision associated + * with the _Scalar template parameter of the PolynomialSolver class as the default value. + * If no real root is found the boolean hasArealRoot is set to false and the real part of + * the root with smallest absolute imaginary part is returned instead. + * + * \param[out] hasArealRoot : boolean true if a real root is found according to the + * absImaginaryThreshold criterion, false otherwise. + * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide + * whether or not a root is real. + */ + inline const RealScalar& absGreatestRealRoot( + bool& hasArealRoot, + const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const + { + std::greater<Scalar> greater; + return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold ); + } + + + /** + * \returns a real root with smallest absolute magnitude. + * A real root is defined as the real part of a complex root with absolute imaginary + * part smallest than absImaginaryThreshold. + * absImaginaryThreshold takes the dummy_precision associated + * with the _Scalar template parameter of the PolynomialSolver class as the default value. + * If no real root is found the boolean hasArealRoot is set to false and the real part of + * the root with smallest absolute imaginary part is returned instead. + * + * \param[out] hasArealRoot : boolean true if a real root is found according to the + * absImaginaryThreshold criterion, false otherwise. + * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide + * whether or not a root is real. + */ + inline const RealScalar& absSmallestRealRoot( + bool& hasArealRoot, + const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const + { + std::less<Scalar> less; + return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold ); + } + + + /** + * \returns the real root with greatest value. + * A real root is defined as the real part of a complex root with absolute imaginary + * part smallest than absImaginaryThreshold. + * absImaginaryThreshold takes the dummy_precision associated + * with the _Scalar template parameter of the PolynomialSolver class as the default value. + * If no real root is found the boolean hasArealRoot is set to false and the real part of + * the root with smallest absolute imaginary part is returned instead. + * + * \param[out] hasArealRoot : boolean true if a real root is found according to the + * absImaginaryThreshold criterion, false otherwise. + * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide + * whether or not a root is real. + */ + inline const RealScalar& greatestRealRoot( + bool& hasArealRoot, + const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const + { + std::greater<Scalar> greater; + return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold ); + } + + + /** + * \returns the real root with smallest value. + * A real root is defined as the real part of a complex root with absolute imaginary + * part smallest than absImaginaryThreshold. + * absImaginaryThreshold takes the dummy_precision associated + * with the _Scalar template parameter of the PolynomialSolver class as the default value. + * If no real root is found the boolean hasArealRoot is set to false and the real part of + * the root with smallest absolute imaginary part is returned instead. + * + * \param[out] hasArealRoot : boolean true if a real root is found according to the + * absImaginaryThreshold criterion, false otherwise. + * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide + * whether or not a root is real. + */ + inline const RealScalar& smallestRealRoot( + bool& hasArealRoot, + const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const + { + std::less<Scalar> less; + return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold ); + } + + protected: + RootsType m_roots; +}; + +#define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \ + typedef typename BASE::Scalar Scalar; \ + typedef typename BASE::RealScalar RealScalar; \ + typedef typename BASE::RootType RootType; \ + typedef typename BASE::RootsType RootsType; + + + +/** \ingroup Polynomials_Module + * + * \class PolynomialSolver + * + * \brief A polynomial solver + * + * Computes the complex roots of a real polynomial. + * + * \param _Scalar the scalar type, i.e., the type of the polynomial coefficients + * \param _Deg the degree of the polynomial, can be a compile time value or Dynamic. + * Notice that the number of polynomial coefficients is _Deg+1. + * + * This class implements a polynomial solver and provides convenient methods such as + * - real roots, + * - greatest, smallest complex roots, + * - real roots with greatest, smallest absolute real value. + * - greatest, smallest real roots. + * + * WARNING: this polynomial solver is experimental, part of the unsuported Eigen modules. + * + * + * Currently a QR algorithm is used to compute the eigenvalues of the companion matrix of + * the polynomial to compute its roots. + * This supposes that the complex moduli of the roots are all distinct: e.g. there should + * be no multiple roots or conjugate roots for instance. + * With 32bit (float) floating types this problem shows up frequently. + * However, almost always, correct accuracy is reached even in these cases for 64bit + * (double) floating types and small polynomial degree (<20). + */ +template< typename _Scalar, int _Deg > +class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> +{ + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) + + typedef PolynomialSolverBase<_Scalar,_Deg> PS_Base; + EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) + + typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType; + typedef EigenSolver<CompanionMatrixType> EigenSolverType; + + public: + /** Computes the complex roots of a new polynomial. */ + template< typename OtherPolynomial > + void compute( const OtherPolynomial& poly ) + { + assert( Scalar(0) != poly[poly.size()-1] ); + ei_companion<Scalar,_Deg> companion( poly ); + companion.balance(); + m_eigenSolver.compute( companion.denseMatrix() ); + m_roots = m_eigenSolver.eigenvalues(); + } + + public: + template< typename OtherPolynomial > + inline PolynomialSolver( const OtherPolynomial& poly ){ + compute( poly ); } + + inline PolynomialSolver(){} + + protected: + using PS_Base::m_roots; + EigenSolverType m_eigenSolver; +}; + + +template< typename _Scalar > +class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1> +{ + public: + typedef PolynomialSolverBase<_Scalar,1> PS_Base; + EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) + + public: + /** Computes the complex roots of a new polynomial. */ + template< typename OtherPolynomial > + void compute( const OtherPolynomial& poly ) + { + assert( Scalar(0) != poly[poly.size()-1] ); + m_roots[0] = -poly[0]/poly[poly.size()-1]; + } + + protected: + using PS_Base::m_roots; +}; + +#endif // EIGEN_POLYNOMIAL_SOLVER_H diff --git a/unsupported/Eigen/src/Polynomials/PolynomialUtils.h b/unsupported/Eigen/src/Polynomials/PolynomialUtils.h new file mode 100644 index 000000000..d78821f90 --- /dev/null +++ b/unsupported/Eigen/src/Polynomials/PolynomialUtils.h @@ -0,0 +1,153 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_POLYNOMIAL_UTILS_H +#define EIGEN_POLYNOMIAL_UTILS_H + +/** \ingroup Polynomials_Module + * \returns the evaluation of the polynomial at x using Horner algorithm. + * + * \param[in] poly : the vector of coefficients of the polynomial ordered + * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial + * e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$. + * \param[in] x : the value to evaluate the polynomial at. + * + * <i><b>Note for stability:</b></i> + * <dd> \f$ |x| \le 1 \f$ </dd> + */ +template <typename Polynomials, typename T> +inline +T poly_eval_horner( const Polynomials& poly, const T& x ) +{ + T val=poly[poly.size()-1]; + for( int i=poly.size()-2; i>=0; --i ){ + val = val*x + poly[i]; } + return val; +} + +/** \ingroup Polynomials_Module + * \returns the evaluation of the polynomial at x using stabilized Horner algorithm. + * + * \param[in] poly : the vector of coefficients of the polynomial ordered + * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial + * e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$. + * \param[in] x : the value to evaluate the polynomial at. + */ +template <typename Polynomials, typename T> +inline +T poly_eval( const Polynomials& poly, const T& x ) +{ + typedef typename NumTraits<T>::Real Real; + + if( ei_abs2( x ) <= Real(1) ){ + return poly_eval_horner( poly, x ); } + else + { + T val=poly[0]; + T inv_x = T(1)/x; + for( int i=1; i<poly.size(); ++i ){ + val = val*inv_x + poly[i]; } + + return std::pow(x,(T)(poly.size()-1)) * val; + } +} + +/** \ingroup Polynomials_Module + * \returns a maximum bound for the absolute value of any root of the polynomial. + * + * \param[in] poly : the vector of coefficients of the polynomial ordered + * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial + * e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$. + * + * <i><b>Precondition:</b></i> + * <dd> the leading coefficient of the input polynomial poly must be non zero </dd> + */ +template <typename Polynomial> +inline +typename NumTraits<typename Polynomial::Scalar>::Real cauchy_max_bound( const Polynomial& poly ) +{ + typedef typename Polynomial::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real Real; + + assert( Scalar(0) != poly[poly.size()-1] ); + const Scalar inv_leading_coeff = Scalar(1)/poly[poly.size()-1]; + Real cb(0); + + for( int i=0; i<poly.size()-1; ++i ){ + cb += ei_abs(poly[i]*inv_leading_coeff); } + return cb + Real(1); +} + +/** \ingroup Polynomials_Module + * \returns a minimum bound for the absolute value of any non zero root of the polynomial. + * \param[in] poly : the vector of coefficients of the polynomial ordered + * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial + * e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$. + */ +template <typename Polynomial> +inline +typename NumTraits<typename Polynomial::Scalar>::Real cauchy_min_bound( const Polynomial& poly ) +{ + typedef typename Polynomial::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real Real; + + int i=0; + while( i<poly.size()-1 && Scalar(0) == poly(i) ){ ++i; } + if( poly.size()-1 == i ){ + return Real(1); } + + const Scalar inv_min_coeff = Scalar(1)/poly[i]; + Real cb(1); + for( int j=i+1; j<poly.size(); ++j ){ + cb += ei_abs(poly[j]*inv_min_coeff); } + return Real(1)/cb; +} + +/** \ingroup Polynomials_Module + * Given the roots of a polynomial compute the coefficients in the + * monomial basis of the monic polynomial with same roots and minimal degree. + * If RootVector is a vector of complexes, Polynomial should also be a vector + * of complexes. + * \param[in] rv : a vector containing the roots of a polynomial. + * \param[out] poly : the vector of coefficients of the polynomial ordered + * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial + * e.g. \f$ 3 + x^2 \f$ is stored as a vector \f$ [ 3, 0, 1 ] \f$. + */ +template <typename RootVector, typename Polynomial> +void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly ) +{ + + typedef typename Polynomial::Scalar Scalar; + + poly.setZero( rv.size()+1 ); + poly[0] = -rv[0]; poly[1] = Scalar(1); + for( int i=1; i<(int)rv.size(); ++i ) + { + for( int j=i+1; j>0; --j ){ poly[j] = poly[j-1] - rv[i]*poly[j]; } + poly[0] = -rv[i]*poly[0]; + } +} + + +#endif // EIGEN_POLYNOMIAL_UTILS_H |