diff --git a/MMVII/ExternalInclude/eigen-3.4.0/unsupported/Eigen/src/Polynomials/PolynomialUtils.h b/MMVII/ExternalInclude/eigen-3.4.0/unsupported/Eigen/src/Polynomials/PolynomialUtils.h index 394e857acf..b490a03588 100644 --- a/MMVII/ExternalInclude/eigen-3.4.0/unsupported/Eigen/src/Polynomials/PolynomialUtils.h +++ b/MMVII/ExternalInclude/eigen-3.4.0/unsupported/Eigen/src/Polynomials/PolynomialUtils.h @@ -53,7 +53,8 @@ T poly_eval( const Polynomials& poly, const T& x ) { T val=poly[0]; T inv_x = T(1)/x; - for( DenseIndex i=1; i class cPolynom Type Value(const Type & aVal) const; + /// return som(|a_k x^k|) , used for some bounding stuffs + Type AbsValue(const Type & aVal) const; cPolynom operator * (const cPolynom & aP2) const; cPolynom operator + (const cPolynom & aP2) const; cPolynom operator - (const cPolynom & aP2) const; cPolynom operator * (const Type & aVal) const; - std::vector RealRoots(const Type & aTol,int ItMax); + cPolynom Deriv() const; + + std::vector RealRoots(const Type & aTol,int ItMax) const; Type& operator [] (size_t aK) {return mVCoeffs[aK];} diff --git a/MMVII/src/MMV1/Numerics.cpp b/MMVII/src/MMV1/Numerics.cpp index 072cfdf883..9546bf2a0f 100755 --- a/MMVII/src/MMV1/Numerics.cpp +++ b/MMVII/src/MMV1/Numerics.cpp @@ -54,16 +54,6 @@ INST_V1ROOTS(tREAL8) INST_V1ROOTS(tREAL16) -/* -template void - RealRootsOfRealPolynome - ( - ElSTDNS vector & Sols, - const ElPolynome &aPol, - Type tol, - INT ItMax - ) - */ }; diff --git a/MMVII/src/UtiMaths/Polynoms.cpp b/MMVII/src/UtiMaths/Polynoms.cpp index 3cd3ddd8e3..3064ccc416 100755 --- a/MMVII/src/UtiMaths/Polynoms.cpp +++ b/MMVII/src/UtiMaths/Polynoms.cpp @@ -1,17 +1,303 @@ + +#if defined(__GNUC__) && !defined(__clang__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wmaybe-uninitialized" +#endif + #include "cMMVII_Appli.h" #include "V1VII.h" +#include +// #include "unsupported/Eigen/Polynomials" + +//#if defined(__GNUC__) && !defined(__clang__) +// # pragma GCC diagnostic pop +// #endif -/* // #include "include/MMVII_nums.h" // #include "include/MMVII_Bench.h" - //#include "Eigen/unsupported/Eigen/Polynomials" - */ +//#include +//#include + +using namespace Eigen; +using namespace std; namespace MMVII { +/** + Class for extraction of roots of polynoms using "Companion matrix method" + + See https://en.wikipedia.org/wiki/Companion_matrix +*/ + + +template class cEigenPolynRoots +{ + public : + cEigenPolynRoots(const cPolynom &) ; + const std::vector & RealRoots() const; + const VectorXcd & ComplexRoots() const; + bool RootIsReal(const std::complex &,std::string * sayWhy=nullptr); + const MatrixXd & CompM() const; + private : + + cPolynom mPol; + cPolynom mDPol; + size_t mDeg; + MatrixXd mCompM; ///< companion matrix + VectorXcd mCEV; ///< Complex eigen values + std::vector mRR; ///< Real roots + +}; + + +template cEigenPolynRoots::cEigenPolynRoots(const cPolynom & aPol) : + mPol (aPol), + mDPol (aPol.Deriv()), + mDeg (mPol.Degree()), + mCompM (mDeg,mDeg) +{ + + if (mDeg==0) + return; + + for (size_t aI = 0; aI < mDeg; ++aI) + for (size_t aJ = 0; aJ < mDeg; ++aJ) + mCompM(aI,aJ) = 0; + + // fill the diagonal up the principal diag + for (size_t aK = 0; aK < mDeg-1; ++aK) + { + mCompM(aK+1, aK ) = 1; + } + + const Type & aHighCoeff = mPol[mDeg]; + // Fill last line with normalized coeff + for (size_t aK = 0; aK < mDeg; ++aK) + { + mCompM(aK, mDeg - 1) = -mPol[aK] / aHighCoeff; + } + + EigenSolver aSolver(mCompM); + mCEV = aSolver.eigenvalues() ; + + for (const auto & aC : mCEV) + if (RootIsReal(aC)) + mRR.push_back(aC.real()); + std::sort(mRR.begin(),mRR.end()); +} + +template const std::vector & cEigenPolynRoots::RealRoots() const {return mRR;} + +template const VectorXcd & cEigenPolynRoots::ComplexRoots() const {return mCEV;} + +template const MatrixXd & cEigenPolynRoots::CompM() const {return mCompM;} + + +/** Also the question seems pretty basic, it becomes more complicated due to numericall approximation */ + +template bool cEigenPolynRoots::RootIsReal(const std::complex & aC,std::string * sayWhy) +{ + // [1] Test is "aC" is a real number + + tREAL8 C_i = aC.imag(); + // [1.1] if absolute value of imaginary part is "big" it's not + if (std::abs(C_i) > 1e-5) + { + if (sayWhy) + *sayWhy = "ABS REAL COMPLEX=" + ToStr(std::abs(C_i)); + return false; + } + + tREAL8 C_r = aC.real(); + // [1.1] if relative imaginary part is "big" + if (std::abs(C_i) > 1e-8 * (std::abs(C_r)+1e-5)) + { + if (sayWhy) + *sayWhy = "RELAT REAL COMPLEX=" + ToStr(std::abs(C_i)/(std::abs(C_r)+1e-5)); + return false; + } + + // [2] Test + tREAL8 aAbsVP = std::abs(mPol.Value(C_r)); + // [2.1] if absolute value of polynom is big + if (aAbsVP > 1e-5) + { + if (sayWhy) + *sayWhy = "ABS VALUE POL " + ToStr(aAbsVP); + return false; + } + + tREAL8 aAVA = mPol.AbsValue(C_r); + // [2.1] if absolute value of polynom is big relatively to norm + if (aAbsVP > 1e-8 * (aAVA+1e-5)) + { + if (sayWhy) + *sayWhy = "RELATIVE VALUE POL " + ToStr(aAbsVP/(aAVA+1e-5)); + return false; + } + + if (sayWhy) + *sayWhy = "is real"; + + return true; +} + + +template class cEigenPolynRoots; + +template void My_Roots(const cPolynom & aPol1) +{ +} +template <> void My_Roots(const cPolynom & aPol1) +{ +// StdOut() << "DDDD " << aPol1.Degree() << "\n"; + // (X2+1)(X-1) = X3-X2+X-1 + int aNb=300; + // vector aCoeffs1 = {-1,1,-1,1,5,-2,0.12}; + // cPolynom aPol1(aCoeffs1); + + cAutoTimerSegm aTimeEigen(GlobAppTS(),"Eigen"); + for (int aK=0 ; aK aEPR(aPol1); + } + + cAutoTimerSegm aTimeV1(GlobAppTS(),"V1"); + for (int aK=0 ; aK aEPR(aPol1); + auto aV2 = aEPR.RealRoots(); + auto aV1 = aPol1.RealRoots(1e-20,60); + StdOut() << " SZzzzZ= " << aV1.size() << " " << aV2.size() << "\n"; + if (aV1.size() != aV2.size()) + { + StdOut() << "Coeffs=" << aPol1.VCoeffs() << "\n"; + StdOut() << "V1=" << aV1 << "\n"; + StdOut() << "V2=" << aV2 << "\n"; + for (const auto & aC : aEPR.ComplexRoots()) + { + std::string strWhy; + bool isR = aEPR.RootIsReal(aC,&strWhy); + StdOut() << "R=" << isR << " C=" << aC << " W=" << strWhy << "\n"; + } + + StdOut() << "DET=" << aEPR.CompM().determinant() << "\n"; + StdOut() << " ------------------ MAT ---------------------\n"; + StdOut() << aEPR.CompM() << "\n"; +getchar(); + } + else + StdOut() << aV1 << aV2 << "\n"; + // (X2+1)(X-1) = X3-X2+X-1 + // vector coeffs = {-1,1,-1,1}; +} + + +#if (0) + +/* +MatrixXd createCompanionMatrix(const vector& coeffs) { +StdOut() << "BEGIN createCompanionMatrixcreateCompanionMatrixcreateCompanionMatrix\n"; + int n = coeffs.size(); + MatrixXd companionMatrix(n - 1, n - 1); + +StdOut() << "LLLLL " << __LINE__ << "\n"; + + // Remplir la matrice compagnon + // for (int i = 0; i < n - 1; ++i) { + for (int i = 0; i < n - 2; ++i) { + companionMatrix(i+1, i ) = 1; // Remplir la diagonale au-dessus de la diagonale principale + } +StdOut() << "LLLLL " << __LINE__ << "\n"; + + // Remplir la dernière colonne de la matrice avec les coefficients du polynôme + for (int i = 0; i < n - 1; ++i) { + // companionMatrix(i, n - 1) = -coeffs[i] / coeffs[n - 1]; + companionMatrix(i, n - 2) = -coeffs[i] / coeffs[n - 1]; + } + +StdOut() << "END createCompanionMatrixcreateCompanionMatrixcreateCompanionMatrix\n"; + return companionMatrix; +} + +// Fonction principale +int My_Roots() { + // (X2+1)(X-1) = X3-X2+X-1 + // Exemple de coefficients d'un polynôme : x^3 - 6x^2 + 11x - 6 + // (x-1) (x2+ 5x +6) + // vector coeffs = {1, -6, 11, -6}; + vector coeffs = {1,-1,1,-1}; + + + + // Créer la matrice compagnon + MatrixXd companionMatrix = createCompanionMatrix(coeffs); + + std::cout << "mmmm " << companionMatrix << "\n"; + + // Calculer les valeurs propres (racines du polynôme) + EigenSolver solver(companionMatrix); + std::cout << "eeeee " << solver.eigenvalues() << "\n"; + + // const auto& theEigenV = solver.eigenvalues(); + VectorXcd eivals = solver.eigenvalues(); + std::cout << "EEEEEE " << eivals << "\n"; + + VectorXd R_roots = solver.eigenvalues().real(); + VectorXd I_roots = solver.eigenvalues().imag(); + + for (int i = 0; i < R_roots.size(); ++i) { + cout << "R=" << R_roots[i] << " I=" << I_roots[i] << endl; + } +*/ + +/* + VectorXd roots = solver.eigenvalues().real(); + + // Afficher les racines + cout << "Les racines du polynôme sont : " << endl; + for (int i = 0; i < roots.size(); ++i) { + cout << roots[i] << endl; + } +*/ + + return 0; +} +#endif + +/* + */ + + + + +#if (0) +void TestPolynEigen() +{ + // cEigenMMVIIPoly aPoly({-1,0,3}); + std::vector aPoly {-1,0,3}; + // bool hasArealRoot; + // tREAL8 aS = greatestRealRoot(hasArealRoot); + + StdOut() << "EEE:V= " << Eigen::poly_eval (aPoly,2) << "\n"; + + // polynomialsolver( internal::random(9,13)); + PolynomialSolver aSolver; + + aSolver.compute(aPoly); + // aSolver.roots(); + // FakeUseIt(aSolver); + +} +#endif + /* ************************************************************************ */ /* */ @@ -24,7 +310,7 @@ namespace MMVII template cPolynom::cPolynom(const tCoeffs & aVCoeffs) : mVCoeffs (aVCoeffs) { - MMVII_INTERNAL_ASSERT_tiny(!mVCoeffs.empty(),"Empty polynom not handled"); + // MMVII_INTERNAL_ASSERT_tiny(!mVCoeffs.empty(),"Empty polynom not handled"); } template cPolynom::cPolynom(const cPolynom & aPol) : cPolynom(aPol.mVCoeffs) @@ -96,8 +382,9 @@ template return aRes; } -template std::vector cPolynom::RealRoots(const Type & aTol,int ItMax) +template std::vector cPolynom::RealRoots(const Type & aTol,int ItMax) const { +// StdOut() << "RealRootsRealRootsRealRootsRealRootsRealRootsRealRootsRealRootsRealRoots \n"; getchar(); return V1RealRoots(mVCoeffs,aTol,ItMax); } @@ -119,6 +406,22 @@ template Type cPolynom::Value(const Type & aVal) const return aResult; } +template Type cPolynom::AbsValue(const Type & aVal) const +{ + Type aResult = 0.0; + Type aPowV = 1.0; + for (const auto & aCoef : mVCoeffs) + { + aResult += std::abs(aCoef * aPowV); + aPowV *= aVal; + } + return aResult; +} + + + + + template cPolynom cPolynom::operator * (const cPolynom & aP2) const { cPolynom aRes(Degree() + aP2.Degree()); @@ -159,6 +462,15 @@ template cPolynom cPolynom::operator - (const cPolynom return aRes; } +template cPolynom cPolynom::Deriv() const +{ + std::vector aVCD; // Vector Coeff Derivates + for (size_t aDeg=1 ; aDeg(aVCD); +} + template cPolynom cPolynom::operator * (const Type & aMul) const { @@ -201,6 +513,10 @@ template void TplBenchPolynome() cPolynom aP1min2 = aPol1 - aPol2; cPolynom aP2min1 = aPol2 - aPol1; + cPolynom aDerP1P2_A = aP1mul2.Deriv(); + cPolynom aDerP1P2_B = aPol1.Deriv() * aPol2 + aPol1*aPol2.Deriv(); + + Type aEps = tElemNumTrait ::Accuracy(); for (int aK=0 ; aK< 20 ; aK++) @@ -210,6 +526,10 @@ template void TplBenchPolynome() Type aChekP = aPol1.Value(aV) + aPol2.Value(aV); Type aChekMin = aPol1.Value(aV) - aPol2.Value(aV); + Type aDerA = aDerP1P2_A.Value(aV) ; + Type aDerB = aDerP1P2_B.Value(aV) ; + MMVII_INTERNAL_ASSERT_bench(std::abs(aDerA-aDerB) void TplBenchPolynome() MMVII_INTERNAL_ASSERT_bench(std::abs(RelativeSafeDifference(-aChekMin,aP2min1.Value(aV))) aVRootsGen; Type aAmpl = 10*RandUnif_NotNull(1e-2); @@ -247,6 +567,7 @@ template void TplBenchPolynome() Type aDif = RelativeSafeDifference(aVRootsGen[aK],aVRootsCalc[aK]); MMVII_INTERNAL_ASSERT_bench(aDif(); TplBenchPolynome(); TplBenchPolynome();