EgtNumKernel/PolynomialRoots.cpp

//----------------------------------------------------------------------------
//  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 ;

//----------------------------------------------------------------------------
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 ;

  // 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 ;
}