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