From 44b5a1007b36469e619a4ff4953f9a26640e7120 Mon Sep 17 00:00:00 2001 From: Dario Sassi Date: Wed, 8 Jan 2014 14:40:04 +0000 Subject: [PATCH] EgtNumKernel 1.5a1 : Calcolo degli zeri di polinomi reali e complessi. --- Complex.cpp | 271 +++++++ DllMain.h | 20 + ENkDllMain.cpp | 76 ++ EgtNumKernel.rc | Bin 0 -> 7626 bytes EgtNumKernel.sln | 20 + EgtNumKernel.vcxproj | 125 +++ EgtNumKernel.vcxproj.filters | 68 ++ JenkinsTraub.cpp | 1390 ++++++++++++++++++++++++++++++++++ JenkinsTraub.h | 90 +++ PolynomialZeros.cpp | 236 ++++++ resource.h | Bin 0 -> 826 bytes stdafx.cpp | 7 + stdafx.h | 31 + 13 files changed, 2334 insertions(+) create mode 100644 Complex.cpp create mode 100644 DllMain.h create mode 100644 ENkDllMain.cpp create mode 100644 EgtNumKernel.rc create mode 100644 EgtNumKernel.sln create mode 100644 EgtNumKernel.vcxproj create mode 100644 EgtNumKernel.vcxproj.filters create mode 100644 JenkinsTraub.cpp create mode 100644 JenkinsTraub.h create mode 100644 PolynomialZeros.cpp create mode 100644 resource.h create mode 100644 stdafx.cpp create mode 100644 stdafx.h 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 0000000000000000000000000000000000000000..a8d3075e67da1772ba8d3c109eb414985516d182 GIT binary patch literal 7626 zcmeI1+ioI95QghIQr=;LbKzBzwefl{a=A8`Sc}0j2PKjv zrpFxYSzk;fgqohIsjja6>(tZp?>~FCXNjHJ&~9vEW1G_YHiOoMHnm5)vJrR3_7lIv z26kunvMsde74l+_#Rm^lbcl`(h=> zt&uy9ukXf5_9r{HmbI;86+5-2_3Xl`P%B(FtZsimY4Lw%M^Nfqhxw&L@qR(O&?i;{ zmn_D8P#W8xtbYbN!lBKof9L;{RUflX;ZcQtK+9TWvkk7W`q%lm82{oS7BuqS@d#G( zq(1A;=Nkd}Ah%e3&@RE!_r%|bwN9M}?Bd%q7l%XE-r?zx`=PyC;sMiCJPjRIKSXhL zYrCL#f=BsyY7sjfBI_8R>fu>+evKSJM1K z@P^qwyK7~LQLfqIOcZ7*qx(nJReqIBSyB7I3B54MV&!$ut8pF49&1_T0b2-z0hrIw z(p%&ClZ0G!>@wp4UQs40e_>aot3%McMoRLAHf+vyl6MP%U{l>?y=DeJvSk{N3-mhx z`&-7o&YtJYJZ9V@ULncP=Zn4ct^6pScbzvaj|J@wD)T((FX+R9O>JdJ|gRAKX6T6A8?(`ujyS@*j0o+&T^p=WhMB! zxvls4_H5r;?+d$xS#{C|<986k74L3mQp7UCS0a zGaIqz(++S9WjCS~;q>g3cgs&NI6d=xIz7{Aneuw}^epBGc8X6YHk~mA>;jg>>VtNn zUv)Abo~u_YgVw4ubUj~w60!Vi#pToK!y>auhR)*i%pv<^zmk;j*X$X6$u4=5xr8@< z(-J+Lhr}9v$f?n5=Sja0#kue$-X`o;OUNI=#q!*=Hd=AI`3Th;<*F|m@j0-@rsWmu kI%pUAZ~ew4FJAwRD_~hH9MalI!__{n%A~KbFO8P}0HRE#PXGV_ literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..961498951d9c0a14b27a87d097b985f99d6d1cc7 GIT binary patch literal 826 zcmb7?%}c{T5XIkF@PAnFC|Ior&oSwtv`M8gMNh$|jZ|s_`S1hrUsr#-vCw)5m0`0p zvv1zKoz3@GPd$BLktRwNXrLQb_Y8=1zAG(qmS! zZj~^tp_|{3t~0z5iUU{MR%|+IreHnL6Kco0=M`eUp{~h{@!P)b9qjDv%+lZ9DK!J# zp&^y!iK$8OeEbw89t;V-`D1;ZJ3Exyhd0DL`URcLMwjaA3cEO1xn+;I!_?ilOSggh af6Cig+SlmjwG-Oswtw!zy>~mM{rCqDb9R*g literal 0 HcmV?d00001 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")