proj/globalhs/libs/pkbgn/zroots.c

// $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/libs/pkbgn/zroots.c,v 1.1 2013/04/28 07:50:36 tplarson Exp $

#include <math.h>
#include "complex.h"

#define EPS 2.0e-6
#define MAXM 100

void zroots(a,m,roots,polish)
fcomplex a[],roots[];
int m,polish;
{
    int jj,j,i;
    fcomplex x,b,c,ad[MAXM];
    void laguer();

    for (j=0;j<=m;j++) ad[j]=a[j];
    for (j=m;j>=1;j--) {
        x=Complex(0.0,0.0);
        laguer(ad,j,&x,EPS,0);
        if (fabs(x.i) <= (2.0*EPS*fabs(x.r))) x.i=0.0;
        roots[j]=x;
        b=ad[j];
        for (jj=j-1;jj>=0;jj--) {
            c=ad[jj];
            ad[jj]=b;
            b=Cadd(Cmul(x,b),c);
        }
    }
    if (polish)
        for (j=1;j<=m;j++)
            laguer(a,m,&roots[j],EPS,1);
    for (j=2;j<=m;j++) {
        x=roots[j];
        for (i=j-1;i>=1;i--) {
            if (roots[i].r <= x.r) break;
            roots[i+1]=roots[i];
        }
        roots[i+1]=x;
    }
}

#undef EPS
#undef MAXM

---