EgtNumKernel/Polynomial.cpp

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

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

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

//----------------------------------------------------------------------------
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  &&  abs( 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 = (int) vRoots.size() ;
   for ( int i = 0 ; i < nZ ;) {
      if ( vRoots[i] < dMin  ||  vRoots[i] > dMax  ||
           ( i >= 1  &&  abs( 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 ;
}