diff --git a/Complex.cpp b/Complex.cpp deleted file mode 100644 index 2652320..0000000 --- a/Complex.cpp +++ /dev/null @@ -1,271 +0,0 @@ -//---------------------------------------------------------------------------- -// EgalTech 2013-2013 -//---------------------------------------------------------------------------- -// File : Complex.cpp Data : 08.01.14 Versione : 1.5a1 -// Contenuto : Implementazione classe dei numeri complessi. -// -// -// -// Modifiche : 08.01.14 DS Creazione modulo. -// -// -//---------------------------------------------------------------------------- - -//--------------------------- Include ---------------------------------------- -#include "stdafx.h" -#include "\EgtDev\Include\ENkComplex.h" - - -//---------------------------- Classe Complex --------------------------------- -Complex -Complex::operator +=( double dVal) -{ - this->re += dVal ; - return *this; -} - -//---------------------------------------------------------------------------- -Complex -Complex::operator +=( Complex& cVal) -{ - this->re += cVal.re; - this->im += cVal.im; - - return *this; -} - -//---------------------------------------------------------------------------- -Complex -Complex::operator -=( double dVal) -{ - this->re -= dVal ; - return *this ; -} - -//---------------------------------------------------------------------------- -Complex -Complex::operator -=( Complex& cVal) -{ - this->re -= cVal.re ; - this->im -= cVal.im ; - - return *this ; -} - -//---------------------------------------------------------------------------- -Complex -Complex::operator *=( double dVal) -{ - this->re *= dVal ; - this->im *= dVal ; - - return *this; -} - -//---------------------------------------------------------------------------- -Complex -Complex::operator /=( double dVal) -{ - double dInv ; - - - dInv = 1.0 / dVal ; - this->re *= dInv ; - this->im *= dInv ; - - return *this ; -} - -//---------------------------------------------------------------------------- -Complex -Complex::operator >>=( int n) -{ - this->re = ldexp( this->re, -n) ; - this->im = ldexp( this->im, -n) ; - - return *this ; -} - -//---------------------------------------------------------------------------- -Complex -Complex::operator <<=( int n) -{ - this->re = ldexp( this->re, +n) ; - this->im = ldexp( this->im, +n) ; - - return *this ; -} - -//---------------------------------------------------------------------------- -Complex -Complex::operator *=( Complex& cVal) -{ - return ( *this = *this * cVal) ; -} - -//---------------------------------------------------------------------------- -Complex -Complex::operator /=( Complex& cVal) -{ - return ( *this *= inv( cVal)) ; -} - - -//------------------------------ Functions ----------------------------------- -// sqrt for Complex -Complex -sqrt( Complex& cVal) -{ - Complex z ; // Power 0.5 simple enough - double m ; // to do separate, faster - // than full exp(0.5*log(z)) - // just like reals have their - m = mod( cVal) ; // sqrt. Ours gives the one - z.re = sqrt( (m + cVal.re) / 2) ; // with -pi/2 < arg <= +pi/2 - z.im = sqrt( (m - cVal.re) / 2) ; // Our log interprets arg - if ( cVal.im < 0.) // as in range -pi to pi, - z.im = - z.im ; // like the atan2 used. - - return z ; -} - -//---------------------------------------------------------------------------- -// log for Complex -Complex -log( Complex& cVal) -{ - Complex z ; - - - z.re = log( m2( cVal)) / 2 ; - z.im = atan2( cVal.im, cVal.re) ; - - return z ; -} - -//---------------------------------------------------------------------------- -Complex -exp( Complex& cVal) -{ - Complex ez ; - double m ; - - - m = exp( cVal.re) ; - ez.re = m * cos( cVal.im) ; - ez.im = m * sin( cVal.im) ; - - return ez ; -} - -//---------------------------------------------------------------------------- -Complex -cosh( Complex& cVal) -{ - Complex ez ; - - - ez = exp( cVal) ; - - return ( ( ez + inv(ez)) >> 1) ; -} - -//---------------------------------------------------------------------------- -Complex -sinh( Complex& cVal) -{ - Complex ez ; - - - ez = exp( cVal) ; - - return ( ( ez - inv(ez)) >> 1) ; -} - -//---------------------------------------------------------------------------- -Complex -tanh( Complex& cVal) -{ - Complex e2z ; - - - e2z = exp( cVal << 1) ; - return ( ( e2z - 1) / ( e2z + 1)) ; -} - -//---------------------------------------------------------------------------- -Complex -cos( Complex& cVal) -{ - return cosh( itimes( cVal)) ; -} - -//---------------------------------------------------------------------------- -Complex -isin( Complex& cVal) -{ - return sinh( itimes( cVal)) ; -} - -//---------------------------------------------------------------------------- -Complex -sin( Complex& cVal) -{ - return -itimes( isin( cVal)) ; -}; - -//---------------------------------------------------------------------------- -Complex -itan( Complex& cVal) -{ - return tanh( itimes( cVal)) ; -} - -//---------------------------------------------------------------------------- -Complex -tan( Complex& cVal) -{ - return -itimes( itan( cVal)) ; -} - -//---------------------------------------------------------------------------- -Complex -acosh( Complex& cVal) -{ - return log( cVal + sqrt( cVal * cVal - 1)) ; -} - -//---------------------------------------------------------------------------- -Complex -asinh( Complex& cVal) -{ - return log( cVal + sqrt( cVal * cVal + 1)) ; -} - -//---------------------------------------------------------------------------- -Complex -atanh( Complex& cVal) -{ - return ( log(( 1 + cVal) / ( 1 - cVal)) >> 1) ; -} - -//---------------------------------------------------------------------------- -Complex -acos( Complex& cVal) -{ - return -itimes( acosh( cVal)) ; -} - -//---------------------------------------------------------------------------- -Complex -asin( Complex& cVal) -{ - return -itimes( asinh( itimes( cVal))) ; -} - -//---------------------------------------------------------------------------- -Complex -atan( Complex& cVal) -{ - return -itimes( atanh( itimes( cVal))) ; -} diff --git a/ENkDllMain.cpp b/ENkDllMain.cpp index 033f2d5..4102861 100644 --- a/ENkDllMain.cpp +++ b/ENkDllMain.cpp @@ -34,7 +34,7 @@ const int STR_DIM = 40 ; //----------------------------------------------------------------------------- -static HINSTANCE s_hModule = NULL ; +static HINSTANCE s_hModule = nullptr ; static char s_szENkNameVer[STR_DIM] ; //----------------------------------------------------------------------------- @@ -59,7 +59,7 @@ DllMain( HMODULE hModule, DWORD dwReason, LPVOID lpReserved) EGT_TRACE( "EgtNumKernel.dll Initializing!\n") ; } else if ( dwReason == DLL_PROCESS_DETACH) { - s_hModule = NULL ; + s_hModule = nullptr ; EGT_TRACE( "EgtNumKernel.dll Terminating!\n") ; } diff --git a/EgtNumKernel.rc b/EgtNumKernel.rc index 6b7d46d..d2bcd1b 100644 Binary files a/EgtNumKernel.rc and b/EgtNumKernel.rc differ diff --git a/EgtNumKernel.vcxproj b/EgtNumKernel.vcxproj index 3659d6b..e3e913e 100644 --- a/EgtNumKernel.vcxproj +++ b/EgtNumKernel.vcxproj @@ -193,7 +193,6 @@ copy $(TargetPath) \EgtProg\Dll64 - @@ -218,7 +217,6 @@ copy $(TargetPath) \EgtProg\Dll64 - diff --git a/EgtNumKernel.vcxproj.filters b/EgtNumKernel.vcxproj.filters index 0443847..6b2b646 100644 --- a/EgtNumKernel.vcxproj.filters +++ b/EgtNumKernel.vcxproj.filters @@ -27,9 +27,6 @@ File di origine - - File di origine\Polynomial - File di origine\Polynomial @@ -89,9 +86,6 @@ File di intestazione - - File di intestazione - File di intestazione diff --git a/JenkinsTraub.cpp b/JenkinsTraub.cpp index 2694ea8..1011c16 100644 --- a/JenkinsTraub.cpp +++ b/JenkinsTraub.cpp @@ -107,7 +107,7 @@ Rpoly::Calculate( const double* op, int degree, double* zeror, double* zeroi) max = 0.0 ; min = infin ; for ( i = 0 ; i <= n ; i++) { - x = fabs( p[i]) ; + x = abs( p[i]) ; if ( x > max) max = x ; if ( x != 0.0 && x < min) @@ -140,7 +140,7 @@ Rpoly::Calculate( const double* op, int degree, double* zeror, double* zeroi) // Compute lower bound on moduli of roots. for ( i = 0 ; i <= n ; i++) { - pt[i] = fabs( p[i]) ; + pt[i] = abs( p[i]) ; } pt[n] = - pt[n] ; // Compute upper estimate of bound. @@ -165,7 +165,7 @@ Rpoly::Calculate( const double* op, int degree, double* zeror, double* zeroi) // Do Newton iteration until x converges to two decimal places. dx = x ; - while ( fabs( dx / x) > 0.005) { + while ( abs( dx / x) > 0.005) { ff = pt[0] ; df = ff ; for ( i = 1 ; i < n ; i++) { @@ -198,7 +198,7 @@ Rpoly::Calculate( const double* op, int degree, double* zeror, double* zeroi) k[j] = t * k[j-1] + p[j] ; } k[0] = p[0] ; - bZerOk = ( fabs( k[n-1]) <= fabs( bb) * eta * 10.0) ; + bZerOk = ( abs( k[n-1]) <= abs( bb) * eta * 10.0) ; } else { // Use unscaled form of recurrence. @@ -317,9 +317,9 @@ Rpoly::fxshfr( int l2, int* nz) // Compute relative measures of convergence of s and v sequences. if ( vv != 0.0) - tv = fabs( ( vv - ovv) / vv) ; + tv = abs( ( vv - ovv) / vv) ; if ( ss != 0.0) - ts = fabs( ( ss - oss) / ss) ; + ts = abs( ( ss - oss) / ss) ; /* If decreasing, multiply two most recent convergence measures. */ tvv = 1.0 ; if ( tv < otv) @@ -436,23 +436,23 @@ Rpoly::quadit( double *uu, double *vv, int *nz) // Return if roots of the quadratic are real and not // close to multiple or nearly equal and of opposite sign. - if ( fabs( fabs( szr) - fabs( lzr)) > 0.01 * fabs( lzr)) + if ( abs( abs( szr) - abs( lzr)) > 0.01 * abs( lzr)) return ; // Evaluate polynomial by quadratic synthetic division. quadsd( n, &u, &v, p, qp, &a, &b) ; - mp = fabs( a - szr * b) + fabs( szi * b) ; + mp = abs( a - szr * b) + abs( szi * b) ; // Compute a rigorous bound on the rounding error in evaluating p. - zm = sqrt( fabs( v)) ; - ee = 2.0 * fabs( qp[0]) ; + zm = sqrt( abs( v)) ; + ee = 2.0 * abs( qp[0]) ; t = -szr * b ; for ( i = 1 ; i < n ; i ++) { - ee = ee * zm + fabs( qp[i]) ; + ee = ee * zm + abs( qp[i]) ; } - ee = ee * zm + fabs( a + t) ; + ee = ee * zm + abs( a + t) ; ee *= (5.0 * mre + 4.0 * are) ; - ee = ee - ( 5.0 * mre + 2.0 * are) * ( fabs( a + t) + fabs( b) * zm) ; - ee = ee + 2.0 * are * fabs( t) ; + ee = ee - ( 5.0 * mre + 2.0 * are) * ( abs( a + t) + abs( b) * zm) ; + ee = ee + 2.0 * are * abs( t) ; // Iteration has converged sufficiently if the polynomial value is less than 20 times this bound. if ( mp <= 20.0 * ee) { *nz = 2 ; @@ -493,7 +493,7 @@ Rpoly::quadit( double *uu, double *vv, int *nz) // If vi is zero the iteration is not converging. if ( vi == 0.0) return ; - relstp = fabs( ( vi - v) / vi) ; + relstp = abs( ( vi - v) / vi) ; u = ui ; v = vi ; } @@ -529,12 +529,12 @@ Rpoly::realit( double sss, int* nz, bool* pbIflag) pv = pv * s + p[i] ; qp[i] = pv ; } - mp = fabs( pv) ; + mp = abs( pv) ; // Compute a rigorous bound on the error in evaluating p. - ms = fabs( s) ; - ee = ( mre / ( are + mre)) * fabs( qp[0]) ; + ms = abs( s) ; + ee = ( mre / ( are + mre)) * abs( qp[0]) ; for ( i = 1 ; i <= n ; i ++) { - ee = ee * ms + fabs( qp[i]) ; + ee = ee * ms + abs( qp[i]) ; } // Iteration has converged sufficiently if the polynomial value is less than 20 times this bound. if ( mp <= 20.0 * (( are + mre) * ee - mre * mp)) { @@ -549,7 +549,7 @@ Rpoly::realit( double sss, int* nz, bool* pbIflag) return ; // A cluster of zeros near the real axis has been encountered. if ( j >= 2 && - ! ( fabs( t) > 0.001 * fabs( s-t) || mp < omp)) { + ! ( abs( t) > 0.001 * abs( s-t) || mp < omp)) { // Return with iflag set to initiate a quadratic iteration. *pbIflag = true ; return ; @@ -565,7 +565,7 @@ Rpoly::realit( double sss, int* nz, bool* pbIflag) kv = kv*s + k[i] ; qk[i] = kv; } - if ( fabs( kv) <= fabs( k[n-1]) * 10.0 * eta) { + if ( abs( kv) <= abs( k[n-1]) * 10.0 * eta) { // Use unscaled form. k[0] = 0.0 ; for ( i = 1 ; i < n ; i ++) { @@ -585,7 +585,7 @@ Rpoly::realit( double sss, int* nz, bool* pbIflag) kv = kv * s + k[i] ; } t = 0.0 ; - if ( fabs( kv) > ( fabs( k[n-1] * 10.0 * eta))) + if ( abs( kv) > ( abs( k[n-1] * 10.0 * eta))) t = - pv / kv ; s += t ; } @@ -605,14 +605,14 @@ Rpoly::calcsc( int *type) quadsd( n-1, &u, &v, k, qk, &c, &d) ; // Type=3 indicates the quadratic is almost a factor of k. - if ( fabs( c) <= fabs( k[n-1] * 100.0 * eta) && - fabs( d) <= fabs( k[n-2] * 100.0 * eta)) { + if ( abs( c) <= abs( k[n-1] * 100.0 * eta) && + abs( d) <= abs( k[n-2] * 100.0 * eta)) { *type = 3 ; return ; } // Type=1 indicates that all formulas are divided by c. - if ( fabs( d) < fabs( c)) { + if ( abs( d) < abs( c)) { *type = 1 ; e = a / c ; f = d / c ; @@ -657,7 +657,7 @@ Rpoly::nextk( int* type) temp = a ; if ( *type == 1) temp = b ; - if ( fabs( a1) <= fabs( temp) * eta * 10.0) { + if ( abs( a1) <= abs( temp) * eta * 10.0) { // If a1 is nearly zero then use a special form of the recurrence. k[0] = 0.0; k[1] = -a7*qp[0] ; @@ -774,22 +774,22 @@ Rpoly::quad( double a, double b1, double c, } /* Compute discriminant avoiding overflow. */ b = b1 / 2.0 ; - if ( fabs( b) < fabs( c)) { + if ( abs( b) < abs( c)) { if ( c < 0.0) e = - a ; else e = a ; - e = b * ( b / fabs( c)) - e ; - d = sqrt( fabs( e)) * sqrt( fabs( c)) ; + e = b * ( b / abs( c)) - e ; + d = sqrt( abs( e)) * sqrt( abs( c)) ; } else { e = 1.0 - ( a / b) *( c / b) ; - d = sqrt( fabs( e)) * fabs( b) ; + d = sqrt( abs( e)) * abs( b) ; } if ( e < 0.0) { /* complex conjugate zeros */ *sr = - b / a ; *lr = *sr ; - *si = fabs( d / a) ; + *si = abs( d / a) ; *li = - ( *si) ; } else { @@ -1266,7 +1266,7 @@ Cpoly::cauchy( const int nn, double pt[], double q[], double* fn_val) dx = x ; // Do Newton iteration until x converges to two decimal places - while ( fabs( dx / x ) > 0.005) { + while ( abs( dx / x ) > 0.005) { q[0] = pt[0] ; for ( i = 1 ; i <= nn ; i++) q[i] = q[i - 1] * x + pt[i] ; @@ -1346,7 +1346,7 @@ Cpoly::cdivid( const double ar, const double ai, const double br, const double b return ; } - if ( fabs( br) < fabs( bi)) { + if ( abs( br) < abs( bi)) { r = br / bi ; dinv = 1.0 / ( bi + r * br) ; *cr = ( ar * r + ai) * dinv ; @@ -1366,11 +1366,8 @@ Cpoly::cdivid( const double ar, const double ai, const double br, const double b double Cpoly::cmod( const double r, const double i) { - double ar, ai ; - - - ar = fabs( r) ; - ai = fabs( i) ; + double ar = abs( r) ; + double ai = abs( i) ; if ( ar < ai) return ( ai * sqrt( 1.0 + ( ar * ar) / ( ai * ai))) ; diff --git a/Polynomial.cpp b/Polynomial.cpp index 574b239..cfefff8 100644 --- a/Polynomial.cpp +++ b/Polynomial.cpp @@ -171,7 +171,7 @@ void Polynomial::AdjustDegree( void) { // se il coefficiente del grado più alto è zero, diminuisco il grado - while ( m_nDegree >= 0 && fabs( m_Coeff[m_nDegree]) < DBL_EPSILON) { + while ( m_nDegree >= 0 && abs( m_Coeff[m_nDegree]) < DBL_EPSILON) { m_Coeff.pop_back() ; -- m_nDegree ; } @@ -206,7 +206,7 @@ FilterMultipleAndOutOfRangeRoots( DBLVECTOR& vRoots, double dMin, double dMax, d int nZ = (int) vRoots.size() ; for ( int i = 0 ; i < nZ ;) { if ( vRoots[i] < dMin || vRoots[i] > dMax || - ( i >= 1 && fabs( vRoots[i]- vRoots[i-1]) < dEps)) { + ( i >= 1 && abs( vRoots[i]- vRoots[i-1]) < dEps)) { nZ -- ; for ( int j = i ; j < nZ ; ++ j) vRoots[j] = vRoots[j+1] ; diff --git a/PolynomialRoots.cpp b/PolynomialRoots.cpp index d5bb26e..e9bfca1 100644 --- a/PolynomialRoots.cpp +++ b/PolynomialRoots.cpp @@ -1,7 +1,7 @@ //---------------------------------------------------------------------------- -// EgalTech 2013-2013 +// EgalTech 2014-2018 //---------------------------------------------------------------------------- -// File : PolynomialRoots.cpp Data : 08.01.14 Versione : 1.5a2 +// File : PolynomialRoots.cpp Data : 07.08.18 Versione : 1.9h1 // Contenuto : Funzioni per il calcolo degli zeri di polinomi. // // @@ -15,28 +15,16 @@ #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[]) ; +#include +using namespace std ; //---------------------------------------------------------------------------- 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) + if ( pnIter != nullptr) *pnIter = 0 ; // controllo che il vettore dei coefficienti sia lungo almeno come il grado + 1 @@ -44,7 +32,7 @@ PolynomialRoots( int nDegree, DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter) return 0 ; // se il coefficiente del grado più alto è zero, diminuisco il grado - while ( nDegree >= 0 && fabs( vdPoly[nDegree]) < DBL_EPSILON) + while ( nDegree >= 0 && abs( vdPoly[nDegree]) < DBL_EPSILON) nDegree -- ; // se il grado è nullo o negativo, errore @@ -56,27 +44,32 @@ PolynomialRoots( int nDegree, DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter) return 0 ; // riordino i coefficienti reali - for ( i = 0 ; i <= nDegree ; i++) + double dPreal[POLY_MAXDEG+1] ; + for ( int i = 0 ; i <= nDegree ; i++) dPreal[i] = vdPoly[nDegree-i] ; // calcolo gli zeri - nZeros = cRpoly.Calculate( dPreal, nDegree, dZreal, dZcplx) ; + 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 ( i = 0 ; i < nZeros ; i++) { - if ( fabs( dZcplx[i]) < 100 * DBL_EPSILON) { + 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 - SortRoots( nZeros, vdRoot.data()) ; + sort( vdRoot.begin(), vdRoot.end(), + []( const double& a, const double& b) { return ( a > b) ; }) ; // assegno il numero di iterazioni - if ( pnIter != NULL) + if ( pnIter != nullptr) *pnIter = cRpoly.itercnt ; return nZeros ; @@ -86,19 +79,8 @@ PolynomialRoots( int nDegree, DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter) 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) + if ( pnIter != nullptr) *pnIter = 0 ; // controllo che il vettore dei coefficienti sia lungo almeno come il grado + 1 @@ -106,7 +88,7 @@ PolynomialRoots( int nDegree, CPLXVECTOR& vcPoly, CPLXVECTOR& vcRoot, int* pnIte return 0 ; // se il coefficiente del grado più alto è zero, diminuisco il grado - while ( nDegree >= 0 && m2( vcPoly[nDegree]) < DBL_EPSILON * DBL_EPSILON) + while ( nDegree >= 0 && norm( vcPoly[nDegree]) < DBL_EPSILON * DBL_EPSILON) nDegree -- ; // se il grado è nullo o negativo, errore @@ -118,128 +100,62 @@ PolynomialRoots( int nDegree, CPLXVECTOR& vcPoly, CPLXVECTOR& vcRoot, int* pnIte return 0 ; // ricavo i coefficienti reali - for ( i = 0 ; i <= nDegree ; i++) - dPreal[i] = vcPoly[nDegree-i].re ; + double dPreal[POLY_MAXDEG+1] ; + for ( int i = 0 ; i <= nDegree ; i++) + dPreal[i] = real( vcPoly[nDegree-i]) ; - // 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 ; + // 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 - if ( bCplx) - nZeros = cCpoly.Calculate( dPreal, dPcplx, nDegree, dZreal, dZcplx) ; - else + 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 ( i = 0 ; i < nZeros ; i++) { - vcRoot.push_back( Complex( dZreal[i], dZcplx[i])) ; + for ( int 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 ; + 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 - SortRoots( nZeros, vcRoot.data()) ; + // ordino le radici in senso decrescente della parte reale e in subordine della immaginaria + sort( vcRoot.begin(), vcRoot.end(), + []( const complex& a, const complex& 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 != NULL) - *pnIter = ( bCplx ? cCpoly.itercnt : cRpoly.itercnt) ; + if ( pnIter != nullptr) + *pnIter = nIterCnt ; 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) ; -}