NullSpaceBasis := proc(A::Matrix, p::integer)
  local M,n,k,i,mult,j,zero,v;

  # The parameter "p" specifies that all arithmetic is to be done in Zp .

  # Given a square matrix A, we return a basis [ v[1], ..., v[k] ] for
  # the null space { v : v.A = 0 } of A.  The algorithm does this by
  # transforming the matrix to triangular idempotent form.

  M := Matrix(A);

  use LinearAlgebra in
    n := RowDimension(M);
    if ColumnDimension(M) <> n then 
      error "expecting a square matrix" 
    end if;

    for k from 1 to n do
      # Search for pivot element.
      for i from k to n while M[k,i] = 0 do i := i+1 end do;

      if i <= n then
        # Normalize column i and interchange this with column k.
        mult := M[k,i]^(-1) mod p;
        M := ColumnOperation(M,i,mult);
        M := map(`mod`,M,p);
        M := ColumnOperation(M,[i,k]);

        # Eliminate rest of row k via column operations.
        for i from 1 to n do
          if i <> k then
            M := ColumnOperation(M,[i,k],-M[k,i]);
            M := map(`mod`,M,p)
          end if
        end do
      end if
    end do;

    # Convert M to I-M.
    M := ScalarMultiply(M,-1);
    for i from 1 to n do M[i,i] := 1 + M[i,i] end do;

    # Read off nonzero rows of M.
    i := 0;  j := 1;
    zero := Vector[row](n);
    while j <= n do
      while j <= n and Equal(Row(M,j),zero) do j := j+1 end do;
      if j <= n then
        i := i+1;
        v[i] := Row(M,j)
      end if;
      j := j+1
    end do
  end use;

  # Return the list of row vectors representing the null space basis.
  [ seq(v[k], k=1..i) ]

end proc: