Dario Sassi
f19e7827fd
- corretta soluzione di equazione quadratica con coefficiente di x^2 molto piccolo.
230 lines
7.1 KiB
C++
230 lines
7.1 KiB
C++
//----------------------------------------------------------------------------
|
|
// EgalTech 2014-2018
|
|
//----------------------------------------------------------------------------
|
|
// File : PolynomialRoots.cpp Data : 07.08.18 Versione : 1.9h1
|
|
// Contenuto : Funzioni per il calcolo degli zeri di polinomi.
|
|
//
|
|
//
|
|
//
|
|
// Modifiche : 08.01.14 DS Creazione modulo.
|
|
//
|
|
//
|
|
//----------------------------------------------------------------------------
|
|
|
|
//--------------------------- Include ----------------------------------------
|
|
#include "stdafx.h"
|
|
#include "JenkinsTraub.h"
|
|
#include "/EgtDev/Include/ENkPolynomialRoots.h"
|
|
#include <algorithm>
|
|
|
|
using namespace std ;
|
|
|
|
//----------------------------------------------------------------------------
|
|
static int LinearPolynomialRoots( DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter) ;
|
|
static int QuadraticPolynomialRoots( DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter) ;
|
|
|
|
//----------------------------------------------------------------------------
|
|
int
|
|
PolynomialRoots( int nDegree, DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter)
|
|
{
|
|
// inizializzo il numero di iterazioni
|
|
if ( pnIter != nullptr)
|
|
*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 && abs( vdPoly[nDegree]) < DBL_EPSILON)
|
|
nDegree -- ;
|
|
|
|
// se il grado è nullo o negativo, errore
|
|
if ( nDegree <= 0)
|
|
return 0 ;
|
|
|
|
// verifico di non superare il massimo grado ammesso
|
|
if ( nDegree > POLY_MAXDEG)
|
|
return 0 ;
|
|
|
|
// se polinomio lineare
|
|
if ( nDegree == 1)
|
|
return LinearPolynomialRoots( vdPoly, vdRoot, pnIter) ;
|
|
|
|
// se polinomi quadratico
|
|
if ( nDegree == 2)
|
|
return QuadraticPolynomialRoots( vdPoly, vdRoot, pnIter) ;
|
|
|
|
// riordino i coefficienti reali
|
|
double dPreal[POLY_MAXDEG+1] ;
|
|
for ( int i = 0 ; i <= nDegree ; i++)
|
|
dPreal[i] = vdPoly[nDegree-i] ;
|
|
|
|
// calcolo gli zeri
|
|
Rpoly cRpoly ;
|
|
double dZreal[POLY_MAXDEG] ;
|
|
double dZcplx[POLY_MAXDEG] ;
|
|
int nZeros = cRpoly.Calculate( dPreal, nDegree, dZreal, dZcplx) ;
|
|
|
|
// assegno gli zeri reali ai parametri di ritorno
|
|
vdRoot.clear() ;
|
|
vdRoot.reserve( nZeros) ;
|
|
for ( int i = 0 ; i < nZeros ; i++) {
|
|
if ( abs( dZcplx[i]) < 100 * DBL_EPSILON) {
|
|
vdRoot.push_back( dZreal[i]) ;
|
|
}
|
|
}
|
|
nZeros = (int) vdRoot.size() ;
|
|
|
|
// ordino le radici in senso decrescente
|
|
sort( vdRoot.begin(), vdRoot.end(),
|
|
[]( const double& a, const double& b) { return ( a > b) ; }) ;
|
|
|
|
// assegno il numero di iterazioni
|
|
if ( pnIter != nullptr)
|
|
*pnIter = cRpoly.itercnt ;
|
|
|
|
return nZeros ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int
|
|
PolynomialRoots( int nDegree, CPLXVECTOR& vcPoly, CPLXVECTOR& vcRoot, int* pnIter)
|
|
{
|
|
// inizializzo il numero di iterazioni
|
|
if ( pnIter != nullptr)
|
|
*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 && norm( vcPoly[nDegree]) < DBL_EPSILON * DBL_EPSILON)
|
|
nDegree -- ;
|
|
|
|
// se il grado è nullo o negativo, errore
|
|
if ( nDegree <= 0)
|
|
return 0 ;
|
|
|
|
// verifico di non superare il massimo grado ammesso
|
|
if ( nDegree > POLY_MAXDEG)
|
|
return 0 ;
|
|
|
|
// ricavo i coefficienti reali
|
|
double dPreal[POLY_MAXDEG+1] ;
|
|
for ( int i = 0 ; i <= nDegree ; i++)
|
|
dPreal[i] = real( vcPoly[nDegree-i]) ;
|
|
|
|
// ricavo i coefficienti immaginari ( e verifico se non nulli)
|
|
bool bImag = false ;
|
|
double dPimag[POLY_MAXDEG+1] ;
|
|
for ( int i = 0 ; i <= nDegree ; i++) {
|
|
dPimag[i] = imag( vcPoly[nDegree-i]) ;
|
|
if ( abs( dPimag[i]) > DBL_EPSILON)
|
|
bImag = true ;
|
|
}
|
|
|
|
// calcolo gli zeri
|
|
int nZeros ;
|
|
int nIterCnt ;
|
|
double dZreal[POLY_MAXDEG] ;
|
|
double dZcplx[POLY_MAXDEG] ;
|
|
if ( bImag) {
|
|
Cpoly cCpoly ;
|
|
nZeros = cCpoly.Calculate( dPreal, dPimag, nDegree, dZreal, dZcplx) ;
|
|
nIterCnt = cCpoly.itercnt ;
|
|
}
|
|
else {
|
|
Rpoly cRpoly ;
|
|
nZeros = cRpoly.Calculate( dPreal, nDegree, dZreal, dZcplx) ;
|
|
nIterCnt = cRpoly.itercnt ;
|
|
}
|
|
|
|
// assegno gli zeri ai parametri di ritorno
|
|
vcRoot.clear() ;
|
|
vcRoot.reserve( nZeros) ;
|
|
for ( int i = 0 ; i < nZeros ; i++) {
|
|
vcRoot.push_back( complex<double>( dZreal[i], dZcplx[i])) ;
|
|
}
|
|
|
|
// annullo le parti reali e immaginarie molto piccole
|
|
for ( int i = 0 ; i < nZeros ; i++) {
|
|
if ( abs( real( vcRoot[i])) < 100 * DBL_EPSILON)
|
|
vcRoot[i].real( 0) ;
|
|
if ( abs( imag( vcRoot[i])) < 100 * DBL_EPSILON)
|
|
vcRoot[i].imag( 0) ;
|
|
}
|
|
|
|
// ordino le radici in senso decrescente della parte reale e in subordine della immaginaria
|
|
sort( vcRoot.begin(), vcRoot.end(),
|
|
[]( const complex<double>& a, const complex<double>& b) {
|
|
// se parti reali identiche, confronto quelle immaginarie
|
|
if ( abs( real( a) - real( b)) < FLT_MIN)
|
|
return ( imag( a) > imag( b)) ;
|
|
else
|
|
return ( real( a) > real( b)) ; }) ;
|
|
|
|
// assegno il numero di iterazioni
|
|
if ( pnIter != nullptr)
|
|
*pnIter = nIterCnt ;
|
|
|
|
return nZeros ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
static int
|
|
LinearPolynomialRoots( DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter)
|
|
{
|
|
// inizializzazioni
|
|
vdRoot.clear() ;
|
|
vdRoot.reserve( 1) ;
|
|
if ( pnIter != nullptr)
|
|
*pnIter = 1 ;
|
|
// il coefficiente del grado più alto vdPoly[1] è sicuramente non nullo
|
|
vdRoot.push_back( - vdPoly[0] / vdPoly[1]) ;
|
|
return 1 ;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
static int
|
|
QuadraticPolynomialRoots( DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter)
|
|
{
|
|
// inizializzazioni
|
|
vdRoot.clear() ;
|
|
vdRoot.reserve( 2) ;
|
|
if ( pnIter != nullptr)
|
|
*pnIter = 1 ;
|
|
// il coefficiente del grado più alto vdPoly[2] è sicuramente non nullo
|
|
// caso omogeneo( termine noto nullo)
|
|
if ( abs( vdPoly[0]) < DBL_EPSILON) {
|
|
vdRoot.push_back( 0) ;
|
|
vdRoot.push_back( - vdPoly[1] / vdPoly[2]) ;
|
|
if ( vdRoot[0] < vdRoot[1])
|
|
swap( vdRoot[0], vdRoot[1]) ;
|
|
return 2 ;
|
|
}
|
|
// caso generale ( x^2 + 2 p x + q = 0)
|
|
double p = vdPoly[1] / ( 2.0 * vdPoly[2]) ;
|
|
double q = vdPoly[0] / vdPoly[2] ;
|
|
double Delta = p * p - q ;
|
|
if ( abs( Delta) < DBL_EPSILON * DBL_EPSILON) {
|
|
vdRoot.clear() ;
|
|
vdRoot.reserve( 1) ;
|
|
vdRoot.push_back( -p) ;
|
|
return 1 ;
|
|
}
|
|
if ( Delta < 0.0) {
|
|
return 0 ;
|
|
}
|
|
// Delta positivo
|
|
double dSqrtD = sqrt( Delta) ;
|
|
double r1 = -( p + copysign( dSqrtD, p)) ;
|
|
double r2 = q / r1 ;
|
|
if ( r1 > r2)
|
|
swap( r1, r2) ;
|
|
vdRoot.push_back( r1) ;
|
|
vdRoot.push_back( r2) ;
|
|
return 2 ;
|
|
}
|