refmod(A, p)

function A = refmod(A, p)
% REFMOD    Row echelon form modulo P
%    A = REFMOD(A, P) puts A in row echelon form using gaussian 
%    elimination mod P

% by David Szotten 2009
% Follows http://en.wikipedia.org/wiki/Gaussian_elimination#Pseudocode
  
  [m, n] = size(A);
  
  i = 1;
  j = 1;
  
  A = mod(A, p);
  
  while (i <= m && j <= n) 
    
    % Find pivot in column j, starting in row i:
    [tmp, maxi] = max( A( i:end, j) );
    
    % we started counting from row i
    maxi = maxi + i - 1;
    
    if( A( maxi, j) ~= 0)
      % swap rows i and maxi, but do not change the value of i
      % Now A[i,j] will contain the old value of A[maxi,j].
      rows = 1:m;
      rows(i) = maxi;
      rows(maxi) = i;
      
      A = A(rows,:);
      
      % divide each entry in row i by A[i,j]
      % Now A[i,j] will have the value 1.
      A(i, :) = mod( A(i, :) * invMod( A(i,j), p), p);
      
      % subtract A[u,j] * row i from row u
      % Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
      A(i+1:end, :) = A(i+1:end, :) - A(i+1:end, j) * A(i,:);
      A(i+1:end, :) = mod(  A(i+1:end, :) , p);
	   
      i = i + 1;
    end
    j = j + 1;
  end
end
  
  
  
function y = invMod(x,p)
% inverse of x mod p  

  [g,c,d] = gcd(x,p); 
  y = mod(c,p);
end