macro( Nullspace = `Matrix/Nullspace` );

#
#--> Nullspace(M,B,'nullity')
#
# Computes the nullspace (kernel) of the matrix A in M over a field
# Outputs the nullspace as a set of vectors and optionally the nullity
#
# Author: MBM 1990
#

Nullspace := proc(M,B,nullity)
local A,N,R,l,m,n,i,j,r,s,t,u,v;

    R := M[CoefficientRing];
    if not hasCategory(R,Field) then
	ERROR(`Coefficient ring must be a field`) fi;

    if nargs > 3 and not type(nullity,name) then
        ERROR(`2nd optional argument must be a name`) fi;

    A := M[GaussJordan](B,r); # reduced row Echelon form
    l := M[Rows](A);
    m := M[Cols](A);
    n := m-r;
    if nargs > 2 then nullity := n fi;
    if n=0 then RETURN( {} ) fi;

    A := M[ArrayCoeffs](A);
    if r > 0 then s := array(1..r) fi;
    u := array(1..n); # vector of column pointers
    v := array(1..n);
    j := 1;
    for i to l while j <= m do
	while j <= m and A[i,j] = R[0] do
	    u[j-i+1] := j;
	    v[j-i+1] := i;
	    j := j+1
	od;
	if i <= r then s[i] := j fi;
	j := j+1
    od;

    N := array(1..n);
    t := array(1..m);
    for i to n do
	# N[i], the i'th nullspace vector, is obtained from col(M,u[i])
	for j to m do t[j] := R[0] od;
	t[u[i]] := R[1];
	for j to v[i]-1 do t[s[j]] := R[`-`](A[j,u[i]]) od;
	N[i] := [seq(t[j],j=1..m)];
    od;

    [ seq(N[i], i=1..n) ]

end:

save `NS.m`;
quit
