Egalware/EgtNumKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

EgtNumKernel 1.5a2 : Aggiunta classe Polynomial.

Browse Source 
Corretto calcolo zeri di polinomi reali.
This commit is contained in:
Dario Sassi
2014-01-12 22:50:16 +00:00
parent 44b5a1007b
commit b6af673186
7 changed files with 295 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files
EgtNumKernel.rc
BIN
View File
Binary file not shown.
EgtNumKernel.vcxproj
+7 -4
View File
@@ -87,8 +87,7 @@ copy $(TargetPath) \EgtProg\Dll</Command>
<OptimizeReferences>true</OptimizeReferences>
</Link>
<PostBuildEvent>
<Command>copy $(TargetDir)$(TargetName).pdb \EgtDev\Lib\
copy $(TargetDir)$(TargetName).lib \EgtDev\Lib\
<Command>copy $(TargetDir)$(TargetName).lib \EgtDev\Lib\
copy $(TargetPath) \EgtProg\Dll</Command>
</PostBuildEvent>
<ResourceCompile>
@@ -99,7 +98,8 @@ copy $(TargetPath) \EgtProg\Dll</Command>
<ClCompile Include="Complex.cpp" />
<ClCompile Include="ENkDllMain.cpp" />
<ClCompile Include="JenkinsTraub.cpp" />
<ClCompile Include="PolynomialZeros.cpp" />
<ClCompile Include="Polynomial.cpp" />
<ClCompile Include="PolynomialRoots.cpp" />
<ClCompile Include="stdafx.cpp">
<PrecompiledHeader Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">Create</PrecompiledHeader>
<PrecompiledHeader Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">Create</PrecompiledHeader>
@@ -109,8 +109,11 @@ copy $(TargetPath) \EgtProg\Dll</Command>
<ClInclude Include="..\Include\EgnGetModuleVer.h" />
<ClInclude Include="..\Include\EgtTargetVer.h" />
<ClInclude Include="..\Include\EgtTrace.h" />
<ClInclude Include="..\Include\ENkComplex.h" />
<ClInclude Include="..\Include\ENkDllMain.h" />
<ClInclude Include="..\Include\ENkPolynomialZeros.h" />
<ClInclude Include="..\Include\ENkNumCollection.h" />
<ClInclude Include="..\Include\ENkPolynomial.h" />
<ClInclude Include="..\Include\ENkPolynomialRoots.h" />
<ClInclude Include="DllMain.h" />
<ClInclude Include="JenkinsTraub.h" />
<ClInclude Include="resource.h" />
EgtNumKernel.vcxproj.filters
+15 -3
View File
@@ -27,7 +27,10 @@
<ClCompile Include="Complex.cpp">
<Filter>File di origine</Filter>
</ClCompile>
<ClCompile Include="PolynomialZeros.cpp">
<ClCompile Include="Polynomial.cpp">
<Filter>File di origine</Filter>
</ClCompile>
<ClCompile Include="PolynomialRoots.cpp">
<Filter>File di origine</Filter>
</ClCompile>
</ItemGroup>
@@ -53,10 +56,19 @@
<ClInclude Include="..\Include\EgtTargetVer.h">
<Filter>File di intestazione</Filter>
</ClInclude>
<ClInclude Include="..\Include\ENkPolynomialZeros.h">
<ClInclude Include="JenkinsTraub.h">
<Filter>File di intestazione</Filter>
</ClInclude>
<ClInclude Include="JenkinsTraub.h">
<ClInclude Include="..\Include\ENkComplex.h">
<Filter>File di intestazione</Filter>
</ClInclude>
<ClInclude Include="..\Include\ENkNumCollection.h">
<Filter>File di intestazione</Filter>
</ClInclude>
<ClInclude Include="..\Include\ENkPolynomialRoots.h">
<Filter>File di intestazione</Filter>
</ClInclude>
<ClInclude Include="..\Include\ENkPolynomial.h">
<Filter>File di intestazione</Filter>
</ClInclude>
</ItemGroup>
JenkinsTraub.cpp
+21 -10
View File
@@ -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 ;
}
JenkinsTraub.h
+1 -1
View File
@@ -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) ;
Polynomial.cpp
+219
View File
@@ -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 ;
}
PolynomialZeros.cpp → PolynomialRoots.cpp
+32 -23
View File
@@ -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 <stdlib.h>
@@ -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)