EgtNumKernel/PolynomialRoots.cpp

//----------------------------------------------------------------------------
//  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 <stdlib.h>


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