//---------------------------------------------------------------------------- // EgalTech 2013-2013 //---------------------------------------------------------------------------- // File : PolynomialRoots.cpp Data : 08.01.14 Versione : 1.5a2 // 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 //--------------------------------- Prototipi locali -------------------------------- static void SortRoots( int nNum, double adRoot[]) ; static void SortRoots( int nNum, Complex acRoot[]) ; //---------------------------------------------------------------------------- int PolynomialRoots( int nDegree, DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter) { int i ; int nZeros ; double dPreal[POLY_MAXDEG+1] ; double dZreal[POLY_MAXDEG] ; double dZcplx[POLY_MAXDEG] ; Rpoly cRpoly ; // inizializzo il numero di iterazioni 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( 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 ; // riordino i coefficienti reali for ( i = 0 ; i <= 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 vdRoot.clear() ; vdRoot.reserve( nZeros) ; for ( i = 0 ; i < nZeros ; i++) { if ( fabs( dZcplx[i]) < 100 * DBL_EPSILON) { vdRoot.push_back( dZreal[i]) ; } } nZeros = (int) vdRoot.size() ; // ordino le radici in senso decrescente SortRoots( nZeros, vdRoot.data()) ; // assegno il numero di iterazioni if ( pnIter != NULL) *pnIter = cRpoly.itercnt ; return nZeros ; } //---------------------------------------------------------------------------- int PolynomialRoots( int nDegree, CPLXVECTOR& vcPoly, CPLXVECTOR& vcRoot, int* pnIter) { bool bCplx ; int i ; int nZeros ; double dPreal[POLY_MAXDEG+1] ; double dPcplx[POLY_MAXDEG+1] ; double dZreal[POLY_MAXDEG] ; double dZcplx[POLY_MAXDEG] ; Rpoly cRpoly ; Cpoly cCpoly ; // inizializzo il numero di iterazioni 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( 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 for ( i = 0 ; i <= nDegree ; i++) 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] = vcPoly[nDegree-i].im ; if ( fabs( dPcplx[i]) > DBL_EPSILON) bCplx = true ; } // calcolo gli zeri if ( bCplx) nZeros = cCpoly.Calculate( dPreal, dPcplx, nDegree, dZreal, dZcplx) ; else 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++) { vcRoot.push_back( Complex( dZreal[i], dZcplx[i])) ; } // annullo le parti reali e immaginarie molto piccole for ( i = 0 ; i < nZeros ; i++) { 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, vcRoot.data()) ; // assegno il numero di iterazioni if ( pnIter != NULL) *pnIter = ( bCplx ? cCpoly.itercnt : cRpoly.itercnt) ; return nZeros ; } //----------------------------------------------------------------------------- // Confronto tra numeri reali per ordinarli secondo l'ordine crescente //----------------------------------------------------------------------------- int CompareRealRoots( const void* pRoot1, const void* pRoot2) { double dRe1 ; double dRe2 ; // valori reali dRe1 = *(double*) pRoot1 ; dRe2 = *(double*) pRoot2 ; // se primo maggiore del secondo if ( dRe1 > dRe2) return - 1 ; // se primo minore del secondo else if ( dRe1 < dRe2) return + 1 ; // altrimenti uguali else return 0 ; } //----------------------------------------------------------------------------- void SortRoots( int nNum, double adRoot[]) { if ( nNum <= 0) return ; qsort( adRoot, size_t( nNum), sizeof( double), CompareRealRoots) ; } //----------------------------------------------------------------------------- // Confronto tra numeri complessi per ordinarli secondo l'ordine crescente // delle parti reali //----------------------------------------------------------------------------- int CompareComplexRoots( const void* pRoot1, const void* pRoot2) { double dRe1 ; double dRe2 ; double dIm1 ; double dIm2 ; // parti reali dRe1 = Re( *(Complex*) pRoot1) ; dRe2 = Re( *(Complex*) pRoot2) ; // se parti reali praticamente uguali if ( fabs( dRe1 - dRe2) < FLT_MIN) { // parti immaginarie dIm1 = Im( *(Complex*) pRoot1) ; dIm2 = Im( *(Complex*) pRoot2) ; if ( dIm1 > dIm2) return - 1 ; else if ( dIm1 < dIm2) return + 1 ; else return 0 ; } // se primo maggiore del secondo else if ( dRe1 > dRe2) return - 1 ; // altrimenti secondo maggiore del primo else return + 1 ; } //----------------------------------------------------------------------------- void SortRoots( int nNum, Complex acRoot[]) { if ( nNum <= 0) return ; qsort( acRoot, size_t( nNum), sizeof( Complex), CompareComplexRoots) ; }