## Computing Roots of Polynomials

``` /****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: POLYROOT                                            */
/*   TITLE: Computing Roots of Polynomials                      */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS: MATH                                                */
/*   PROCS: IML                                                 */
/*    DATA: Some polynomials selected from                      */
/*                                                              */
/* Woodfield, Terry J. (1988).  Simulating Stationary Gaussian  */
/*    ARMA Time Series.  Computer Science and Statistics:       */
/*    Proceedings of the 20th Symposium on the Interface,       */
/*    Reston, Virginia.                                         */
/*                                                              */
/* SUPPORT: TJW                         UPDATE: 18Sep89         */
/*     REF:                                                     */
/*                                                              */
/****************************************************************/
options ls=80 ps=60;
proc iml;
reset fuzz noname fw=8;
/*--------------------------------------------*/
/*----  Multiply two complex polynomials  ----*/
/*--------------------------------------------*/
start prodpoly(ain,bin,c);
k=nrow(ain);
m=nrow(bin);
if ncol(ain)^=2 then do;
if ncol(ain)=1 then a=ain||j(k,1,0);
else do;
print 'ERROR: Factor matrix must be N by 1 or N by 2.';
return;
end;
end;
else a=ain;
if ncol(bin)^=2 then do;
if ncol(bin)=1 then b=bin||j(m,1,0);
else do;
print 'ERROR: Factor matrix must be N by 1 or N by 2.';
return;
end;
end;
else b=bin;
n=k+m-1;
c=j(n,2,0);
do i=1 to n;
sumr=0;
sumi=0;
mn=max(1,i-m+1);
mx=min(i,k);
do j=mn to mx;
sumr=sumr+a[j,1]*b[i+1-j,1]-a[j,2]*b[i+1-j,2];
sumi=sumi+a[j,2]*b[i+1-j,1]+a[j,1]*b[i+1-j,2];
end;
c[i,1]=sumr;
c[i,2]=sumi;
end;
free i j k m n mn mx sumr sumi a b;
finish;  /*----  End of prodpoly  ----*/
/*--------------------------------------*/
/*----  Form polynomial from roots  ----*/
/*--------------------------------------*/
start formpoly(roots,poly,tol);
if nrow(tol)=0 then tol=1e-10;
if ncol(roots)^=2 then do;
if ncol(roots)=1 then rts=roots||j(nrow(roots),1,0);
else do;
print 'ERROR: ROOTS matrix must be N by 2.';
return;
end;
end;
else rts=roots;
n=nrow(rts);
poly1=j(2,2,0);
poly2=j(2,2,0);
poly2[1,1]=1;
poly2[1,2]=0;
poly2[2,1]=-rts[1,1];
poly2[2,2]=-rts[1,2];
do i=2 to n;
poly1[1,1]=1;
poly1[1,2]=0;
poly1[2,1]=-rts[i,1];
poly1[2,2]=-rts[i,2];
run prodpoly(poly1,poly2,c);
poly2=c;
end;
n=n+1;
do i=1 to n;
if abs(poly2[i,2])>tol then do;
print
'ERROR: Complex ROOTS must come in conjugate pairs.';
return;
end;
end;
poly=(poly2[,1])`;
free i n poly1 poly2 tol rts;
finish;  /*----  End of formpoly  ----*/
/*-------------------------------*/
/*----  Evaluate polynomial  ----*/
/*-------------------------------*/
start evalpl(poly,arg,result);
n=max(nrow(poly),ncol(poly));
nargs=nrow(arg);
result=j(nargs,2,0);
do i=1 to nargs;
real=arg[i,1];
imag=arg[i,2];
preal=poly[n];
pimag=0.;
preal=preal+real*poly[n-1] ;
pimag=pimag+imag*poly[n-1] ;
do j=n-2 to 1 by -1;
creal=real;
real=real*arg[i,1]-imag*arg[i,2] ;
imag=imag*arg[i,1]+creal*arg[i,2] ;
preal=preal+real*poly[j] ;
pimag=pimag+imag*poly[j] ;
end;
result[i,1]=preal;
result[i,2]=pimag;
end;
print 'Results of evaluating polynomial with coefficients';
if nrow(poly)>ncol(poly) then print (poly`);
else print poly;
tabl=arg||result||j(nargs,2,0);
do i=1 to nargs;
tabl[i,5]=arg[i,1]*arg[i,1]+arg[i,2]*arg[i,2];
tabl[i,6]=result[i,1]*result[i,1]+result[i,2]*result[i,2];
end;
cname={'R(ARG)' 'I(ARG)' 'R(P(ARG))' 'I(P(ARG))'
'|ARG|**2' '|P(ARG)|**2'};
print tabl [colname=cname];
free i j real imag preal pimag creal tabl cname;
finish;  /*----  End of evalpl  ----*/
/*------------------------------------------*/
/*----  Reverse order of vector values  ----*/
/*------------------------------------------*/
start flip(a,b);
n=ncol(a);
b=a[-(-n:-1)];
if nrow(b)>1 then b=b`;
free n;
finish;  /*----  End of flip  ----*/
/*---------------------------------*/
/*----  MAIN Calculation area  ----*/
/*---------------------------------*/
/*----  Roots: 0.9, -1.3  ----*/
roots={0.9 0, -1.3 0};
run formpoly(roots,phi,tol);
b=polyroot(phi);
run evalpl(phi,b,result);
/*----  Roots: 2.4 +/- 5i  ----*/
roots={2.4 5, 2.4 -5};
run formpoly(roots,theta,tol);
theta=3*theta;
b=polyroot(theta);
run evalpl(theta,b,result);
/*----  Roots: 0.9, -1.3, 2.4 +/- 5i  ----*/
roots={0.9 0, -1.3 0, 2.4 5, 2.4 -5};
run formpoly(roots,phi,tol);
phi=3*phi;
b=polyroot(phi);
run evalpl(phi,b,result);
/*---------------------------------------------------------*/
/*----  The following six models are taken from the    ----*/
/*----  Woodfield(1988) paper, "Simulating Univariate  ----*/
/*----  Gaussian ARMA Time Series".                    ----*/
/*---------------------------------------------------------*/
/*----  Model 1  ----*/
phi={1 -0.8 0.65};
run flip(phi,a);
b=polyroot(a);
run evalpl(a,b,r);
/*----  Model 2  ----*/
theta={1 1.25 0.35};
run flip(theta,a);
b=polyroot(a);
run evalpl(a,b,r);
/*----  Model 3  ----*/
phi={1 -0.95};
run flip(phi,a);
b=polyroot(a);
run evalpl(a,b,r);
theta={1 0.85};
run flip(theta,a);
b=polyroot(a);
run evalpl(a,b,r);
/*----  Model 4  ----*/
phi={1 -1.5 1.21 -0.455};
run flip(phi,a);
b=polyroot(a);
run evalpl(a,b,r);
theta={1 0.2 0.9};
run flip(theta,a);
b=polyroot(a);
run evalpl(a,b,r);
/*----  Model 5  ----*/
phi={1 -1 0.24};
run flip(phi,a);
b=polyroot(a);
run evalpl(a,b,r);
theta={1 0.4 0.2 0.1};
run flip(theta,a);
b=polyroot(a);
run evalpl(a,b,r);
/*----  Model 6  ----*/
phi={1 -0.3357 0.0821 0.1570 0.2567};
run flip(phi,a);
b=polyroot(a);
run evalpl(a,b,r);
theta={1 -0.6077 0.0831 0.1903};
run flip(theta,a);
b=polyroot(a);
run evalpl(a,b,r);
quit;
/*-------------------------------*/
/*----  End of POLYROOT.SAS  ----*/
/*-------------------------------*/

```