Hi Gael, You finally modified PolynomialSolver.h code for finding roots of polynomial with complex coefficients and wrote unit test. This enhancement is very important for signal processing application. This modification of PolynomialSolver.h and Companion.h are rather trivial and localized. So it is safe to modify existing code. Would it be possible to include them to the next release?
Best regards, Oleg On Wed, Nov 23, 2016 at 5:09 PM Gael Guennebaud <[email protected]> wrote: > Done: > > https://bitbucket.org/eigen/eigen/commits/3842186ae759/ > https://bitbucket.org/eigen/eigen/commits/eff1bfef71a2/ > https://bitbucket.org/eigen/eigen/commits/df6c723906fc/ > > gael > > On Wed, Nov 23, 2016 at 3:39 PM, Gael Guennebaud < > [email protected]> wrote: > >> Sorry for late reply, but thanks for the patches. >> >> Regarding the eigen-solver choice issue, it's possible to automatically >> use the right one based on NumTraits<Scalar>::IsComplex and conditional<>. >> >> We would also need to extended the unit test to cover complexes. >> >> gael >> >> On Tue, Nov 8, 2016 at 3:41 PM, Oleg Shirokobrod < >> [email protected]> wrote: >> >>> Dear list, >>> >>> Polynomial solver currently finds roots only for real polynomial >>> coefficients. Having Eigen::ComplexEigenSolver lets Polynomial module to >>> compute roots for complex polynomial coefficients. I have made necessary >>> modifications. I have attached corresponding patch files. Unfortunately >>> module cannot deduce type of the solver from partial specialization of >>> corresponding matrix type, such that >>> >>> template<typename RealScalar> class EigenSolver {}; >>> >>> template<typenameRealScalar> class >>> EigenSolver<Eigen::Matrix<RealScalar,_Deg,_Deg> > {// EigenSolver }; >>> >>> template<typename RealScalar> class >>> EigenSolver<Eigen::Matrix<std::complex<RealScalar>,_Deg,_Deg > > { >>> //ComplexEigenSolver}; >>> >>> So I have to add additional template parameter EigenSolverType with >>> default value for real coefficients: >>> >>> template< typename _Scalar, int _Deg , typename EigenSolverType = >>> EigenSolver<Matrix<_Scalar,_Deg,_Deg> > > >>> class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>{}; >>> >>> I have to replace in number of places in the file Companion.h (mainly in >>> functions balanced() and balance()) Scalar with RealScalar, where variables >>> are really real. >>> >>> With this code I have run test against Matlab and it gave similar >>> results. >>> >>> Test >>> VectorXcd roots = VectorXcd::Random(4); >>> VectorXcd polynomialCoefficients; >>> roots_to_monicPolynomial(roots, polynomialCoefficients); >>> >>> cout << "roots : " << endl; >>> cout << setprecision(14) << roots << endl; >>> cout << "polynomialCoefficients : " << endl; >>> cout << setprecision(14) << polynomialCoefficients << endl; >>> >>> PolynomialSolver<std::complex<double>, Dynamic, >>> ComplexEigenSolver<MatrixXcd> > psolve(polynomialCoefficients); >>> >>> VectorXcd computedRoots = psolve.roots(); >>> >>> cout << "computedRoots : " << endl; >>> cout << setprecision(14) << computedRoots << endl; >>> >>> for(int index = 0; index < computedRoots.size(); ++index) >>> { >>> cout << setprecision(14) << "polynom at computedRoots * 1.0 e-11: " >>> << poly_eval(polynomialCoefficients, computedRoots(index))*1.0e11 << endl; >>> } >>> >>> Output: >>> roots : >>> (0.12717062898648,-0.99749748222297) >>> (0.61748100222785,-0.61339152195807) >>> (-0.04025391399884,0.17001861629078) >>> (0.79192480239265,-0.29941709646901) >>> >>> polynomialCoefficients : >>> (0.091649023983597,-0.091441416775918) >>> (0.24020598856634,0.37401934653925) >>> (-0.16301627948124,-1.8544616197629) >>> (-1.4963225196081,1.7402874843593) >>> (1,0) >>> >>> computedRoots : >>> (-0.04025391399884,0.17001861629078) >>> (0.79192480239265,-0.29941709646901) >>> (0.61748100222785,-0.61339152195807) >>> (0.12717062898648,-0.99749748222297) >>> >>> polynom at computedRoots * 1.0 e-11: >>> (8.3266726846887e-006,-1.2490009027033e-005) >>> polynom at computedRoots * 1.0 e-11: >>> (4.7184478546569e-005,-1.6653345369377e-005) >>> polynom at computedRoots * 1.0 e-11: >>> (2.2204460492503e-005,-1.3877787807814e-005) >>> polynom at computedRoots * 1.0 e-11: >>> (1.5163809063462e-005,-2.7286889655471e-005) >>> >>> Matlab for the same coefficients gives following roots (remember that >>> matlab array of coefficients and Eigen one are reversed to each other) >>> >>> computedRoots : >>> 0.127170628986480 - 0.997497482222969i >>> 0.617481002227849 - 0.613391521958071i >>> 0.791924802392650 - 0.299417096469009i >>> -0.040253913998840 + 0.170018616290780i >>> >>> Best regards, >>> >>> Oleg Shirokobrod >>> >>> >>> >> >
