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  ----*/
 /*-------------------------------*/