#include 
#include 
/*  compile  with a command like
cc  wgs.c -lm
or
cc -O3 wgs.s -lm
*/

main()
/* Sample driver program to call wgs */
{
int i,n,m;
double ar[11],ai[11],rr[11],ri[11];
/* solve  x^10-x^2+1=0 for x */
n = 10;
for (i=0; i<=n; i++) {
 ar[i] = 0.;
 ai[i] = 0.; 
 rr[i] = 0.;
 ri[i] = 0.; 
}
ar[10] = 1.;
ar[2] = -1.;
ar[0] = 1.;
printf(" The polynomial has degree %d\n",n);
printf(" The coefficients are:\n");
printf("  degree        real part                   imaginary part\n");
for (i=0; i<=n; i++) {
 printf("%8d  %25.16e  +  %25.16e * i\n",i,ar[i],ai[i]);
}
m = wgs(n,ar,ai,rr,ri);
printf("\n");
printf("\n");
printf(" The roots of the polynomial are:\n");
printf("  root #        real part                   imaginary part\n");
for (i=1; i<=n; i++) {
 printf("%8d  %25.16e  +  %25.16e * i\n",i,rr[i],ri[i]);
}
exit(0);
}

/* wgs
   Input:
    nn  an integer which is the degree of the polynomial
    ar  an array of doubles (the real part of the coefficients)
          indexed from 0 to nn  (of size nn+1);
    ai  an array of doubles (the imaginary part of the coefficients)
          indexed from 0 to nn  (of size nn+1);
   Output:
         the function value is an int (the number of non-zero roots)
    rr  an array of doubles (the real part of the roots)
          indexed from 0 to nn (of size nn+1) rr[0] is not used
    ri  an array of doubles (the imaginary part of the roots)
          indexed from 0 to nn (of size nn+1) ri[0] is not used
*/
int wgs(nn,ar,ai,rr,ri)
int nn;
double ar[],ai[],rr[],ri[];
{
int i,k,n,Mit,it;
double t,tr,ti,pr,pi,dr,di,xkr,xki,eps,eps1,err,err1;
for (i=0; ieps && err1>eps1; it++) {
  err = 0;
  err1 = 0;
  for (k=1; k<=n; k++) {
    xkr = rr[k]; xki = ri[k];
    pr = xkr+ar[n-1]; pi = xki+ai[n-1];
    for (i=n-2; i>=0; i--) {
      t = pr*xkr-pi*xki+ar[i];
      pi = pr*xki+pi*xkr+ai[i];
      pr = t;
    } 
    dr = 1.; di = 0.;
    for (i=1; ierr) err = t;
    rr[k] = xkr-tr; ri[k] = xki-ti;
    ti = fabs(rr[k])+fabs(ri[k]);
    tr = ti+fabs(ar[n-1])+fabs(ai[n-1]);
    for (i=n-2; i>=0; i--) {
      tr = tr*ti+fabs(ar[i])+fabs(ai[i]);
    } 
    ti = (fabs(pr)+fabs(pi))/tr;
    if (ti>err1) err1 = ti;
  }
}
return n;
}