diff --git a/EgtNumKernel.rc b/EgtNumKernel.rc
index a8d3075..739fc70 100644
Binary files a/EgtNumKernel.rc and b/EgtNumKernel.rc differ
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)