procedure SS(num:array,den:array,out A[1,1],out B[1],out C[1], out D)
var i:integer,j:integer,n:integer,m:integer,bn;
  n=cols(den)-1; m=cols(num)-1;
  for (i=1,n) begin B[i]=0; C[i]=0; for (j=1,n) A[i,j]=0; end;
  if (m=n) then bn=num[n+1] else bn=0;
  D=den[n+1];
  for (i=1,n-1) A[i,i+1]=1;
  for (i=1,n) A[n,i]=-den[i]/D;
  for (i=1,min(n,m+1)) C[i]=(num[i]+A[n,i]*bn)/D;
  B[n]=1; D=bn/D;
end;