Egalware/EgtNumKernel
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
master
EgtNumKernel/JenkinsTraub.cpp
T
DarioS eeab41af4e EgtNumKernel 2.3g1 : 
- versione x64 compilata con Clang-cl/LLVM
- modifiche varie per eliminare warning più gravi di questo compilatore.
2021-07-20 18:13:21 +02:00

1400 lines
41 KiB
C++
Raw Permalink Blame History

//----------------------------------------------------------------------------
// 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"
//--------------------------- Constants --------------------------------------
#define RADIX_DOUBLE 2
#define EPS_DOUBLE 2.22e-16
#define SMALL_DOUBLE 2.3e-308
#define BIG_DOUBLE 1.7e+308
//--------------------------- 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 = RADIX_DOUBLE ;
eta = EPS_DOUBLE ;
infin = BIG_DOUBLE ;
smalno = SMALL_DOUBLE ;
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 = abs( 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] = abs( 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 ( abs( 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 = ( abs( k[n-1]) <= abs( 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 bFflag ;
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
bIflag = true ;
*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 ++) {
bFflag = true ;
// 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 = abs( ( vv - ovv) / vv) ;
if ( ss != 0.0)
ts = abs( ( 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) {
if ( bFflag && 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.
bIflag = true ;
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] ;
}
}
bFflag = false ;
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 ( abs( abs( szr) - abs( lzr)) > 0.01 * abs( lzr))
return ;
// Evaluate polynomial by quadratic synthetic division.
quadsd( n, &u, &v, p, qp, &a, &b) ;
mp = abs( a - szr * b) + abs( szi * b) ;
// Compute a rigorous bound on the rounding error in evaluating p.
zm = sqrt( abs( v)) ;
ee = 2.0 * abs( qp[0]) ;
t = -szr * b ;
for ( i = 1 ; i < n ; i ++) {
ee = ee * zm + abs( qp[i]) ;
}
ee = ee * zm + abs( a + t) ;
ee *= (5.0 * mre + 4.0 * are) ;
ee = ee - ( 5.0 * mre + 2.0 * are) * ( abs( a + t) + abs( b) * zm) ;
ee = ee + 2.0 * are * abs( t) ;
// Iteration has converged sufficiently if the polynomial value is less than 20 times this bound.
if ( mp <= 20.0 * ee) {
*nz = 2 ;
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 = abs( ( 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 = abs( pv) ;
// Compute a rigorous bound on the error in evaluating p.
ms = abs( s) ;
ee = ( mre / ( are + mre)) * abs( qp[0]) ;
for ( i = 1 ; i <= n ; i ++) {
ee = ee * ms + abs( qp[i]) ;
}
// Iteration has converged sufficiently if the polynomial value is less than 20 times this bound.
if ( mp <= 20.0 * (( are + mre) * ee - mre * mp)) {
*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 &&
! ( abs( t) > 0.001 * abs( 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 ( abs( kv) <= abs( 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 ( abs( kv) > ( abs( 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 ( abs( c) <= abs( k[n-1] * 100.0 * eta) &&
abs( d) <= abs( k[n-2] * 100.0 * eta)) {
*type = 3 ;
return ;
}
// Type=1 indicates that all formulas are divided by c.
if ( abs( d) < abs( 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 ( abs( a1) <= abs( temp) * eta * 10.0) {
// If a1 is nearly zero then use a special form of the recurrence.
k[0] = 0.0;
k[1] = -a7*qp[0] ;
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 ( abs( b) < abs( c)) {
if ( c < 0.0)
e = - a ;
else
e = a ;
e = b * ( b / abs( c)) - e ;
d = sqrt( abs( e)) * sqrt( abs( c)) ;
}
else {
e = 1.0 - ( a / b) *( c / b) ;
d = sqrt( abs( e)) * abs( b) ;
}
if ( e < 0.0) { /* complex conjugate zeros */
*sr = - b / a ;
*lr = *sr ;
*si = abs( 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 ( abs( 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 ( abs( br) < abs( 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 = abs( r) ;
double ai = abs( 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 &gt; 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 = RADIX_DOUBLE ;
*eta = EPS_DOUBLE ;
*infiny = BIG_DOUBLE ;
*smalno = SMALL_DOUBLE ;
}
Reference in New Issue View Git Blame Copy Permalink