function [D,M,N] = pquick2(G,maxop)

% PQUICK2 Equivalence of matrices of polynomials.
%	[D,M,N] = PQUICK2(A) returns unimodular matrices M and N 
% 	and diagonal matrix D equivalent to A such that M*A*N = D.
%
%	A second input argument can be used to specify the
% 	maximum number of operations (default: maxop=2000).  
%
%	See also PALGO2, PALGO
  
%	Adapted from Harold M. Edwards, "Linear Algebra."
%	Written by Serge Tchikanda 06/18/96.
%       Last revised by Serge Tchikanda 07/03/96.

%=============Possible OPERATIONS===============
% 	1 corresponds to a row operation.
% 	2 corresponds to a column operation.
%===============================================
	
	
if nargin < 2
   maxop = 2000;  
end

[m,n,OPERATIONS] = ppriori(G);  	        % Determine priority
						% entry in G.
[mrows2,ncols2] = psize(G);
M = peye(1,mrows2);			
N = peye(1,ncols2);	
		
steps = 0;
while (~isempty(m)) & (~isempty(n)) & (steps <= maxop)						                
                                                                
	if OPERATIONS == 2	 
	    rc = 2;
	    k = m;
            l = n;		
	elseif OPERATIONS == 1
	    rc = 1; 
    	    k = n;
	    l = m;
	end  
	[pet,d,lc] = gpet(G,k,k);		        
        degl = gdeg(G,k,k); degr = gdeg(G,m,n);

	% ------ Case 1: Diagonal entry is zero ------- 
              
	    if  pet == 0		     	
		
		if rc == 2
  		   G = prcop(G,rc,[1 1],l,k);
		   N = prcop(N,rc,[1 1],l,k);
 		else
 		   G = prcop(G,rc,[1 1],l,k);
		   M = prcop(M,rc,[1 1],l,k);
 		end

	% ------ Case 2: Equal degree and diagonal entry is monic ------
	    elseif (degr == degl) & (ismonic(pet,d))
	
		[pet,d,lc] = gpet(G,m,n);
		if rc == 2
  		   G = prcop(G,rc,[-lc,d],k,l);
		   N = prcop(N,rc,[-lc,d],k,l);
 		else
 		   G = prcop(G,rc,[-lc,d],k,l);
		   M = prcop(M,rc,[-lc,d],k,l);
 		end

	% ------ Case 3: Equal degree and diagonal entry is not monic ------
	    elseif (degr == degl) & (~ismonic(pet,d))

		[pet,d,lc] = gpet(G,[k m],[k n]);
		if rc == 2
  		   G = prcop(G,rc,[(d-lc(1)),lc(2)],l,k);
		   N = prcop(N,rc,[(d-lc(1)),lc(2)],l,k);
 		else
 		   G = prcop(G,rc,[(d-lc(1)),lc(2)],l,k);
		   M = prcop(M,rc,[(d-lc(1)),lc(2)],l,k);
 		end

	% ------ Case 4: Different degrees ------

	    elseif (degl ~= degr)
		
		[pet,d,lc] = gpet(G,[k m],[k n]);
		if degl < degr
		   pk = pet(1,:);
		   pl = pet(2,:)*lc(1);
		   mul = [-lc(2) zeros(1,degr-degl)];
		   pl1 = conv(mul,pk);
		   npl = length(pl);
		   npl1 = length(pl1);
		   if npl < npl1
			pl = [0*(1:npl1-npl) pl];
		   elseif npl > npl1
			pl1 = [0*(1:npl-npl1) pl1];
		   end
		   pl = pl + pl1;
                  
		   if all(pl == 0)
 		      if rc == 2
  		         G = prcop(G,rc,[-sign(lc(1))*lc(2) zeros(1,degr-degl-1) abs([1 1]*lc(1))],k,l);
			 N = prcop(N,rc,[-sign(lc(1))*lc(2) zeros(1,degr-degl-1) abs([1 1]*lc(1))],k,l);
 		      else
 		         G = prcop(G,rc,[-sign(lc(1))*lc(2) zeros(1,degr-degl-1) abs([1 1]*lc(1))],k,l);
			 M = prcop(M,rc,[-sign(lc(1))*lc(2) zeros(1,degr-degl-1) abs([1 1]*lc(1))],k,l);
 		      end
		   else
 		      if rc == 2
  		         G = prcop(G,rc,[-lc(2),zeros(1,degr-degl),lc(1)],k,l);
			 N = prcop(N,rc,[-lc(2),zeros(1,degr-degl),lc(1)],k,l);
 		     else
 		         G = prcop(G,rc,[-lc(2),zeros(1,degr-degl),lc(1)],k,l);
			 M = prcop(M,rc,[-lc(2),zeros(1,degr-degl),lc(1)],k,l);
 		     end
		   end

             	elseif degl > degr
		   pk = pet(1,:)*lc(2);
		   pl = pet(2,:);
		   mul = [-lc(1) zeros(1,degl-degr)];
		   pk1 = conv(mul,pl);
		   npk = length(pk);
		   npk1 = length(pk1);
		   if npk < npk1
			pk = [0*(1:npk1-npk) pk];
		   elseif npk > npk1
			pk1 = [0*(1:npk-npk1) pk1];
		   end
		   pk = pk + pk1;

		   if all(pk == 0)
 		      if rc == 2
  		         G = prcop(G,rc,[-sign(lc(2))*lc(1) zeros(1,degl-degr-1) abs([1 1]*lc(2))],l,k);
			 N = prcop(N,rc,[-sign(lc(2))*lc(1) zeros(1,degl-degr-1) abs([1 1]*lc(2))],l,k);
 		      else
 		         G = prcop(G,rc,[-sign(lc(2))*lc(1) zeros(1,degl-degr-1) abs([1 1]*lc(2))],l,k);
			 M = prcop(M,rc,[-sign(lc(2))*lc(1) zeros(1,degl-degr-1) abs([1 1]*lc(2))],l,k);
 		      end
		   else
 		      if rc == 2
  		         G = prcop(G,rc,[-lc(1),zeros(1,degl-degr),lc(2)],l,k);
			 N = prcop(N,rc,[-lc(1),zeros(1,degl-degr),lc(2)],l,k);
 		      else
 		         G = prcop(G,rc,[-lc(1),zeros(1,degl-degr),lc(2)],l,k);
			 M = prcop(M,rc,[-lc(1),zeros(1,degl-degr),lc(2)],l,k);
		      end
 		   end
                end
        end
	if any(any(abs(G) > 99999999999999)) 
           disp(' ')
           disp('   Warning: Numbers are getting too large for')
	   disp('   this example.')
	   return
        end
   
	[m,n,OPERATIONS] = ppriori(G);
	steps = steps + 1;

        if steps > maxop
            disp(' ')
            disp('    Maximum number of operations exceeded. ')
            disp(['    ',int2str(steps),' operations performed.'])
        end
end     
	     
D = G;