% --------------------------

% A -> Gauss column-reduction of A + index vector
% A = array of its columns
% A must have at least one column
% modifies A itself

/gauss { 1 dict begin 
/A exch def
/r 0 def
/c 0 def
/C A length def         % C = number of columns
/R A 0 get length def   % R = number of rows
/index [ 0 1 C 1 sub {
} for
] def
% index = index array [ 0 1 2 ... C-1 ]
/PR C array def
% PR[c] = pivot row for column c
{
  r R ge 
  c C ge 
  or { exit } if

  % r < R and c < C
  % find pivot entry in row r
  /pc c def  % pc = pivot column
  /max A c get r get abs def
  c 1 add 1 C 1 sub { /j exch def
    /x A j get r get abs def
    x max gt {
      /pc j def
      /max x def
    } if
  } for
  % now pc = the pivot column
  max 0 gt {
    % swap entries pc and c in i 
    /tmp index c get def
    index c index pc get put
    index pc tmp put
    % and columns pc and c in A
    /tmp A c get def
    A c A pc get put
    A pc tmp put

    % pivot
    /c0 A c get def  % the pivoting column
    /p0 c0 r get def % the pivot value
    c 1 add 1 C 1 sub { /j exch def
      /c1 A j get def
      /p c1 r get p0 div def
      c1 r p put
      % subtract (p1/p0)*column c from column j
      r 1 add 1 R 1 sub { /i exch def
        c1 
        dup i exch % c1 i c1
        i get c0 i get p mul sub put
      } for
    } for
    PR c r put
    /r r 1 add def
    /c c 1 add def
  }{
    % all zero entries
    /c c 1 add def
    % continue loop with same r
  } ifelse
} loop
/B [
  0 1 c 1 sub { /i exch def
    0 1 r 1 sub { /j exch def
      
    } for
  } for
] def
/B [

] def
[
  /c 0 def
  PR {
    % r(c)
    
    /c c 1 add def
  } forall
] ==
A index PR
% the matrix is singular 
% when the last entry in PR is null
end } def

% --------------------------

% A -> copy 
/clone { 1 dict begin 
[ exch {
  [ exch {
  } forall
  ]
} forall
]
end } def

/matrix-mul { 1 dict begin
/B exch def
/A exch def
/r A 0 get length def
/c A length def
[
% replace each column of B by a column of length r
B { /b exch def
  [
  0 1 r 1 sub { /i exch def
    0 
    0 1 c 1 sub { /j exch def
      b j get 
      A j get i get mul 
      add
    } for
  } for
  ]
} forall
]
end } def


% ----------------------------------------------------

/A [[1 0][1 0]] def

A gauss [ 4 1 roll ] ==

quit