Code covered by the BSD License

# gfnull

### Mark Wilde (view profile)

simple routine to find the null space of a matrix over gf (2)

gfrref(A,tol)
```% this is a modification of the standard file "rref.m" that puts a matrix
% into row-reduced echelon form. This modification acts on matrices over
% gf(2).

function [A,jb] = gfrref(A,tol)
%RREF   Reduced row echelon form.
%   R = RREF(A) produces the reduced row echelon form of A.
%
%   [R,jb] = RREF(A) also returns a vector, jb, so that:
%       r = length(jb) is this algorithm's idea of the rank of A,
%       x(jb) are the bound variables in a linear system, Ax = b,
%       A(:,jb) is a basis for the range of A,
%       R(1:r,jb) is the r-by-r identity matrix.
%
%   [R,jb] = RREF(A,TOL) uses the given tolerance in the rank tests.
%
%   Roundoff errors may cause this algorithm to compute a different
%   value for the rank than RANK, ORTH and NULL.
%
%   Class support for input A:
%      float: double, single
%

%   Copyright 1984-2004 The MathWorks, Inc.
%   \$Revision: 5.9.4.2 \$  \$Date: 2004/04/10 23:30:09 \$

[m,n] = size(A);

% Does it appear that elements of A are ratios of small integers?
% [num, den] = rat(A);
% rats = isequal(A,num./den);

% Compute the default tolerance if none was provided.
% if (nargin < 2), tol = max(m,n)*eps(class(A))*norm(A,'inf'); end
tol = .00000000000001;
% Loop over the entire matrix.
i = 1;
j = 1;
jb = [];
while (i <= m) && (j <= n)
% Find value and index of largest element in the remainder of column j.
tmp1 = A(i:m,j);
[p,k] = max(tmp1.x); k = k+i-1;
if (p <= tol)
% The column is negligible, zero it out.
A(i:m,j) = zeros(m-i+1,1);
j = j + 1;
else
% Remember column index
jb = [jb j];
% Swap i-th and k-th rows.
A([i k],j:n) = A([k i],j:n);
% Divide the pivot row by the pivot element.
A(i,j:n) = A(i,j:n)/A(i,j);
% Subtract multiples of the pivot row from all the other rows.
for k = [1:i-1 i+1:m]
A(k,j:n) = A(k,j:n) - A(k,j)*A(i,j:n);
end
i = i + 1;
j = j + 1;
end
end

A = gf(A,2);

% Return "rational" numbers if appropriate.
% if rats
%   [num,den] = rat(A);
%   A=num./den;
% end
```