`mod/rtable/Gausselim` := proc(A::Matrix,rank::{integer,name},
                               det::{integer,name})
  local    B, d, nrows, ncols, right_margin, i,j,
      col, row, s, t, p, Q;

  Q := type( A, 'Matrix( rational )' );
  if not Q and not type( A, 'Matrix( ratpoly( algnum ) )' ) then
      error "first argument must be a matrix over a finite field"
  end if;
  ( nrows, ncols ) := LinearAlgebra:-Dimensions( A );

  p := args[-1];
  right_margin := ncols;
  if type( args[-2], 'integer' ) then
      right_margin := args[-2]
  end if;

  if nargs > 3 and type( det, 'name' ) 
               and nrows <> right_margin then
    error "cannot compute determinant of a non-square matrix"
  end if;

  if Q then
      B := map( s -> s mod p,  A );
  else
      B := map( s -> 'Normal'( s ) mod p, A );
  end if;

  row := 1; # initialise the row
  d := 1; # initialise the determinant

  for col to right_margin while row <= nrows do
      # Find the first nonzero entry in this column
      for i from row to nrows while B[ i, col ] = 0 do
          # empty
      end do;
      if i > nrows then
          # found a zero column so determinant is zero
          d := 0;
          next
      end if;
      if i <> row then
          # interchange row i with row row
          B[ i, 1..-1 ], B[ row, 1..-1 ] := 
              B[ row, 1..-1 ], B[ i, 1 .. -1 ];
          d := -d mod p; # keep track of determinant
      end if;

      if nargs = 4 then
          d := 'Normal'( d * B[ row, col ] ) mod p
      end if;

      # Make all the nonzero entries below this one 
      # are equal to 1
      if Q then
          s := 1/B[ row, col ] mod p;
          for i from row + 1 to nrows do
              # do i-th row
              if B[ i, col ] = 0 then
                  next
              end if;
              t := s * B[ i, col ] mod p;
              for j from col + 1 to ncols do
                  B[ i, j ] := B[ i, j ] - 
                               t * B[ row, j ] mod p
              end do;
              B[ i, col ] := 0
          end do
      else
          s := 'Normal'( 1 / B[ row, col ] ) mod p;
          for i from row + 1 to nrows do
              if B[ i, col ] = 0 then
                  next
              end if;
              t := 'Normal'( s * B[ i, col ] ) mod p;
              for j from col + 1 to ncols do
                B[i,j] := 'Normal'( B[ i, j ] - 
                                    t * B[ row, j ] ) mod p
              end do;
              B[ i, col ] := 0
          end do
      end if;

      row := 1 + row # go to next row
  end do; # go to next column

  if nargs > 2 and type( rank, 'name' ) then
      rank := row-1
  end if;
  if nargs > 3 and type( det, 'name' ) then
      det := d
  end if;
  B
end proc: