/* Translated from F77 to C, rjrw 10/04/2000 */
/* replaced 'bool' by 'boolvar' to get it to compile on my
linux machine, DJB Aug 02 2000 */
/* algorithm 419 collected algorithms from acm.
algorithm appeared in comm. acm, vol. 15, no. 02, p. 097. */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/*
#if !defined(WIN32) && !defined(_WIN32) && !defined(__APPLE__) && !defined(__CYGWIN__)
#include <values.h>
#endif
*/
#include <float.h>
/* #define DEBUGMAIN */ /* Set up debugging main, etc. */
#include "cpoly.h"
/* Internal routines */
static void noshft(int l1);
static int fxshft(int l2, double *zr, double *zi);
static int vrshft(int l3, double *zr, double *zi);
static int calct(void);
static void nexth(int boolvar);
static void polyev(int nn, double sr, double si, double pr[], double pi[],
double qr[], double qi[], double *pvr, double *pvi);
static double errev(int nn, double qr[], double qi[], double ms, double mp);
static double cauchy(int nn, double pt[], double q[]);
static double scale(int nn, double pt[]);
static void cdivid(double ar, double ai, double br, double bi,
double *cr, double *ci);
static double cmod(double r, double i);
static void mcon(void);
static int init(int nncr);
/* Internal global variables */
static double *pr,*pi,*hr,*hi,*qpr,*qpi,*qhr,*qhi,*shr,*shi;
static double sr,si,tr,ti,pvr,pvi,are,mre,eta,infin,smalno,base;
static int nn;
#ifdef DEBUGMAIN
/* driver to test cpoly */
int main()
{
int fail;
double p[50],pi[50],zr[50],zi[50];
int i;
printf("Example 1. polynomial with zeros 1,2,...,10.\n");
p[0]=1L;
p[1]=-55L;
p[2]=1320L;
p[3]=-18150L;
p[4]=157773L;
p[5]=-902055L;
p[6] = 3416930L;
p[7]=-8409500L;
p[8]=12753576L;
p[9]=-10628640L;
p[10]=3628800L;
for (i=0;i<11;i++)
pi[i]=0;
prtc(11,p,pi);
fail = cpoly(p,pi,10,zr,zi);
if(fail)
printf("cpoly has failed on this example\n");
prtz (10,zr,zi);
printf("Example 2. zeros on imaginary axis degree 3.\n");
p[0]=1;
p[1]=0;
p[2]=-10001.0001L;
p[3]=0;
pi[0]=0;
pi[1]=-10001.0001L;
pi[2]=0;
pi[3]=1;
prtc(4,p,pi);
fail = cpoly(p,pi,3,zr,zi);
if (fail)
printf("cpoly has failed on this example\n");
prtz(3,zr,zi);
printf("Example 3. zeros at 1+i,1/2*(1+i)....1/(2**-9)*(1+i)\n");
p[0]=1.0;
p[1]=-1.998046875L;
p[2]=0.0;
p[3]=.7567065954208374L;
p[4]=-.2002119533717632L;
p[5]=1.271507365163416e-2L;
p[6]=0;
p[7]=-1.154642632172909e-5L;
p[8]=1.584803612786345e-7L;
p[9]=-4.652065399568528e-10L;
p[10]=0;
pi[0]=0;
pi[1]=p[1];
pi[2]=2.658859252929688L;
pi[3]=-7.567065954208374e-1L;
pi[4]=0;
pi[5]=p[5];
pi[6]=-7.820779428584501e-4L;
pi[7]=-p[7];
pi[8]=0;
pi[9]=p[9];
pi[10]=9.094947017729282e-13L;
prtc(11,p,pi);
fail = cpoly(p,pi,10,zr,zi);
if (fail)
printf("cpoly has failed on this example\n");
prtz(10,zr,zi);
printf("Example 4. multiple zeros\n");
p[0]=1L;
p[1]=-10L;
p[2]=3L;
p[3]=284L;
p[4]=-1293L;
p[5]=2374L;
p[6]=-1587L;
p[7]=-920L;
p[8]=2204L;
p[9]=-1344L;
p[10]=288L;
pi[0]=0;
pi[1]=-10L;
pi[2]=100L;
pi[3]=-334L;
pi[4]=200L;
pi[5]=1394L;
pi[6] =-3836L;
pi[7]=4334L;
pi[8]=-2352L;
pi[9]=504L;
pi[10]=0;
prtc(11,p,pi);
fail = cpoly(p,pi,10,zr,zi);
if (fail)
printf("cpoly has failed on this example\n");
prtz(10,zr,zi);
printf("Example 5. 12 zeros evenly distributed on a circle of radius 1. centered at 0+2i.\n");
p[0]=1L;
p[1]=0;
p[2]=-264L;
p[3]=0;
p[4]=7920L;
p[5]=0;
p[6]=-59136L;
p[7]=0;
p[8]=126720L;
p[9]=0;
p[10]=-67584L;
p[11]=0;
p[12]=4095L;
pi[0]=0;
pi[1]=-24L;
pi[2]=0;
pi[3]=1760L;
pi[4]=0;
pi[5]=-25344L;
pi[6]=0;
pi[7]=101376L;
pi[8]=0;
pi[9]=-112640L;
pi[10]=0;
pi[11]=24576L;
pi[12]=0;
prtc(13,p,pi);
fail = cpoly(p,pi,12,zr,zi);
if(fail)
printf("cpoly has failed on this example\n");
prtz(12,zr,zi);
return 0;
}
void prtc(int n, double p[], double q[])
{
int i;
printf("Coefficients\n");
for (i=0;i<n;i++)
printf("%26.16g %26.16g\n",p[i],q[i]);
}
void prtz(int n,double zr[], double zi[])
{
int i;
printf("Zeroes\n");
for (i=0;i<n;i++)
printf("%26.16g %26.16g\n",zr[i],zi[i]);
}
#endif
/* rjrw 10/04/2000: fix for cos 94: was -.060756474L */
#define COSR (-.069756474L)
#define SINR (.99756405L)
int cpoly(double opr[], double opi[], int degree,
double zeror[], double zeroi[])
{
/* Finds the zeros of a complex polynomial.
opr, opi - double precision vectors of real and imaginary parts
of the coefficients in order of decreasing powers.
degree - integer degree of polynomial
zeror, zeroi
- output double precision vectors of real and imaginary
parts of the zeros.
fail - output logical parameter, TRUE if leading coefficient
is zero, if cpoly has found fewer than degree zeros,
or if there is another internal error.
The program has been written to reduce the chance of overflow
occurring. If it does occur, there is still a possibility that
the zerofinder will work provided the overflowed quantity is
replaced by a large number. */
double xx,yy,xxx,zr,zi,bnd;
int fail,conv;
int cnt1,cnt2,i,idnn2;
/* initialization of constants */
nn = degree+1;
if (!init(nn)) {
fail = TRUE;
return fail;
}
xx = .70710678L;
yy = -xx;
fail = FALSE;
/* algorithm fails if the leading coefficient is zero. */
if (opr[0] == 0.0 && opi[0] == 0.0) {
fail = TRUE;
return fail;
}
/* Remove the zeros at the origin if any */
while (opr[nn-1] == 0.0 && opi[nn-1] == 0.0) {
idnn2 = degree+1-nn;
zeror[idnn2] = 0.0;
zeroi[idnn2] = 0.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);
if (bnd != 1.0) {
for (i=0;i<nn;i++) {
pr[i] *= bnd;
pi[i] *= bnd;
}
}
while (!fail) {
/* Start the algorithm for one zero */
if (nn < 3) {
/* Calculate the final zero and return */
cdivid(-pr[1],-pi[1],pr[0],pi[0],&(zeror[degree-1]),&(zeroi[degree-1]));
return fail;
}
/* Calculate bnd, a lower bound on the modulus of the zeros */
for (i=0;i<nn;i++) {
shr[i] = cmod(pr[i],pi[i]);
}
bnd = cauchy(nn,shr,shi);
/* Outer loop to control 2 major passes with different sequences
of shifts */
fail = TRUE;
for(cnt1=1;fail && (cnt1<=2);cnt1++) {
/* First stage calculation, no shift */
noshft(5);
/* Inner loop to select a shift. */
for (cnt2=1;fail && (cnt2<10);cnt2++) {
/* Shift is chosen with 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;
/* Second stage calculation, fixed shift */
conv = fxshft(10*cnt2,&zr,&zi);
if (conv) {
/* The second stage jumps directly to the third stage iteration
If successful the zero is stored and the polynomial deflated */
idnn2 = degree+1-nn;
zeror[idnn2] = zr;
zeroi[idnn2] = zi;
nn--;
for(i=0;i<nn;i++) {
pr[i] = qpr[i];
pi[i] = qpi[i];
}
fail = FALSE;
}
/* 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 */
return fail;
}
static void noshft(int l1)
{
/* Computes the derivative polynomial as the initial h
polynomial and computes l1 no-shift h polynomials. */
double xni,t1,t2;
int i,j,jj,n = nn-1,nm1 = n-1,nm2=nm1-1;
for (i=0;i<n;i++) {
xni = n-i;
hr[i] = xni*pr[i]/((double)(n));
hi[i] = xni*pi[i]/((double)(n));
}
for (jj=0;jj<l1;jj++) {
if (cmod(hr[nm2],hi[nm2]) > eta*10.0*cmod(pr[nm2],pi[nm2])) {
cdivid(-pr[n],-pi[n],hr[nm1],hi[nm1],&tr,&ti);
for (i=0;i<nm1;i++) {
j = nm1-i;
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 = nm1-i;
hr[j] = hr[j-1];
hi[j] = hi[j-1];
}
hr[0] = 0.0;
hi[0] = 0.0;
}
}
}
static int fxshft(int l2, double *zr, double *zi)
/* 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 - Flag indicating convergence of stage 3 iteration
*/
{
double otr,oti,svsr,svsi;
int conv,test,pasd,boolvar;
int i,j,n = nn-1;
/* Evaluate p at s */
polyev(nn,sr,si,pr,pi,qpr,qpi,&pvr,&pvi);
test = TRUE;
pasd = FALSE;
/* Calculate first t = -p(s)/h(s) */
boolvar = calct();
/* Main loop for one second stage step */
for (j=0;j<l2;j++) {
otr = tr;
oti = ti;
/* Compute next h polynomial and new t */
nexth(boolvar);
boolvar = calct();
*zr = sr+tr;
*zi = si+ti;
/* Test for convergence unless stage 3 has failed once or
this is the last h polynomial */
if (!boolvar && test && j != l2) {
if (cmod(tr-otr,ti-oti) < .5*cmod(*zr,*zi)) {
if (pasd) {
/* The weak convergence test has been passed 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;
conv = vrshft(10,zr,zi);
if (conv)
return conv;
/* The iteration failed to converge
Turn off testing and restore h,s,pv and t */
test = 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);
boolvar = calct();
} else {
pasd = TRUE;
}
}
} else {
pasd = FALSE;
}
}
/* Attempt an iteration with final h polynomial from second stage */
conv = vrshft(10,zr,zi);
return conv;
}
static int vrshft(int l3, double *zr, double *zi)
/* Carries out the third stage iteration
l3 - Limit of steps in stage 3
zr,zi - On entry contains the initial iterate,
On exit, it contains the final iterate (if it converges).
conv - TRUE if iteration converges
*/
{
double mp,ms,omp,relstp,r1,r2,tp;
int i,j,conv,b,boolvar;
conv = FALSE;
b = FALSE;
sr = *zr;
si = *zi;
/* Main loop for stage three */
for (i=0; i<l3;i++) {
/* 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.0L*errev(nn,qpr,qpi,ms,mp)) {
/* Polynomial value is smaller in value than a bound on the error
in evaluating p, terminate the iteration */
conv = TRUE;
*zr = sr;
*zi = si;
return conv;
} else {
if (i!=0) {
if (!b && mp>=omp && relstp < .05L) {
/* Iteration has stalled, probably a cluster of zeros
Do 5 fixed shift steps into the cluster to force one zero
to dominate */
b = TRUE;
if (relstp < eta)
tp = eta;
else
tp = relstp;
r1 = sqrt(tp);
r2 = sr*(1.0L+r1)-si*r1;
si = sr*r1+si*(1.0L+r1);
sr = r2;
polyev(nn,sr,si,pr,pi,qpr,qpi,&pvr,&pvi);
for (j=0;j<5;j++) {
boolvar = calct();
nexth(boolvar);
}
omp = infin;
} else {
/* Exit if polynomial value increases significantly */
if (mp*0.1L > omp)
return conv;
omp = mp;
}
} else {
omp = mp;
}
}
/* Calculate next iterate. */
boolvar = calct();
nexth(boolvar);
boolvar = calct();
if (!boolvar) {
relstp = cmod(tr,ti)/cmod(sr,si);
sr += tr;
si += ti;
}
}
return conv;
}
static int calct(void)
/* Computes t = -p(s)/h(s)
Returns TRUE if h(s) is essentially zero
*/
{
double hvr,hvi;
int n = nn-1, boolvar;
/* Evaluate h(s) */
polyev(n,sr,si,hr,hi,qhr,qhi,&hvr,&hvi);
boolvar = (cmod(hvr,hvi) <= are*10.0*cmod(hr[n-1],hi[n-1]));
if (!boolvar) {
cdivid(-pvr,-pvi,hvr,hvi,&tr,&ti);
} else {
tr = 0.0;
ti = 0.0;
}
return boolvar;
}
static void nexth(int boolvar)
/* Calculates the next shifted h polynomial
boolvar - TRUE if h(s) is essentially zero
*/
{
double t1,t2;
int j,n = nn-1;
if (!boolvar) {
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];
} else {
/* 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.0;
hi[0] = 0.0;
}
}
static void polyev(int nn, double sr, double si, double pr[], double pi[],
double qr[], double qi[], double *tvr, double *tvi)
/* Evaluates a polynomial p at s by the Horner recurrence,
placing the partial sums in q and the computed value in pv
*/
{
double t, vr, vi;
int i;
qr[0] = pr[0];
qi[0] = pi[0];
vr = qr[0];
vi = qi[0];
for (i=1;i<nn;i++) {
t = vr*sr-vi*si+pr[i];
vi = vr*si+vi*sr+pi[i];
vr = t;
qr[i] = vr;
qi[i] = vi;
}
*tvr = vr;
*tvi = vi;
}
static double errev(int nn, double qr[], double qi[], double ms, double mp)
/* 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
*/
{
double e;
int i;
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;
}
static double cauchy(int nn, double pt[], double q[])
/* Cauchy computes a lower bound on the moduli of the zeros of a
polynomial - pt is the modulus of the coefficients
*/
{
double x,xm,f,dx,df;
int n=nn-1, nm=nn-2, i;
pt[n] = -pt[n];
/* Compute upper estimate of bound */
xm = exp( (log(-pt[n]) - log(pt[0]))/((double)n) );
if (pt[nm] != 0.0) {
/* If Newton step at the origin is better, use it */
x = -pt[n]/pt[nm];
if (x < xm)
xm = x;
}
/* Chop the interval (0,x) until f <= 0 */
do {
x = xm;
xm *= .1;
f = pt[0];
for (i=1;i<nn;i++)
f = f*xm+pt[i];
} while (f > 0.);
dx = x;
/* Do Newton iteration until x converges to two decimal places */
while (fabs(dx/x) > .005L) {
q[0] = pt[0];
for(i=1;i<nn;i++)
q[i] = q[i-1]*x+pt[i];
f = q[n];
df = q[0];
for (i=1;i<n;i++)
df = df*x+q[i];
dx = f/df;
x -= dx;
}
return x;
}
static double scale(int nn, double pt[])
/* 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
*/
{
double hi,lo,max,min,x,sc;
int i,l;
/* Find largest and smallest moduli of coefficients */
hi = sqrt(infin);
lo = smalno/eta;
max = 0.0;
min = infin;
for (i=0;i<nn;i++) {
x = pt[i];
if (x > max)
max = x;
if (x != 0.0 && x < min)
min = x;
}
/* Scale only if there are very large or very small components */
if (min >= lo && max <= hi)
return 1.0;
x = lo/min;
if (x <= 1.0L) {
sc = 1.0L/(sqrt(max)*sqrt(min));
} else {
sc = x;
if (infin/sc > max)
sc = 1.0;
}
l = log(sc)/log(base) + .500;
return pow(base,l);
}
static void cdivid(double ar, double ai, double br, double bi,
double *cr, double *ci)
/* Complex division c = a/b, avoiding overflow */
{
double r,d;
if (br == 0.0 && bi == 0.0) {
/* division by zero, c = infinity. */
*cr = infin;
*ci = infin;
} else if (fabs(br) < fabs(bi)) {
r = br/bi;
d = bi+r*br;
*cr = (ar*r+ai)/d;
*ci = (ai*r-ar)/d;
} else {
r = bi/br;
d = br+r*bi;
*cr = (ar+ai*r)/d;
*ci = (ai-ar*r)/d;
}
return;
}
static double cmod(double r, double i)
/* Modulus of a complex number avoiding overflow */
{
double ar,ai,f;
ar = fabs(r);
ai = fabs(i);
if (ar < ai) {
f = ar/ai;
return ai*sqrt(1.0+f*f);
} else if (ar > ai) {
f = ai/ar;
return ar*sqrt(1.0+f*f);
} else {
return ar*sqrt(2.0);
}
}
static void mcon()
/* 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.0d0 + eta is
greater than 1.0d0.
infin the largest floating-point number
smalno the smallest positive floating-point number
base the base of the floating-point number system used
Let t be the number of base-digits in each floating-point
number (double precision). Then eta is either .5*b**(1-t)
or b**(1-t) depending on whether rounding or truncation
is used.
Let m be the largest exponent and n the smallest exponent
in the number system. Then infiny is (1-base**(-t))*base**m
and smalno is base**n.
*/
{
/*
#if !defined(WIN32) && !defined(_WIN32) && !defined(__APPLE__) && !defined(__CYGWIN__)
base = 2;
eta = DBL_EPSILON;
smalno = MINDOUBLE;
infin = MAXDOUBLE;
#else
*/
base = 2;
eta = DBL_EPSILON;
smalno = DBL_MIN;
infin = DBL_MAX;
/* #endif */
#ifdef IBM360
/* These values for base,t,m,n correspond to the ibm/360. */
int m,n,t;
base = 16.0;
t = 14;
m = 63;
n = -65;
eta = pow(base,1-t);
infin = (base)*(1.0-pow(base,-t))*pow(base,m-1);
smalno = pow(base,n+3)/pow(base,3);
#endif
}
static int init(int nncr)
{
static int nmax=0;
if (nmax == 0) {
/* Set up once-off constants */
mcon();
/* are, mre - Error bounds on complex addition and multiplication,
cf e.g. errev() above */
are = eta;
mre = 2.0L*sqrt(2.0L)*eta;
} else if (nmax >= nncr) {
return TRUE; /* Present arrays are big enough */
} else {
/* Free old arrays (no need to preserve contents */
free(shi); free(shr); free(qhi); free(qhr);
free(qpi); free(qpr); free(hi); free(hr); free(pi); free(pr);
}
nmax = nncr;
pr = (double *) malloc(nmax*sizeof(double));
pi = (double *) malloc(nmax*sizeof(double));
hr = (double *) malloc(nmax*sizeof(double));
hi = (double *) malloc(nmax*sizeof(double));
qpr = (double *) malloc(nmax*sizeof(double));
qpi = (double *) malloc(nmax*sizeof(double));
qhr = (double *) malloc(nmax*sizeof(double));
qhi = (double *) malloc(nmax*sizeof(double));
shr = (double *) malloc(nmax*sizeof(double));
shi = (double *) malloc(nmax*sizeof(double));
if (!(pr && pi && hr && hi && qpr && qpi && qhr && qhi && shr && shi)) {
fprintf(stderr,"Couldn't allocate space for cpoly\n");
return FALSE;
} else {
return TRUE;
}
}