diff --git a/Complex.cpp b/Complex.cpp new file mode 100644 index 0000000..2652320 --- /dev/null +++ b/Complex.cpp @@ -0,0 +1,271 @@ +//---------------------------------------------------------------------------- +// 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/DllMain.h b/DllMain.h new file mode 100644 index 0000000..b2aa4b4 --- /dev/null +++ b/DllMain.h @@ -0,0 +1,20 @@ +//---------------------------------------------------------------------------- +// EgalTech 2013-2014 +//---------------------------------------------------------------------------- +// File : DllMain.h Data : 08.01.14 Versione : 1.5a1 +// Contenuto : Prototipi funzioni per uso locale della DLL. +// +// +// +// Modifiche : 08.01.14 DS Creazione modulo. +// +// +//---------------------------------------------------------------------------- + +#pragma once + + +#include "/EgtDev/Include/EgtILogger.h" + +//----------------------------------------------------------------------------- +ILogger* GetENkLogger( void) ; diff --git a/ENkDllMain.cpp b/ENkDllMain.cpp new file mode 100644 index 0000000..85c71a2 --- /dev/null +++ b/ENkDllMain.cpp @@ -0,0 +1,76 @@ +//---------------------------------------------------------------------------- +// EgalTech 2013-2013 +//---------------------------------------------------------------------------- +// File : ENkDllMain.cpp Data : 08.01.14 Versione : 1.5a1 +// Contenuto : Inizializzazione della DLL. +// +// +// +// Modifiche : 08.01.14 DS Creazione modulo. +// +// +//---------------------------------------------------------------------------- + +//--------------------------- Include ---------------------------------------- +#include "stdafx.h" +#include "\EgtDev\Include\ENkDllMain.h" +#include "\EgtDev\Include\EgnGetModuleVer.h" +#include "\EgtDev\Include\EgtTrace.h" + +//--------------------------- Costanti ---------------------------------------- +#if defined( _DEBUG) + const char* ENK_STR = "EgtNumKernelD32.dll ver. " ; +#else + const char* ENK_STR = "EgtNumKernelR32.dll ver. " ; +#endif +const int STR_DIM = 40 ; + +//----------------------------------------------------------------------------- +static HINSTANCE s_hModule = NULL ; +static char s_szENkNameVer[STR_DIM] ; + +//----------------------------------------------------------------------------- +extern "C" int APIENTRY +DllMain( HMODULE hModule, DWORD dwReason, LPVOID lpReserved) +{ + + if ( dwReason == DLL_PROCESS_ATTACH) { + s_hModule = hModule ; + EGT_TRACE( "EgtNumKernel.dll Initializing!\n") ; + } + else if ( dwReason == DLL_PROCESS_DETACH) { + s_hModule = NULL ; + EGT_TRACE( "EgtNumKernel.dll Terminating!\n") ; + } + + return 1 ; +} + +//----------------------------------------------------------------------------- +const char* +GetENkVersion( void) +{ + std::string sVer ; + + GetModuleVersion( s_hModule, sVer) ; + sprintf_s( s_szENkNameVer, STR_DIM, "%s%s", ENK_STR, sVer.c_str()) ; + + return s_szENkNameVer ; +} + +//----------------------------------------------------------------------------- +static ILogger* s_pLogger = nullptr ; + +//----------------------------------------------------------------------------- +void +SetENkLogger( ILogger* pLogger) +{ + s_pLogger = pLogger ; +} + +//----------------------------------------------------------------------------- +ILogger* +GetENkLogger( void) +{ + return s_pLogger ; +} diff --git a/EgtNumKernel.rc b/EgtNumKernel.rc new file mode 100644 index 0000000..a8d3075 Binary files /dev/null and b/EgtNumKernel.rc differ diff --git a/EgtNumKernel.sln b/EgtNumKernel.sln new file mode 100644 index 0000000..dec073f --- /dev/null +++ b/EgtNumKernel.sln @@ -0,0 +1,20 @@ + +Microsoft Visual Studio Solution File, Format Version 11.00 +# Visual Studio 2010 +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "EgtNumKernel", "EgtNumKernel.vcxproj", "{E47BBFD0-36EA-4EC1-9D6D-B05CA9C092C6}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|Win32 = Debug|Win32 + Release|Win32 = Release|Win32 + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {E47BBFD0-36EA-4EC1-9D6D-B05CA9C092C6}.Debug|Win32.ActiveCfg = Debug|Win32 + {E47BBFD0-36EA-4EC1-9D6D-B05CA9C092C6}.Debug|Win32.Build.0 = Debug|Win32 + {E47BBFD0-36EA-4EC1-9D6D-B05CA9C092C6}.Release|Win32.ActiveCfg = Release|Win32 + {E47BBFD0-36EA-4EC1-9D6D-B05CA9C092C6}.Release|Win32.Build.0 = Release|Win32 + EndGlobalSection + GlobalSection(SolutionProperties) = preSolution + HideSolutionNode = FALSE + EndGlobalSection +EndGlobal diff --git a/EgtNumKernel.vcxproj b/EgtNumKernel.vcxproj new file mode 100644 index 0000000..c844e3f --- /dev/null +++ b/EgtNumKernel.vcxproj @@ -0,0 +1,125 @@ + + + + + Debug + Win32 + + + Release + Win32 + + + + {E47BBFD0-36EA-4EC1-9D6D-B05CA9C092C6} + Win32Proj + EgtNumKernel + + + + DynamicLibrary + true + Unicode + + + DynamicLibrary + false + true + Unicode + + + + + + + + + + + + + true + $(SolutionDir)$(Configuration)$(PlatformArchitecture)\ + $(Configuration)$(PlatformArchitecture)\ + $(ProjectName)D$(PlatformArchitecture) + + + false + $(SolutionDir)$(Configuration)$(PlatformArchitecture)\ + $(Configuration)$(PlatformArchitecture)\ + $(ProjectName)R$(PlatformArchitecture) + + + + Use + Level3 + Disabled + WIN32;_DEBUG;_WINDOWS;_USRDLL;I_AM_ENK;%(PreprocessorDefinitions) + true + + + Windows + true + + + copy $(TargetDir)$(TargetName).pdb \EgtDev\Lib\ +copy $(TargetDir)$(TargetName).lib \EgtDev\Lib\ +copy $(TargetPath) \EgtProg\Dll + + + _UNICODE;UNICODE;_DEBUG;%(PreprocessorDefinitions) + + + + + Level3 + Use + MaxSpeed + true + true + WIN32;NDEBUG;_WINDOWS;_USRDLL;I_AM_ENK;%(PreprocessorDefinitions) + true + + + Windows + false + true + true + + + copy $(TargetDir)$(TargetName).pdb \EgtDev\Lib\ +copy $(TargetDir)$(TargetName).lib \EgtDev\Lib\ +copy $(TargetPath) \EgtProg\Dll + + + _UNICODE;UNICODE;NDEBUG;%(PreprocessorDefinitions) + + + + + + + + + Create + Create + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/EgtNumKernel.vcxproj.filters b/EgtNumKernel.vcxproj.filters new file mode 100644 index 0000000..d516aaa --- /dev/null +++ b/EgtNumKernel.vcxproj.filters @@ -0,0 +1,68 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hpp;hxx;hm;inl;inc;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms + + + + + File di origine + + + File di origine + + + File di origine + + + File di origine + + + File di origine + + + + + File di intestazione + + + File di intestazione + + + File di intestazione + + + File di intestazione + + + File di intestazione + + + File di intestazione + + + File di intestazione + + + File di intestazione + + + File di intestazione + + + + + File di risorse + + + \ No newline at end of file diff --git a/JenkinsTraub.cpp b/JenkinsTraub.cpp new file mode 100644 index 0000000..32c95b8 --- /dev/null +++ b/JenkinsTraub.cpp @@ -0,0 +1,1390 @@ +//---------------------------------------------------------------------------- +// EgalTech 2013-2014 +//---------------------------------------------------------------------------- +// File : JenkinsTraub.cpp Data : 08.01.14 Versione : 1.5a1 +// Contenuto : Implementazione calcolo degli zeri di polinomi a coefficienti +// reali o complessi con il metodo di Jenkins e Traub. +// Rpoly deriva da TOMS493. Cpoly deriva da TOMS419. +// +// +// Modifiche : 08.01.14 DS Creazione modulo. +// +// +//---------------------------------------------------------------------------- + +//--------------------------- Include ---------------------------------------- +#include "stdafx.h" +#include "JenkinsTraub.h" + + +//--------------------------- Class Rpoly -------------------------------------- +//------------------------------------------------------------------------------ +// IN: op - double precision vector of coefficients in order of decreasing powers. +// degree - integer degree of polynomial +// OUT: zeror,zeroi - output double precision vectors of the real and imaginary parts of the zeros. +// RET: -1 if leading coefficient is zero, otherwise number of roots found. +//------------------------------------------------------------------------------ +int +Rpoly::Calculate( const double* op, int degree, double* zeror, double* zeroi) +{ + bool bZerOk ; + bool bScale ; + bool bContinue ; + int cnt, nz, i, j, jj, l, nm1 ; + double t, aa, bb, cc, factor, rot ; + double lo, max, min, xx, yy, cosr, sinr, xxx, x, sc, bnd ; + double xm, ff, df, dx, infin, smalno, base ; + double temp[POLY_MAXDEG+1] ; + double pt[POLY_MAXDEG+1] ; + + + // The following statements set machine constants. + base = _DBL_RADIX ; + eta = DBL_EPSILON ; + infin = DBL_MAX ; + smalno = DBL_MIN ; + are = eta ; + mre = eta ; + lo = smalno / eta ; + + // Initialization of constants for shift rotation. + xx = sqrt( 0.5) ; + yy = - xx ; + rot = 94.0 ; + rot *= 0.017453293 ; + cosr = cos( rot) ; + sinr = sin( rot) ; + n = degree ; + + // Inizializzo numero iterazioni + itercnt = 0 ; + + // Algorithm fails if the leading coefficient is zero. + if ( op[0] == 0.0) + return -1 ; + + // Remove the zeros at the origin, if any. + while ( op[n] == 0.0) { + j = degree - n ; + zeror[j] = 0.0 ; + zeroi[j] = 0.0 ; + n -- ; + } + if ( n < 1) + return degree ; + + // Make a copy of the coefficients. + for ( i = 0 ; i <= n ; i++) + p[i] = op[i] ; + + // Start the algorithm for one zero. + bContinue = true ; + while ( bContinue) { + + // Calculate the final zero + if ( n == 1) { + zeror[degree-1] = - p[1] / p[0] ; + zeroi[degree-1] = 0.0 ; + n -= 1 ; + return ( degree - n) ; + } + + // Calculate a pair of zeros. + if ( n == 2) { + quad( p[0], p[1], p[2], &zeror[degree-2], &zeroi[degree-2], + &zeror[degree-1], &zeroi[degree-1]) ; + n -= 2 ; + return ( degree - n) ; + } + + // Find largest and smallest moduli of coefficients. + max = 0.0 ; + min = infin ; + for ( i = 0 ; i <= n ; i++) { + x = fabs( p[i]) ; + if ( x > max) + max = x ; + if ( x != 0.0 && x < min) + min = x ; + } + // Scale if there are large or very small coefficients. + // Computes a scale factor to multiply the coefficients of the + // polynomial. The scaling si done to avoid overflow and to + // avoid undetected underflow interfering with the convergence + // criterion. The factor is a power of the base. + bScale = true ; + sc = lo / min ; + if ( sc > 1.0 && ( infin / sc) < max) + bScale = false ; + if ( sc <= 1.0) { + if ( max < 10.0) + bScale = false ; + if ( sc == 0.0) + sc = smalno ; + } + if ( bScale) { + // Scale polynomial. + l = (int)( log( sc) / log( base) + 0.5) ; + factor = pow( base * 1.0, l) ; + if ( factor != 1.0) { + for ( i = 0 ; i <= n ; i ++) + p[i] = factor * p[i] ; + } + } + + // Compute lower bound on moduli of roots. + for ( i = 0 ; i <= n ; i++) { + pt[i] = fabs( p[i]) ; + } + pt[n] = - pt[n] ; + // Compute upper estimate of bound. + x = exp( ( log( - pt[n]) - log( pt[0])) / (double) n) ; + // If Newton step at the origin is better, use it. + if ( pt[n-1] != 0.0) { + xm = - pt[n] / pt[n-1] ; + if ( xm < x) + x = xm ; + } + + // Chop the interval (0,x) until ff <= 0 + while ( true) { + xm = x * 0.1 ; + ff = pt[0] ; + for ( i = 1 ; i <= n ; i ++) + ff = ff * xm + pt[i] ; + if ( ff <= 0.0) + break ; + x = xm ; + } + + // Do Newton iteration until x converges to two decimal places. + dx = x ; + while ( fabs( dx / x) > 0.005) { + ff = pt[0] ; + df = ff ; + for ( i = 1 ; i < n ; i++) { + ff = ff * x + pt[i] ; + df = df * x + ff ; + } + ff = ff * x + pt[n] ; + dx = ff / df ; + x -= dx ; + itercnt ++ ; + } + bnd = x ; + + // Compute the derivative as the initial k polynomial and do 5 steps with no shift. + nm1 = n - 1 ; + for ( i = 1 ; i < n ; i++) + k[i] = (double)( n - i) * p[i] / (double) n ; + k[0] = p[0] ; + aa = p[n] ; + bb = p[n-1] ; + bZerOk = ( k[n-1] == 0) ; + for ( jj = 0 ; jj < 5 ; jj ++) { + itercnt ++ ; + cc = k[n-1] ; + // Use a scaled form of recurrence if value of k at 0 is nonzero. + if ( ! bZerOk) { + t = - aa / cc ; + for ( i = 0 ; i < nm1 ; i ++) { + j = n - i - 1 ; + k[j] = t * k[j-1] + p[j] ; + } + k[0] = p[0] ; + bZerOk = ( fabs( k[n-1]) <= fabs( bb) * eta * 10.0) ; + } + else { + // Use unscaled form of recurrence. + for ( i = 0 ; i < nm1 ; i ++) { + j = n - i - 1 ; + k[j] = k[j-1] ; + } + k[0] = 0.0 ; + bZerOk = ( k[n-1] == 0.0) ; + } + } + + // Save k for restarts with new shifts. + for ( i = 0 ; i < n ; i ++) + temp[i] = k[i] ; + + // Loop to select the quadratic corresponding to each new shift. + bContinue = false ; + for ( cnt = 0 ; cnt < 20 ; cnt ++) { + // Quadratic corresponds to a double shift to a non-real point and its complex conjugate. + // The point has modulus bnd and amplitude rotated by 94 degrees from the previous shift. + xxx = cosr * xx - sinr * yy ; + yy = sinr * xx + cosr * yy ; + xx = xxx ; + sr = bnd * xx ; + si = bnd * yy ; + u = -2.0 * sr ; + v = bnd ; + fxshfr( 20 * ( cnt + 1), &nz); + // The second stage jumps directly to one of the third stage iterations and returns here if successful. + // Deflate the polynomial, store the zero or zeros and return to the main algorithm. + if ( nz != 0) { + j = degree - n ; + zeror[j] = szr ; + zeroi[j] = szi ; + n -= nz ; + for ( i = 0 ; i <= n ; i ++) + p[i] = qp[i] ; + if ( nz != 1) { + zeror[j+1] = lzr ; + zeroi[j+1] = lzi ; + } + // interrompo il loop e riparto a cercare un nuovo zero + bContinue = true ; + break ; + } + // If the iteration is unsuccessful another quadratic is chosen after restoring k. + else { + for ( i = 0 ; i < n ; i ++) + k[i] = temp[i] ; + } + } + } + + // Return with failure if no convergence after 20 shifts. + return ( degree - n) ; +} + +//------------------------------------------------------------------------------ +// Computes up to L2 fixed shift k-polynomials, +// testing for convergence in the linear or quadratic +// case. Initiates one of the variable shift +// iterations and returns with the number of zeros found. +//------------------------------------------------------------------------------ +void +Rpoly::fxshfr( int l2, int* nz) +{ + bool bVpass ; + bool bSpass ; + bool bVtry ; + bool bStry ; + bool bIflag ; + int type, i, j ; + double svu, svv, ui, vi, s ; + double betas, betav, oss, ovv, ss, vv, ts, tv ; + double ots, otv, tvv, tss ; + double svk[POLY_MAXDEG+1] ; + + + // Inizializzazioni + *nz = 0 ; + betav = 0.25 ; + betas = 0.25 ; + oss = sr ; + ovv = v ; + + // Evaluate polynomial by synthetic division. + quadsd( n, &u, &v, p, qp, &a, &b) ; + calcsc( &type) ; + + for ( j = 0 ; j < l2 ; j ++) { + + // Calculate next k polynomial and estimate v. + nextk( &type) ; + calcsc( &type) ; + newest( type, &ui, &vi) ; + vv = vi ; + + // Estimate s. + ss = 0.0 ; + if ( k[n-1] != 0.0) + ss = - p[n] / k[n-1] ; + tv = 1.0 ; + ts = 1.0 ; + if ( j == 0 || type == 3) { + ovv = vv ; + oss = ss ; + otv = tv ; + ots = ts ; + continue ; + } + + // Compute relative measures of convergence of s and v sequences. + if ( vv != 0.0) + tv = fabs( ( vv - ovv) / vv) ; + if ( ss != 0.0) + ts = fabs( ( ss - oss) / ss) ; + /* If decreasing, multiply two most recent convergence measures. */ + tvv = 1.0 ; + if ( tv < otv) + tvv = tv * otv ; + tss = 1.0 ; + if ( ts < ots) + tss = ts * ots ; + // Compare with convergence criteria. + bVpass = ( tvv < betav) ; + bSpass = ( tss < betas) ; + if ( ! ( bSpass || bVpass)) { + ovv = vv ; + oss = ss ; + otv = tv ; + ots = ts ; + continue ; + } + + // At least one sequence has passed the convergence test. Store variables before iterating. + svu = u; + svv = v; + for ( i = 0 ; i < n ; i ++) + svk[i] = k[i] ; + s = ss ; + + // Choose iteration according to the fastest converging sequence. + bVtry = false ; + bStry = false ; + while ( true) { + bIflag = true ; + if ( bSpass && ! bVpass || tss < tvv) + ; + else { + quadit( &ui, &vi, nz) ; + if ( *nz > 0) + return ; + // Quadratic iteration has failed. Flag that it has been tried and decrease the convergence criterion. + bVtry = true ; + betav *= 0.25; + // Try linear iteration if it has not been tried and the S sequence is converging. + if ( bStry || ! bSpass) + bIflag = false ; + else { + for ( i = 0 ; i < n ; i ++) + k[i] = svk[i] ; + } + } + if ( bIflag) { + realit( s, nz, &bIflag) ; + if ( *nz > 0) + return ; + // Linear iteration has failed. Flag that it has been tried and decrease the convergence criterion. + bStry = true ; + betas *= 0.25 ; + // If linear iteration signals an almost double real zero attempt quadratic iteration. + if ( bIflag) { + ui = -( s + s) ; + vi = s * s ; + break ; + } + } + // Restore variables + u = svu; + v = svv; + for ( i = 0 ; i < n ; i ++) { + k[i] = svk[i] ; + } + // Try quadratic iteration if it has not been tried and the V sequence is convergin. + if ( ! bVpass || bVtry) + break ; + } + + // Recompute QP and scalar values to continue the second stage. + quadsd( n, &u, &v, p, qp, &a, &b) ; + calcsc( &type) ; + + // Salvo valori come precedenti + ovv = vv ; + oss = ss ; + otv = tv ; + ots = ts ; + } +} + +//------------------------------------------------------------------------------ +// Variable-shift k-polynomial iteration for a +// quadratic factor converges only if the zeros are +// equimodular or nearly so. +// uu, vv - coefficients of starting quadratic. +// nz - number of zeros found. +//------------------------------------------------------------------------------ +void +Rpoly::quadit( double *uu, double *vv, int *nz) +{ + bool bTried ; + int type, i, j ; + double ui, vi ; + double mp, omp, ee, relstp, t, zm ; + + + // Inizializzazioni + *nz = 0 ; + bTried = false ; + u = *uu ; + v = *vv ; + j = 0 ; + + // Main loop. + while ( true) { + itercnt ++ ; + + quad( 1.0, u, v, &szr, &szi, &lzr, &lzi) ; + + // 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)) + return ; + + // Evaluate polynomial by quadratic synthetic division. + quadsd( n, &u, &v, p, qp, &a, &b) ; + mp = fabs( a - szr * b) + fabs( szi * b) ; + // Compute a rigorous bound on the rounding error in evaluating p. + zm = sqrt( fabs( v)) ; + ee = 2.0 * fabs( qp[0]) ; + t = -szr * b ; + for ( i = 1 ; i < n ; i ++) { + ee = ee * zm + fabs( qp[i]) ; + } + ee = ee * zm + fabs( 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) ; + // Iteration has converged sufficiently if the polynomial value is less than 20 times this bound. + if ( mp <= 20.0 * ee) { + *nz = 2 ; + return ; + } + j ++ ; + + // Stop iteration after 20 steps. + if ( j > 20) + return ; + + // A cluster appears to be stalling the convergence. + // Five fixed shift steps are taken with a u,v close to the cluster. + if ( j >= 2 && + ! ( relstp > 0.01 || mp < omp || bTried)) { + if ( relstp < eta) + relstp = eta; + relstp = sqrt( relstp) ; + u = u - u * relstp ; + v = v + v * relstp ; + quadsd( n, &u, &v, p, qp, &a, &b) ; + for ( i = 0 ; i < 5 ; i ++) { + calcsc( &type) ; + nextk( &type) ; + } + bTried = true ; + j = 0 ; + } + + // Salvo valore + omp = mp ; + + // Calculate next k polynomial and new u and v. + calcsc( &type) ; + nextk( &type) ; + calcsc( &type) ; + newest( type, &ui, &vi) ; + // If vi is zero the iteration is not converging. + if ( vi == 0.0) + return ; + relstp = fabs( ( vi - v) / vi) ; + u = ui ; + v = vi ; + } +} + +//------------------------------------------------------------------------------ +// Variable-shift H polynomial iteration for a real zero. +// sss - starting iterate +// nz - number of zeros found +// iflag - flag to indicate a pair of zeros near real axis. +//------------------------------------------------------------------------------ +void +Rpoly::realit( double sss, int* nz, bool* pbIflag) +{ + int i, j ; + double pv, kv, t, s ; + double ms, mp, omp, ee ; + + + // Inizializzazioni + *nz = 0 ; + s = sss ; + *pbIflag = false ; + j = 0 ; + + // Main loop + while ( true) { + itercnt ++ ; + pv = p[0] ; + // Evaluate p at s. + qp[0] = pv ; + for ( i = 1 ; i <= n ; i++) { + pv = pv * s + p[i] ; + qp[i] = pv ; + } + mp = fabs( pv) ; + // Compute a rigorous bound on the error in evaluating p. + ms = fabs( s) ; + ee = ( mre / ( are + mre)) * fabs( qp[0]) ; + for ( i = 1 ; i <= n ; i ++) { + ee = ee * ms + fabs( 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)) { + *nz = 1 ; + szr = s ; + szi = 0.0 ; + return ; + } + j ++ ; + // Stop iteration after 10 steps. + if ( j > 10) + return ; + // A cluster of zeros near the real axis has been encountered. + if ( j >= 2 && + ! ( fabs( t) > 0.001 * fabs( s-t) || mp < omp)) { + // Return with iflag set to initiate a quadratic iteration. + *pbIflag = true ; + return ; + } + + // Return if the polynomial value has increased significantly. + omp = mp ; + + // Compute t, the next polynomial, and the new iterate. + kv = k[0] ; + qk[0] = kv ; + for ( i = 1 ; i < n ; i ++) { + kv = kv*s + k[i] ; + qk[i] = kv; + } + if ( fabs( kv) <= fabs( k[n-1]) * 10.0 * eta) { + // Use unscaled form. + k[0] = 0.0 ; + for ( i = 1 ; i < n ; i ++) { + k[i] = qk[i-1] ; + } + } + else { + // Use the scaled form of the recurrence if the value of k at s is nonzero. + t = - pv / kv ; + k[0] = qp[0] ; + for ( i = 1 ; i < n ; i++) { + k[i] = t * qk[i-1] + qp[i] ; + } + } + kv = k[0] ; + for ( i = 1 ; i < n ; i ++) { + kv = kv * s + k[i] ; + } + t = 0.0 ; + if ( fabs( kv) > ( fabs( k[n-1] * 10.0 * eta))) + t = - pv / kv ; + s += t ; + } +} + +//------------------------------------------------------------------------------ +// This routine calculates scalar quantities used to +// compute the next k polynomial and new estimates of +// the quadratic coefficients. +// type - integer variable set here indicating how the +// calculations are normalized to avoid overflow. +//------------------------------------------------------------------------------ +void +Rpoly::calcsc( int *type) +{ + // Synthetic division of k by the quadratic 1,u,v + 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)) { + *type = 3 ; + return ; + } + + // Type=1 indicates that all formulas are divided by c. + if ( fabs( d) < fabs( c)) { + *type = 1 ; + e = a / c ; + f = d / c ; + g = u * e ; + h = v * b ; + a3 = a * e + ( h / c + g) * b ; + a1 = b - a * ( d / c) ; + a7 = a + g * d + h * f ; + return ; + } + + // Type=2 indicates that all formulas are divided by d. + *type = 2 ; + e = a / d ; + f = c / d ; + g = u * b ; + h = v * b ; + a3 = ( a + g) * e + h * ( b / d) ; + a1 = b * f - a ; + a7 = ( f + u)*a + h ; +} + +//------------------------------------------------------------------------------ +// Computes the next k polynomials using scalars computed in calcsc. +//------------------------------------------------------------------------------ +void +Rpoly::nextk( int* type) +{ + double temp ; + int i ; + + + if ( *type == 3) { + /* Use unscaled form of the recurrence if type is 3. */ + k[0] = 0.0 ; + k[1] = 0.0 ; + for ( i = 2 ; i < n ; i ++) { + k[i] = qk[i-2] ; + } + return ; + } + temp = a ; + if ( *type == 1) + temp = b ; + if ( fabs( a1) <= fabs( 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] ; + for ( i = 2 ; i < n ; i ++) { + k[i] = a3 * qk[i-2] - a7 * qp[i-1] ; + } + return ; // HVE return added + } + /* Use scaled form of the recurrence. */ + a7 /= a1 ; + a3 /= a1 ; + k[0] = qp[0] ; + k[1] = qp[1] - a7 * qp[0] ; + for ( i = 2 ; i < n ; i ++) { + k[i] = a3 * qk[i-2] - a7 * qp[i-1] + qp[i] ; + } +} + +//------------------------------------------------------------------------------ +// Compute new estimates of the quadratic coefficients using the scalars computed in calcsc. +//------------------------------------------------------------------------------ +void +Rpoly::newest( int type, double *uu, double *vv) +{ + double a4,a5,b1,b2,c1,c2,c3,c4,temp; + + + /* Use formulas appropriate to setting of type. */ + if ( type == 3) { + /* If type=3 the quadratic is zeroed. */ + *uu = 0.0 ; + *vv = 0.0 ; + return ; + } + if ( type == 2) { + a4 = ( a + g) * f + h ; + a5 = ( f + u) * c + v * d ; + } + else { + a4 = a + u * b + h * f ; + a5 = c + ( u + v * f) * d ; + } + /* Evaluate new quadratic coefficients. */ + b1 = -k[n-1] / p[n] ; + b2 = -( k[n-2] + b1 * p[n-1]) / p[n] ; + c1 = v * b2 * a1 ; + c2 = b1 * a7 ; + c3 = b1 * b1 * a3 ; + c4 = c1 - c2 - c3 ; + temp = a5 + b1 * a4 - c4 ; + if ( temp == 0.0) { + *uu = 0.0 ; + *vv = 0.0 ; + return ; + } + *uu = u - ( u * ( c3 + c2) + v * ( b1 * a1 + b2 * a7)) / temp ; + *vv = v * ( 1.0 + c4 / temp) ; + return ; +} + +//------------------------------------------------------------------------------ +// Divides p by the quadratic 1,u,v placing the quotient in q and the remainder in a,b. +//------------------------------------------------------------------------------ +void +Rpoly::quadsd( int nn, double *u, double *v, double *p, double *q, + double *a, double *b) +{ + int i ; + double c ; + + + *b = p[0] ; + q[0] = *b ; + *a = p[1] - (*b) * (*u) ; + q[1] = *a ; + for ( i = 2 ; i <= nn ; i++) { + c = p[i] - (*a) * (*u) - (*b) * (*v) ; + q[i] = c ; + *b = *a ; + *a = c ; + } +} + +//------------------------------------------------------------------------------ +// Calculate the zeros of the quadratic a*z^2 + b1*z + c. +// The quadratic formula, modified to avoid overflow, is used +// to find the larger zero if the zeros are real and both +// are complex. The smaller real zero is found directly from +// the product of the zeros c/a. +//------------------------------------------------------------------------------ +void +Rpoly::quad( double a, double b1, double c, + double* sr, double* si, double* lr, double* li) +{ + double b, d, e ; + + + if ( a == 0.0) { /* less than two roots */ + if ( b1 != 0.0) + *sr = - c / b1 ; + else + *sr = 0.0 ; + *lr = 0.0 ; + *si = 0.0 ; + *li = 0.0 ; + return; + } + if ( c == 0.0) { /* one real root, one zero root */ + *sr = 0.0 ; + *lr = - b1 / a ; + *si = 0.0 ; + *li = 0.0 ; + return; + } + /* Compute discriminant avoiding overflow. */ + b = b1 / 2.0 ; + if ( fabs( b) < fabs( c)) { + if ( c < 0.0) + e = - a ; + else + e = a ; + e = b * ( b / fabs( c)) - e ; + d = sqrt( fabs( e)) * sqrt( fabs( c)) ; + } + else { + e = 1.0 - ( a / b) *( c / b) ; + d = sqrt( fabs( e)) * fabs( b) ; + } + if ( e < 0.0) { /* complex conjugate zeros */ + *sr = - b / a ; + *lr = *sr ; + *si = fabs( d / a) ; + *li = - ( *si) ; + } + else { + if ( b >= 0.0) + d = - d ; /* real zeros. */ + *lr = ( - b + d) / a ; + *sr = 0.0 ; + if ( *lr != 0.0) + *sr = ( c / *lr) / a ; + *si = 0.0 ; + *li = 0.0 ; + } +} + + +//--------------------------- Class Cpoly -------------------------------------- +//------------------------------------------------------------------------------ +// IN: opr, opi - double precision vector of real and imaginary coefficients in order of decreasing powers. +// degree - integer degree of polynomial +// OUT: zeror,zeroi - output double precision vectors of the real and imaginary parts of the zeros. +// RET: -1 if leading coefficient is zero, otherwise number of roots found. +//------------------------------------------------------------------------------ +int +Cpoly::Calculate( const double* opr, const double* opi, int degree, double* zeror, double* zeroi) +{ + bool bContinue ; + bool bConv ; + int cnt1, cnt2, idnn2, i ; + double xx, yy, cosr, sinr, smalno, base, xxx, zr, zi, bnd ; + + + // The following statements set machine constants. + mcon( &eta, &infin, &smalno, &base) ; + are = eta ; + mre = 2.0 * sqrt( 2.0 ) * eta ; + + // Initialization of constants for shift rotation. + xx = 0.70710678 ; + yy = -xx ; + cosr = -0.060756474 ; + sinr = -0.99756405 ; + nn = degree ; + + // Inizializzo numero iterazioni + itercnt = 0 ; + + // Algorithm fails if the leading coefficient is zero + if ( opr[0] == 0 && opi[0] == 0) + return - 1 ; + + // Remove the zeros at the origin if any + while ( opr[nn] == 0 && opi[nn] == 0) { + idnn2 = degree - nn ; + zeror[idnn2] = 0 ; + zeroi[idnn2] = 0 ; + nn -- ; + } + + // Make a copy of the coefficients + for ( i = 0 ; i <= nn ; i++) { + pr[i] = opr[i] ; + pi[i] = opi[i] ; + shr[i] = cmod( pr[i], pi[i]) ; + } + + // Scale the polynomial + bnd = scale( nn, shr, eta, infin, smalno, base) ; + if ( bnd != 1) + for ( i = 0 ; i <= nn ; i++) { + pr[i] *= bnd ; + pi[i] *= bnd ; + } + + // Main loop + bContinue = true ; + while ( bContinue) { + + if ( nn <= 1) { + cdivid( -pr[1], -pi[1], pr[0], pi[0], &zeror[degree-1], &zeroi[degree-1]) ; + return degree ; + } + + // Calculate bnd, alower bound on the modulus of the zeros + for ( i = 0 ; i <= nn ; i++) + shr[i] = cmod( pr[i], pi[i]) ; + + cauchy( nn, shr, shi, &bnd) ; + + // Outer loop to control 2 Major passes with different sequences of shifts + bContinue = false ; + for ( cnt1 = 1 ; cnt1 <= 2 && ! bContinue ; cnt1++) { + // First stage calculation , no shift + noshft( 5) ; + + // Inner loop to select a shift + for ( cnt2 = 1 ; cnt2 <= 9 && ! bContinue ; cnt2++) { + // Shift is chosen with modulus bnd and amplitude rotated by 94 degree from the previous shif + xxx = cosr * xx - sinr * yy ; + yy = sinr * xx + cosr * yy ; + xx = xxx ; + sr = bnd * xx ; + si = bnd * yy ; + + // Second stage calculation, fixed shift + fxshft( 10 * cnt2, &zr, &zi, &bConv) ; + if ( bConv) { + // The second stage jumps directly to the third stage ieration + // If successful the zero is stored and the polynomial deflated + idnn2 = degree - nn ; + zeror[idnn2] = zr ; + zeroi[idnn2] = zi ; + nn -- ; + for ( i = 0 ; i <= nn ; i++) { + pr[i] = qpr[i] ; + pi[i] = qpi[i] ; + } + bContinue = true ; + } + // If the iteration is unsuccessful another shift is chosen + } + // if 9 shifts fail, the outer loop is repeated with another sequence of shifts + } + } + + // The zerofinder has failed on two major passes + // return empty handed with the number of roots found (less than the original degree) + degree -= nn ; + + return degree ; +} + +//------------------------------------------------------------------------------ +// COMPUTES THE DERIVATIVE POLYNOMIAL AS THE INITIAL H +// POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS. +//------------------------------------------------------------------------------ +void +Cpoly::noshft( const int l1) +{ + int i, j, jj, n, nm1 ; + double xni, t1, t2 ; + + + n = nn ; + nm1 = n - 1 ; + for ( i = 0 ; i < n ; i++) { + xni = nn - i ; + hr[i] = xni * pr[i] / n ; + hi[i] = xni * pi[i] / n ; + } + for ( jj = 1 ; jj <= l1 ; jj++) { + itercnt ++ ; + if ( cmod( hr[n - 1], hi[n - 1]) > eta * 10 * cmod( pr[n - 1], pi[n - 1])) { + cdivid( -pr[nn], -pi[nn], hr[n - 1], hi[n - 1], &tr, &ti) ; + for ( i = 0 ; i < nm1 ; i++) { + j = nn - i - 1 ; + t1 = hr[j - 1] ; + t2 = hi[j - 1] ; + hr[j] = tr * t1 - ti * t2 + pr[j] ; + hi[j] = tr * t2 + ti * t1 + pi[j] ; + } + hr[0] = pr[0] ; + hi[0] = pi[0] ; + } + else { + // If the constant term is essentially zero, shift H coefficients + for ( i = 0 ; i < nm1 ; i++) { + j = nn - i - 1 ; + hr[j] = hr[j - 1] ; + hi[j] = hi[j - 1] ; + } + hr[0] = 0 ; + hi[0] = 0 ; + } + } +} + +//------------------------------------------------------------------------------ +// COMPUTES L2 FIXED-SHIFT H POLYNOMIALS AND TESTS FOR CONVERGENCE. +// INITIATES A VARIABLE-SHIFT ITERATION AND RETURNS WITH THE +// APPROXIMATE ZERO IF SUCCESSFUL. +// L2 - LIMIT OF FIXED SHIFT STEPS +// ZR,ZI - APPROXIMATE ZERO IF CONV IS .TRUE. +// CONV - LOGICAL INDICATING CONVERGENCE OF STAGE 3 ITERATION +//------------------------------------------------------------------------------ +void +Cpoly::fxshft( const int l2, double* zr, double* zi, bool* pbConv) +{ + bool bBol ; + bool bPasd ; + bool bTest ; + int i, j, n ; + double otr, oti, svsr, svsi ; + + + n = nn ; + polyev( nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi) ; + bTest = true ; + bPasd = false ; + + // Calculate first T = -P(S)/H(S) + calct( &bBol) ; + + // Main loop for second stage + for ( j = 1 ; j <= l2 ; j++) { + itercnt ++ ; + + otr = tr ; + oti = ti ; + + // Compute the next H Polynomial and new t + nexth( bBol) ; + calct( &bBol) ; + *zr = sr + tr ; + *zi = si + ti ; + + // Test for convergence unless stage 3 has failed once or this + // is the last H Polynomial + if ( ! ( bBol || ! bTest || j == 12)) + if ( cmod( tr - otr, ti - oti) < 0.5 * cmod( *zr, *zi)) { + if ( bPasd) { + // The weak convergence test has been passwed twice, start the third stage + // Iteration, after saving the current H polynomial and shift + for ( i = 0; i < n; i++ ) { + shr[i] = hr[i] ; + shi[i] = hi[i] ; + } + svsr = sr ; + svsi = si ; + vrshft( 10, zr, zi, pbConv) ; + if ( *pbConv) + return ; + + //The iteration failed to converge. Turn off testing and restore h,s,pv and T + bTest = false ; + for ( i = 0 ; i < n ; i++) { + hr[i] = shr[i] ; + hi[i] = shi[i] ; + } + sr = svsr ; + si = svsi ; + polyev( nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi) ; + calct( &bBol) ; + continue ; + } + bPasd = true ; + } + else + bPasd = false ; + } + + // Attempt an iteration with final H polynomial from second stage + vrshft( 10, zr, zi, pbConv) ; +} + +//------------------------------------------------------------------------------ +// CARRIES OUT THE THIRD STAGE ITERATION. +// L3 - LIMIT OF STEPS IN STAGE 3. +// ZR,ZI - ON ENTRY CONTAINS THE INITIAL ITERATE, IF THE +// ITERATION CONVERGES IT CONTAINS THE FINAL ITERATE ON EXIT. +// CONV - .TRUE. IF ITERATION CONVERGES +//------------------------------------------------------------------------------ +void +Cpoly::vrshft( const int l3, double* zr, double* zi, bool* pbConv) +{ + bool bBol ; + bool bFlag ; + bool bNext ; + int i, j ; + double mp, ms, omp, relstp, r1, r2, tp ; + + + *pbConv = false ; + bFlag = false ; + sr = *zr ; + si = *zi ; + + // Main loop for stage three + for ( i = 1 ; i <= l3 ; i++) { + itercnt ++ ; + // Evaluate P at S and test for convergence + polyev( nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi) ; + mp = cmod( pvr, pvi) ; + ms = cmod( sr, si) ; + if ( mp <= 20 * errev( nn, qpr, qpi, ms, mp, are, mre)) { + // Polynomial value is smaller in value than a bound on the error + // in evaluationg P, terminate the iteration + *pbConv = true ; + *zr = sr ; + *zi = si ; + return ; + } + bNext = false ; + if ( i != 1) { + if ( ! ( bFlag || mp < omp || relstp >= 0.05)) { + // Iteration has stalled. Probably a cluster of zeros. Do 5 fixed + // shift steps into the cluster to force one zero to dominate + tp = relstp ; + bFlag = true ; + if ( relstp < eta) + tp = eta ; + r1 = sqrt( tp) ; + r2 = sr * ( 1 + r1 ) - si * r1 ; + si = sr * r1 + si * ( 1 + r1) ; + sr = r2 ; + polyev( nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi) ; + for ( j = 1 ; j <= 5 ; j++) { + calct( &bBol) ; + nexth( bBol) ; + } + omp = infin ; + bNext = true ; + } + + // Exit if polynomial value increase significantly + if ( ! bNext && mp * 0.1 > omp) + return ; + } + + // eventuale salvataggio dato + if ( ! bNext) + omp = mp ; + + // Calculate next iterate + calct( &bBol) ; + nexth( bBol) ; + calct( &bBol) ; + if ( ! bBol) { + relstp = cmod( tr, ti) / cmod( sr, si) ; + sr += tr ; + si += ti ; + } + } +} + +//------------------------------------------------------------------------------ +// COMPUTES T = -P(S)/H(S). +// bool - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO. +//------------------------------------------------------------------------------ +void +Cpoly::calct( bool* pbBol) +{ + int n ; + double hvr, hvi ; + + + n = nn ; + + // evaluate h(s) + polyev( n - 1, sr, si, hr, hi, qhr, qhi, &hvr, &hvi) ; + *pbBol = ( cmod( hvr, hvi ) <= are * 10 * cmod( hr[n - 1], hi[n - 1])) ; + if ( ! *pbBol) { + cdivid( -pvr, -pvi, hvr, hvi, &tr, &ti) ; + return ; + } + + tr = 0 ; + ti = 0 ; +} + +//------------------------------------------------------------------------------ +// CALCULATES THE NEXT SHIFTED H POLYNOMIAL. +// bool - LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO +//------------------------------------------------------------------------------ +void +Cpoly::nexth( bool bBol) +{ + int j, n ; + double t1, t2 ; + + + n = nn ; + if ( ! bBol) { + for ( j = 1 ; j < n ; j++) { + t1 = qhr[j - 1] ; + t2 = qhi[j - 1] ; + hr[j] = tr * t1 - ti * t2 + qpr[j] ; + hi[j] = tr * t2 + ti * t1 + qpi[j] ; + } + hr[0] = qpr[0] ; + hi[0] = qpi[0] ; + return ; + } + + // If h[s] is zero replace H with qh + for ( j = 1 ; j < n ; j++) { + hr[j] = qhr[j - 1] ; + hi[j] = qhi[j - 1] ; + } + hr[0] = 0 ; + hi[0] = 0 ; +} + +//------------------------------------------------------------------------------ +// EVALUATES A POLYNOMIAL P AT S BY THE HORNER RECURRENCE +// PLACING THE PARTIAL SUMS IN Q AND THE COMPUTED VALUE IN PV. +//------------------------------------------------------------------------------ +void +Cpoly::polyev( const int nn, const double sr, const double si, const double pr[], const double pi[], + double qr[], double qi[], double *pvr, double *pvi ) +{ + int i ; + double t ; + + + qr[0] = pr[0] ; + qi[0] = pi[0] ; + *pvr = qr[0] ; + *pvi = qi[0] ; + + for ( i = 1 ; i <= nn ; i++) { + t = ( *pvr) * sr - ( *pvi) * si + pr[i] ; + *pvi = ( *pvr) * si + ( *pvi) * sr + pi[i] ; + *pvr = t ; + qr[i] = *pvr ; + qi[i] = *pvi ; + } +} + +//------------------------------------------------------------------------------ +// BOUNDS THE ERROR IN EVALUATING THE POLYNOMIAL BY THE HORNER RECURRENCE. +// QR,QI - THE PARTIAL SUMS +// MS -MODULUS OF THE POINT +// MP -MODULUS OF POLYNOMIAL VALUE +// ARE, MRE -ERROR BOUNDS ON COMPLEX ADDITION AND MULTIPLICATION +//------------------------------------------------------------------------------ +double +Cpoly::errev( const int nn, const double qr[], const double qi[], const double ms, const double mp, + const double are, const double mre ) +{ + int i ; + double e ; + + + e = cmod( qr[0], qi[0]) * mre / ( are + mre) ; + for ( i = 0 ; i <= nn ; i++) + e = e * ms + cmod( qr[i], qi[i]) ; + + return ( e * ( are + mre ) - mp * mre) ; +} + +//------------------------------------------------------------------------------ +// CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A +// POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS. +//------------------------------------------------------------------------------ +void +Cpoly::cauchy( const int nn, double pt[], double q[], double* fn_val) +{ + int i, n ; + double x, xm, f, dx, df ; + + + pt[nn] = -pt[nn] ; + + // Compute upper estimate bound + n = nn ; + x = exp( ( log( - pt[n]) - log( pt[0])) / (double) n) ; + if ( pt[n - 1] != 0) { + // Newton step at the origin is better, use it + xm = -pt[nn] / pt[n - 1] ; + if ( xm < x) + x = xm ; + } + + // Chop the interval (0,x) until f < 0 + while ( true) { + xm = x * 0.1 ; + f = pt[0] ; + for ( i = 1 ; i <= nn ; i++) + f = f * xm + pt[i] ; + if ( f <= 0) + break ; + x = xm ; + } + dx = x ; + + // Do Newton iteration until x converges to two decimal places + while ( fabs( dx / x ) > 0.005) { + q[0] = pt[0] ; + for ( i = 1 ; i <= nn ; i++) + q[i] = q[i - 1] * x + pt[i] ; + f = q[nn] ; + df = q[0] ; + for ( i = 1 ; i < n ; i++) + df = df * x + q[i] ; + dx = f / df ; + x -= dx ; + itercnt ++ ; + } + + *fn_val = x ; +} + +//------------------------------------------------------------------------------ +// RETURNS A SCALE FACTOR TO MULTIPLY THE COEFFICIENTS OF THE POLYNOMIAL. +// THE SCALING IS DONE TO AVOID OVERFLOW AND TO AVOID UNDETECTED UNDERFLOW +// INTERFERING WITH THE CONVERGENCE CRITERION. THE FACTOR IS A POWER OF THE BASE. +// PT - MODULUS OF COEFFICIENTS OF P +// ETA, INFIN, SMALNO, BASE - CONSTANTS DESCRIBING THE FLOATING POINT ARITHMETIC. +//------------------------------------------------------------------------------ +double +Cpoly::scale( const int nn, const double pt[], const double eta, + const double infin, const double smalno, const double base) +{ + int i, l ; + double hi, lo, max, min, x, sc ; + double fn_val ; + + + // Find largest and smallest moduli of coefficients + hi = sqrt( infin) ; + lo = smalno / eta ; + max = 0 ; + min = infin ; + + for ( i = 0 ; i <= nn ; i++) { + x = pt[i] ; + if ( x > max) + max = x ; + if ( x != 0 && x < min) + min = x ; + } + + // Scale only if there are very large or very small components + fn_val = 1 ; + if ( min >= lo && max <= hi) + return fn_val ; + x = lo / min ; + if ( x <= 1) + sc = 1 / ( sqrt( max)* sqrt( min)) ; + else { + sc = x; + if ( infin / sc > max) + sc = 1 ; + } + l = (int)( log( sc) / log( base) + 0.5) ; + fn_val = pow( base, l) ; + return fn_val ; +} + +//------------------------------------------------------------------------------ +// COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW. +//------------------------------------------------------------------------------ +void +Cpoly::cdivid( const double ar, const double ai, const double br, const double bi, double* cr, double* ci) +{ + double r, dinv, t, infin ; + + + if ( br == 0 && bi == 0) { + // Division by zero, c = infinity + mcon( &t, &infin, &t, &t) ; + *cr = infin ; + *ci = infin ; + return ; + } + + if ( fabs( br) < fabs( bi)) { + r = br / bi ; + dinv = 1.0 / ( bi + r * br) ; + *cr = ( ar * r + ai) * dinv ; + *ci = ( ai * r - ar) * dinv ; + return ; + } + + r = bi / br ; + dinv = 1.0 / ( br + r * bi) ; + *cr = ( ar + ai * r) * dinv ; + *ci = ( ai - ar * r) * dinv ; +} + +//------------------------------------------------------------------------------ +// MODULUS OF A COMPLEX NUMBER AVOIDING OVERFLOW. +//------------------------------------------------------------------------------ +double +Cpoly::cmod( const double r, const double i) +{ + double ar, ai ; + + + ar = fabs( r) ; + ai = fabs( i) ; + if ( ar < ai) + return ( ai * sqrt( 1.0 + ( ar * ar) / ( ai * ai))) ; + + if ( ar > ai) + return ( ar * sqrt( 1.0 + ( ai * ai) / ( ar * ar))) ; + + return ( ar * sqrt( 2.0)) ; +} + +//------------------------------------------------------------------------------ +// MCON PROVIDES MACHINE CONSTANTS USED IN VARIOUS PARTS OF THE PROGRAM. +// THE USER MAY EITHER SET THEM DIRECTLY OR USE THE STATEMENTS BELOW TO +// COMPUTE THEM. THE MEANING OF THE FOUR CONSTANTS ARE - +// ETA THE MAXIMUM RELATIVE REPRESENTATION ERROR WHICH CAN BE DESCRIBED +// AS THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT +// 1.0_dp + ETA > 1.0. +// INFINY THE LARGEST FLOATING-POINT NUMBER +// SMALNO THE SMALLEST POSITIVE FLOATING-POINT NUMBER +// BASE THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED +//------------------------------------------------------------------------------ +void +Cpoly::mcon( double* eta, double* infiny, double* smalno, double* base) +{ + *base = _DBL_RADIX ; + *eta = DBL_EPSILON ; + *infiny = DBL_MAX ; + *smalno = DBL_MIN ; +} diff --git a/JenkinsTraub.h b/JenkinsTraub.h new file mode 100644 index 0000000..a7b567f --- /dev/null +++ b/JenkinsTraub.h @@ -0,0 +1,90 @@ +//---------------------------------------------------------------------------- +// EgalTech 2013-2014 +//---------------------------------------------------------------------------- +// File : JenkinsTraub.h Data : 08.01.14 Versione : 1.5a1 +// Contenuto : Dichiarazione classi per il calcolo degli zeri di polinomi. +// +// +// +// Modifiche : 08.01.14 DS Creazione modulo. +// +// +//---------------------------------------------------------------------------- + +#pragma once + + +//---------------------------------------------------------------------------- +const int POLY_MAXDEG = 32 ; + + +//--------------------------- Class Rpoly ------------------------------------ +class Rpoly { + public : // methods + int Calculate( const double* op, int degree, double* zeror, double* zeroi) ; + + private : // methods + void quad( double a, double b1, double c, double* sr, double* si, + double* lr, double* li) ; + void fxshfr( int l2, int* nz) ; + void quadit( double* uu, double* vv, int* nz) ; + void realit( double sss, int* nz, bool* pIflag) ; + void calcsc( int* type) ; + void nextk( int* type) ; + void newest( int type, double* uu,double* vv) ; + void quadsd( int n, double* u, double* v, double* p, double* q, + double* a, double* b) ; + + public : // members + int itercnt ; + + private : // members + int n, nn, nmi ; + double sr, si, u, v, a, b, c, d, a1, a2 ; + double a3, a6, a7, e, f, g, h, szr, szi, lzr, lzi ; + double eta, are, mre ; + double p[POLY_MAXDEG+1] ; + double qp[POLY_MAXDEG+1] ; + double k[POLY_MAXDEG+1] ; + double qk[POLY_MAXDEG+1] ; +} ; + +//--------------------------- Class Cpoly ------------------------------------ +class Cpoly { + public : // methods + int Calculate( const double* opr, const double* opi, int degree, double* zeror, double* zeroi) ; + + private : // methods + void noshft( const int l1) ; + void fxshft( const int l2, double* zr, double* zi, bool* pbConv) ; + void vrshft( const int l3, double* zr, double* zi, bool* pbConv) ; + void calct( bool* pbBol) ; + void nexth( bool bBol) ; + void polyev( const int nn, const double sr, const double si, const double pr[], const double pi[], + double qr[], double qi[], double *pvr, double *pvi) ; + double errev( const int nn, const double qr[], const double qi[], + const double ms, const double mp, const double are, const double mre) ; + void cauchy( const int nn, double pt[], double q[], double *fn_val) ; + double scale( const int nn, const double pt[], const double eta, const double infin, + const double smalno, const double base) ; + void cdivid( const double ar, const double ai, const double br, const double bi, double *cr, double *ci) ; + double cmod( const double r, const double i) ; + void mcon( double *eta, double *infiny, double *smalno, double *base) ; + + public : // members + int itercnt ; + + private : // members + int nn ; + double sr, si, tr, ti, pvr, pvi, are, mre, eta, infin ; + double pr[POLY_MAXDEG+1] ; + double pi[POLY_MAXDEG+1] ; + double hr[POLY_MAXDEG+1] ; + double hi[POLY_MAXDEG+1] ; + double qpr[POLY_MAXDEG+1] ; + double qpi[POLY_MAXDEG+1] ; + double qhr[POLY_MAXDEG+1] ; + double qhi[POLY_MAXDEG+1] ; + double shr[POLY_MAXDEG+1] ; + double shi[POLY_MAXDEG+1] ; +} ; diff --git a/PolynomialZeros.cpp b/PolynomialZeros.cpp new file mode 100644 index 0000000..1dc4541 --- /dev/null +++ b/PolynomialZeros.cpp @@ -0,0 +1,236 @@ +//---------------------------------------------------------------------------- +// EgalTech 2013-2013 +//---------------------------------------------------------------------------- +// File : PolynomialZeros.cpp Data : 08.01.14 Versione : 1.5a1 +// Contenuto : Funzione per il calcolo degli zeri di polinomi. +// +// +// +// Modifiche : 08.01.14 DS Creazione modulo. +// +// +//---------------------------------------------------------------------------- + +//--------------------------- Include ---------------------------------------- +#include "stdafx.h" +#include "JenkinsTraub.h" +#include "\EgtDev\Include\ENkPolynomialZeros.h" +#include + + +//--------------------------------- Prototipi locali -------------------------------- +static void SortRoots( int nNum, double adRoot[]) ; +static void SortRoots( int nNum, Complex acRoot[]) ; + + +//---------------------------------------------------------------------------- +int +PolynomialZeros( int nDegree, double adPoly[], double adRoot[], int* pnIter) +{ + int i ; + int j ; + 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 ; + + // se il coefficiente del grado più alto è zero, diminuisco il grado + while ( nDegree >= 0 && fabs( adPoly[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] = adPoly[nDegree-i] ; + + // calcolo gli zeri + nZeros = cRpoly.Calculate( dPreal, nDegree, dZreal, dZcplx) ; + + // assegno gli zeri reali ai parametri di ritorno + for ( i = 0, j = 0 ; i < nZeros ; i++) { + if ( fabs( dZcplx[i]) < 100 * DBL_EPSILON) { + adRoot[j] = dZreal[i] ; + j ++ ; + } + } + nZeros = j ; + + // ordino le radici in senso decrescente + SortRoots( nZeros, adRoot) ; + + // assegno il numero di iterazioni + if ( pnIter != NULL) + *pnIter = cRpoly.itercnt ; + + return nZeros ; +} + +//---------------------------------------------------------------------------- +int +PolynomialZeros( int nDegree, Complex acPoly[], Complex acRoot[], 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 ; + + // se il coefficiente del grado più alto è zero, diminuisco il grado + while ( nDegree >= 0 && m2( acPoly[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] = acPoly[nDegree-i].re ; + + // ricavo i coefficienti complessi ( e verifico se non nulli) + bCplx = false ; + for ( i = 0 ; i <= nDegree ; i++) { + dPcplx[i] = acPoly[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 + for ( i = 0 ; i < nZeros ; i++) { + acRoot[i].re = dZreal[i] ; + acRoot[i].im = dZcplx[i] ; + } + + // annullo le parti reali e immaginarie molto piccole + for ( i = 0 ; i < nZeros ; i++) { + if ( fabs( acRoot[i].re) < 100 * DBL_EPSILON) + acRoot[i].re = 0 ; + if ( fabs( acRoot[i].im) < 100 * DBL_EPSILON) + acRoot[i].im = 0 ; + } + + // ordino le radici in senso decrescente della parte reale + SortRoots( nZeros, acRoot) ; + + // 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) ; +} diff --git a/resource.h b/resource.h new file mode 100644 index 0000000..9614989 Binary files /dev/null and b/resource.h differ diff --git a/stdafx.cpp b/stdafx.cpp new file mode 100644 index 0000000..6826075 --- /dev/null +++ b/stdafx.cpp @@ -0,0 +1,7 @@ +// stdafx.cpp : file di origine che include solo le inclusioni standard +// EgtGeometry.pch sarà l'intestazione precompilata +// stdafx.obj conterrà le informazioni sui tipi precompilati + +#include "stdafx.h" + + diff --git a/stdafx.h b/stdafx.h new file mode 100644 index 0000000..06eabfd --- /dev/null +++ b/stdafx.h @@ -0,0 +1,31 @@ +// stdafx.h : file di inclusione per file di inclusione di sistema standard +// o file di inclusione specifici del progetto utilizzati di frequente, ma +// modificati raramente +// + +#pragma once + +#include "/EgtDev/Include/EgtTargetVer.h" + +#include +#include +#include +#include + +// in Debug riconoscimento memory leakage +#if defined( _DEBUG) + #define _CRTDBG_MAP_ALLOC + #include + #include +#endif + +// in Debug controllo iteratori +#if defined( _DEBUG) + #define _SECURE_SCL 1 +#else + #define _SECURE_SCL 0 +#endif + +#include "/EgtDev/Include/EgtLibVer.h" + +#pragma comment(lib, EGTLIBDIR "EgtGeneral" EGTLIBVER ".lib")