From b6af673186ad7db517b99f8bb0ef2a60b492739c Mon Sep 17 00:00:00 2001 From: Dario Sassi Date: Sun, 12 Jan 2014 22:50:16 +0000 Subject: [PATCH] EgtNumKernel 1.5a2 : Aggiunta classe Polynomial. Corretto calcolo zeri di polinomi reali. --- EgtNumKernel.rc | Bin 7626 -> 7626 bytes EgtNumKernel.vcxproj | 11 +- EgtNumKernel.vcxproj.filters | 18 +- JenkinsTraub.cpp | 31 ++- JenkinsTraub.h | 2 +- Polynomial.cpp | 219 +++++++++++++++++++++ PolynomialZeros.cpp => PolynomialRoots.cpp | 55 +++--- 7 files changed, 295 insertions(+), 41 deletions(-) create mode 100644 Polynomial.cpp rename PolynomialZeros.cpp => PolynomialRoots.cpp (76%) diff --git a/EgtNumKernel.rc b/EgtNumKernel.rc index a8d3075e67da1772ba8d3c109eb414985516d182..739fc7040a177f9a6dbe640fb8e70db619c22c95 100644 GIT binary patch delta 55 zcmX?Qead>nH#SD2&8PWJGfn2=(V5)9$uapGryQfv{p&8r& DH(e5j delta 55 zcmX?Qead>nH#SDY&8PWJGfn2=(V5)9$uapGryQf<{p&8r& DHun;O diff --git a/EgtNumKernel.vcxproj b/EgtNumKernel.vcxproj index c844e3f..39488c2 100644 --- a/EgtNumKernel.vcxproj +++ b/EgtNumKernel.vcxproj @@ -87,8 +87,7 @@ copy $(TargetPath) \EgtProg\Dll true - copy $(TargetDir)$(TargetName).pdb \EgtDev\Lib\ -copy $(TargetDir)$(TargetName).lib \EgtDev\Lib\ + copy $(TargetDir)$(TargetName).lib \EgtDev\Lib\ copy $(TargetPath) \EgtProg\Dll @@ -99,7 +98,8 @@ copy $(TargetPath) \EgtProg\Dll - + + Create Create @@ -109,8 +109,11 @@ copy $(TargetPath) \EgtProg\Dll + - + + + diff --git a/EgtNumKernel.vcxproj.filters b/EgtNumKernel.vcxproj.filters index d516aaa..1e50e26 100644 --- a/EgtNumKernel.vcxproj.filters +++ b/EgtNumKernel.vcxproj.filters @@ -27,7 +27,10 @@ File di origine - + + File di origine + + File di origine @@ -53,10 +56,19 @@ File di intestazione - + File di intestazione - + + File di intestazione + + + File di intestazione + + + File di intestazione + + File di intestazione diff --git a/JenkinsTraub.cpp b/JenkinsTraub.cpp index 32c95b8..2694ea8 100644 --- a/JenkinsTraub.cpp +++ b/JenkinsTraub.cpp @@ -16,6 +16,12 @@ #include "stdafx.h" #include "JenkinsTraub.h" +//--------------------------- Constants -------------------------------------- +#define RADIX_DOUBLE 2 +#define EPS_DOUBLE 2.22e-16 +#define SMALL_DOUBLE 2.3e-308 +#define BIG_DOUBLE 1.7e+308 + //--------------------------- Class Rpoly -------------------------------------- //------------------------------------------------------------------------------ @@ -39,10 +45,10 @@ Rpoly::Calculate( const double* op, int degree, double* zeror, double* zeroi) // The following statements set machine constants. - base = _DBL_RADIX ; - eta = DBL_EPSILON ; - infin = DBL_MAX ; - smalno = DBL_MIN ; + base = RADIX_DOUBLE ; + eta = EPS_DOUBLE ; + infin = BIG_DOUBLE ; + smalno = SMALL_DOUBLE ; are = eta ; mre = eta ; lo = smalno / eta ; @@ -264,6 +270,7 @@ Rpoly::fxshfr( int l2, int* nz) bool bSpass ; bool bVtry ; bool bStry ; + bool bFflag ; bool bIflag ; int type, i, j ; double svu, svv, ui, vi, s ; @@ -273,6 +280,7 @@ Rpoly::fxshfr( int l2, int* nz) // Inizializzazioni + bIflag = true ; *nz = 0 ; betav = 0.25 ; betas = 0.25 ; @@ -285,6 +293,8 @@ Rpoly::fxshfr( int l2, int* nz) for ( j = 0 ; j < l2 ; j ++) { + bFflag = true ; + // Calculate next k polynomial and estimate v. nextk( &type) ; calcsc( &type) ; @@ -339,14 +349,14 @@ Rpoly::fxshfr( int l2, int* nz) bVtry = false ; bStry = false ; while ( true) { - bIflag = true ; - if ( bSpass && ! bVpass || tss < tvv) + if ( bFflag && bSpass && ( ! bVpass || tss < tvv)) ; else { quadit( &ui, &vi, nz) ; if ( *nz > 0) return ; // Quadratic iteration has failed. Flag that it has been tried and decrease the convergence criterion. + bIflag = true ; bVtry = true ; betav *= 0.25; // Try linear iteration if it has not been tried and the S sequence is converging. @@ -357,6 +367,7 @@ Rpoly::fxshfr( int l2, int* nz) k[i] = svk[i] ; } } + bFflag = false ; if ( bIflag) { realit( s, nz, &bIflag) ; if ( *nz > 0) @@ -1383,8 +1394,8 @@ Cpoly::cmod( const double r, const double i) void Cpoly::mcon( double* eta, double* infiny, double* smalno, double* base) { - *base = _DBL_RADIX ; - *eta = DBL_EPSILON ; - *infiny = DBL_MAX ; - *smalno = DBL_MIN ; + *base = RADIX_DOUBLE ; + *eta = EPS_DOUBLE ; + *infiny = BIG_DOUBLE ; + *smalno = SMALL_DOUBLE ; } diff --git a/JenkinsTraub.h b/JenkinsTraub.h index a7b567f..49b84a7 100644 --- a/JenkinsTraub.h +++ b/JenkinsTraub.h @@ -31,7 +31,7 @@ class Rpoly { void realit( double sss, int* nz, bool* pIflag) ; void calcsc( int* type) ; void nextk( int* type) ; - void newest( int type, double* uu,double* vv) ; + void newest( int type, double* uu, double* vv) ; void quadsd( int n, double* u, double* v, double* p, double* q, double* a, double* b) ; diff --git a/Polynomial.cpp b/Polynomial.cpp new file mode 100644 index 0000000..1149f59 --- /dev/null +++ b/Polynomial.cpp @@ -0,0 +1,219 @@ +//---------------------------------------------------------------------------- +// EgalTech 2013-2014 +//---------------------------------------------------------------------------- +// File : PolynomialPoint3d.cpp Data : 12.01.14 Versione : 1.5a2 +// Contenuto : Implementazione classe polinomio con coefficienti Point3d. +// +// +// +// Modifiche : 12.01.14 DS Creazione modulo. +// +// +//---------------------------------------------------------------------------- + +//--------------------------- Include ---------------------------------------- +#include "stdafx.h" +#include "/EgtDev/Include/ENkPolynomial.h" +#include "/EgtDev/Include/ENkPolynomialRoots.h" + + +//---------------------------------------------------------------------------- +bool +Polynomial::SetDegree( int nDegree) +{ + // gradi negativi non hanno senso + if ( nDegree < 0) + return false ; + // pulisco, alloco e inizializzo a 0 + try { + m_Coeff.clear() ; + m_Coeff.reserve( nDegree + 1) ; + m_nDegree = nDegree ; + for ( int i = 0 ; i <= m_nDegree ; ++ i) + m_Coeff.push_back( 0) ; + return true ; + } + catch (...) { + return false ; + } +} + +//---------------------------------------------------------------------------- +bool +Polynomial::EnsureDegree( int nDegree) +{ + // se il grado è già adeguato non devo fare alcunché + if ( nDegree <= m_nDegree) + return true ; + // alloco e inizializzo a 0 i nuovi coefficienti + try { + m_Coeff.reserve( nDegree + 1) ; + for ( int i = m_nDegree + 1 ; i <= nDegree ; ++ i) + m_Coeff.push_back( 0) ; + m_nDegree = nDegree ; + return true ; + } + catch (...) { + return false ; + } +} + +//---------------------------------------------------------------------------- +bool +Polynomial::SetCoeff( int nPower, double dC) +{ + if ( nPower < 0 || nPower > m_nDegree) + return false ; + + // In posizione 0 il termine noto, in 1 il coefficiente della prima potenza,... + m_Coeff[nPower] = dC ; + return true ; +} + +//---------------------------------------------------------------------------- +bool +Polynomial::Set( int nDegree, const DBLVECTOR& vC) +{ + if ( ! SetDegree( nDegree)) + return false ; + + for ( int i = 0 ; i <= nDegree ; i ++) + m_Coeff[i] = vC[i] ; + return true ; +} + +//---------------------------------------------------------------------------- +bool +Polynomial::SetToConstant( double dC) +{ + if ( ! SetDegree( 0)) + return false ; + m_Coeff[0] = dC ; + return true ; +} + +//---------------------------------------------------------------------------- +const Polynomial& +Polynomial::operator +=( const Polynomial& polP) +{ + // mi assicuro che il polinomio risultante abbia grado sufficiente + EnsureDegree( polP.GetDegree()) ; + + // eseguo la somma + for ( int i = 0 ; i <= polP.GetDegree() ; ++ i) + m_Coeff[i] += polP.m_Coeff[i] ; + + return *this ; +} + +//---------------------------------------------------------------------------- +const Polynomial& +Polynomial::operator -=( const Polynomial& polP) +{ + // mi assicuro che il polinomio risultante abbia grado sufficiente + EnsureDegree( polP.GetDegree()) ; + + // eseguo la somma + for ( int i = 0 ; i <= polP.GetDegree() ; ++ i) + m_Coeff[i] += ( - polP.m_Coeff[i]) ; + + return *this ; +} + +//---------------------------------------------------------------------------- +const Polynomial& +Polynomial::operator *=( const Polynomial& polP) +{ + // copio il polinomio corrente + Polynomial polC = *this ; + + // pulisco e imposto il grado del polinomio risultante + SetDegree( polC.GetDegree() + polP.GetDegree()) ; + + // eseguo il prodotto + for ( int i = 0 ; i <= polC.GetDegree() ; ++ i) { + for ( int j = 0 ; j <= polP.GetDegree() ; ++ j) + m_Coeff[i+j] += polC.m_Coeff[i] * polP.m_Coeff[j] ; + } + + return *this ; +} + +//---------------------------------------------------------------------------- +void +Polynomial::Derive( void) +{ + // polinomio non inizializzato + if ( m_nDegree < 0) + return ; + // polinomio costante + if ( m_nDegree == 0) { + m_Coeff[0] = 0 ; + return ; + } + // caso normale + for ( int i = 0 ; i < m_nDegree ; ++ i) + m_Coeff[i] = ( i + 1) * m_Coeff[i+1] ; + m_Coeff.pop_back() ; + -- m_nDegree ; +} + +//---------------------------------------------------------------------------- +void +Polynomial::Derive( const Polynomial& polP) +{ + operator=( polP) ; + Derive() ; +} + +//---------------------------------------------------------------------------- +void +Polynomial::AdjustDegree( void) +{ + // se il coefficiente del grado più alto è zero, diminuisco il grado + while ( m_nDegree >= 0 && fabs( m_Coeff[m_nDegree]) < DBL_EPSILON) { + m_Coeff.pop_back() ; + -- m_nDegree ; + } +} + +//---------------------------------------------------------------------------- +double +Polynomial::Evaluate( double dVal) +{ + // polinomio non inizializzato + if ( m_nDegree < 0) + return 0 ; + // caso normale + double dRes = m_Coeff[m_nDegree] ; + for ( int i = m_nDegree - 1 ; i >= 0 ; -- i) + dRes = dRes * dVal + m_Coeff[i] ; + + return dRes ; +} + +//---------------------------------------------------------------------------- +int +Polynomial::FindRoots( DBLVECTOR& vdRoot) +{ + return PolynomialRoots( m_nDegree, m_Coeff, vdRoot) ; +} + +//---------------------------------------------------------------------------- +int +FilterMultipleAndOutOfRangeRoots( DBLVECTOR& vRoots, double dMin, double dMax, double dEps) +{ + int nZ = vRoots.size() ; + for ( int i = 0 ; i < nZ ;) { + if ( vRoots[i] < dMin || vRoots[i] > dMax || + ( i >= 1 && fabs( vRoots[i]- vRoots[i-1]) < dEps)) { + nZ -- ; + for ( int j = i ; j < nZ ; ++ j) + vRoots[j] = vRoots[j+1] ; + vRoots.pop_back() ; + } + else + ++ i ; + } + return nZ ; +} diff --git a/PolynomialZeros.cpp b/PolynomialRoots.cpp similarity index 76% rename from PolynomialZeros.cpp rename to PolynomialRoots.cpp index 1dc4541..02c5dfe 100644 --- a/PolynomialZeros.cpp +++ b/PolynomialRoots.cpp @@ -1,8 +1,8 @@ //---------------------------------------------------------------------------- // EgalTech 2013-2013 //---------------------------------------------------------------------------- -// File : PolynomialZeros.cpp Data : 08.01.14 Versione : 1.5a1 -// Contenuto : Funzione per il calcolo degli zeri di polinomi. +// File : PolynomialRoots.cpp Data : 08.01.14 Versione : 1.5a2 +// Contenuto : Funzioni per il calcolo degli zeri di polinomi. // // // @@ -14,7 +14,7 @@ //--------------------------- Include ---------------------------------------- #include "stdafx.h" #include "JenkinsTraub.h" -#include "\EgtDev\Include\ENkPolynomialZeros.h" +#include "/EgtDev/Include/ENkPolynomialRoots.h" #include @@ -25,10 +25,9 @@ static void SortRoots( int nNum, Complex acRoot[]) ; //---------------------------------------------------------------------------- int -PolynomialZeros( int nDegree, double adPoly[], double adRoot[], int* pnIter) +PolynomialRoots( int nDegree, DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter) { int i ; - int j ; int nZeros ; double dPreal[POLY_MAXDEG+1] ; double dZreal[POLY_MAXDEG] ; @@ -40,8 +39,12 @@ PolynomialZeros( int nDegree, double adPoly[], double adRoot[], int* pnIter) if ( pnIter != NULL) *pnIter = 0 ; + // controllo che il vettore dei coefficienti sia lungo almeno come il grado + 1 + if ( (int) vdPoly.size() <= nDegree) + return 0 ; + // se il coefficiente del grado più alto è zero, diminuisco il grado - while ( nDegree >= 0 && fabs( adPoly[nDegree]) < DBL_EPSILON) + while ( nDegree >= 0 && fabs( vdPoly[nDegree]) < DBL_EPSILON) nDegree -- ; // se il grado è nullo o negativo, errore @@ -54,22 +57,23 @@ PolynomialZeros( int nDegree, double adPoly[], double adRoot[], int* pnIter) // riordino i coefficienti reali for ( i = 0 ; i <= nDegree ; i++) - dPreal[i] = adPoly[nDegree-i] ; + dPreal[i] = vdPoly[nDegree-i] ; // calcolo gli zeri nZeros = cRpoly.Calculate( dPreal, nDegree, dZreal, dZcplx) ; // assegno gli zeri reali ai parametri di ritorno - for ( i = 0, j = 0 ; i < nZeros ; i++) { + vdRoot.clear() ; + vdRoot.reserve( nZeros) ; + for ( i = 0 ; i < nZeros ; i++) { if ( fabs( dZcplx[i]) < 100 * DBL_EPSILON) { - adRoot[j] = dZreal[i] ; - j ++ ; + vdRoot.push_back( dZreal[i]) ; } } - nZeros = j ; + nZeros = vdRoot.size() ; // ordino le radici in senso decrescente - SortRoots( nZeros, adRoot) ; + SortRoots( nZeros, vdRoot.data()) ; // assegno il numero di iterazioni if ( pnIter != NULL) @@ -80,7 +84,7 @@ PolynomialZeros( int nDegree, double adPoly[], double adRoot[], int* pnIter) //---------------------------------------------------------------------------- int -PolynomialZeros( int nDegree, Complex acPoly[], Complex acRoot[], int* pnIter) +PolynomialRoots( int nDegree, CPLXVECTOR& vcPoly, CPLXVECTOR& vcRoot, int* pnIter) { bool bCplx ; int i ; @@ -97,8 +101,12 @@ PolynomialZeros( int nDegree, Complex acPoly[], Complex acRoot[], int* pnIter) if ( pnIter != NULL) *pnIter = 0 ; + // controllo che il vettore dei coefficienti sia lungo almeno come il grado + 1 + if ( (int) vcPoly.size() <= nDegree) + return 0 ; + // se il coefficiente del grado più alto è zero, diminuisco il grado - while ( nDegree >= 0 && m2( acPoly[nDegree]) < DBL_EPSILON * DBL_EPSILON) + while ( nDegree >= 0 && m2( vcPoly[nDegree]) < DBL_EPSILON * DBL_EPSILON) nDegree -- ; // se il grado è nullo o negativo, errore @@ -111,12 +119,12 @@ PolynomialZeros( int nDegree, Complex acPoly[], Complex acRoot[], int* pnIter) // ricavo i coefficienti reali for ( i = 0 ; i <= nDegree ; i++) - dPreal[i] = acPoly[nDegree-i].re ; + dPreal[i] = vcPoly[nDegree-i].re ; // ricavo i coefficienti complessi ( e verifico se non nulli) bCplx = false ; for ( i = 0 ; i <= nDegree ; i++) { - dPcplx[i] = acPoly[nDegree-i].im ; + dPcplx[i] = vcPoly[nDegree-i].im ; if ( fabs( dPcplx[i]) > DBL_EPSILON) bCplx = true ; } @@ -128,21 +136,22 @@ PolynomialZeros( int nDegree, Complex acPoly[], Complex acRoot[], int* pnIter) nZeros = cRpoly.Calculate( dPreal, nDegree, dZreal, dZcplx) ; // assegno gli zeri ai parametri di ritorno + vcRoot.clear() ; + vcRoot.reserve( nZeros) ; for ( i = 0 ; i < nZeros ; i++) { - acRoot[i].re = dZreal[i] ; - acRoot[i].im = dZcplx[i] ; + vcRoot.push_back( Complex( dZreal[i], dZcplx[i])) ; } // annullo le parti reali e immaginarie molto piccole for ( i = 0 ; i < nZeros ; i++) { - if ( fabs( acRoot[i].re) < 100 * DBL_EPSILON) - acRoot[i].re = 0 ; - if ( fabs( acRoot[i].im) < 100 * DBL_EPSILON) - acRoot[i].im = 0 ; + if ( fabs( vcRoot[i].re) < 100 * DBL_EPSILON) + vcRoot[i].re = 0 ; + if ( fabs( vcRoot[i].im) < 100 * DBL_EPSILON) + vcRoot[i].im = 0 ; } // ordino le radici in senso decrescente della parte reale - SortRoots( nZeros, acRoot) ; + SortRoots( nZeros, vcRoot.data()) ; // assegno il numero di iterazioni if ( pnIter != NULL)