poly_gcd(p,q)

function g = poly_gcd(p,q)
% *** Find polymonial GCD by "Leading-coefficient-Elimination" ***
%     F C Chang    revised   10/05/09

      n = length(p)-1;  nc = max(find(p))-1;
      m = length(q)-1;  mc = max(find(q))-1;
      nz = zeros(1,min(n-nc,m-mc));
   if nc*mc == 0, g = [1,nz];  return,  end;
      g1 = p(1:nc+1);   g2 = q(1:mc+1);
  while(1),
      g3 = [g1,zeros(1,length(g2)-length(g1))]-(g1(1)/g2(1))*  ...
           [g2,zeros(1,length(g1)-length(g2))];
   if norm(g3,inf)/norm(g1,inf) < 1.e-3,  break;  end;
   if length(g1) > length(g2),  g1 = g2;    end;
      g2 = g3(min(find(abs(g3)>1.e-10)):max(find(abs(g3)>1.e-10))); 
  end;
      g = [g1/g1(1),nz];