/* Speicher loeschen */
kill(all)$

/* Definition der Funktion rref */
rref(M):=block([r,c,i,j],
  M:echelon(M),
  [r,c]:matrix_size(M),
  if r > c then h:c else h:r,
  for i:r-1 step -1 thru 1 do
    for j:h step -1 thru i+1 do M:rowop(M,i,j,M[i,j]),
  M
)$

/* Definition der Funktion rowm (Multiplikation einer Zeile einer Matrix mit c)*/
rowm(A,n,c):=rowop(A,n,n,1-c)$


/* Definition der Funktion GtoM (Gleichungssystem zu Matrix)*/
GtoM(G,V):=block([r,c,i],
  M:augcoefmatrix(G,V),
  [r,c]:matrix_size(M),
  for i:1 thru r do M[i,c]:M[i,c]*(-1),
  M
)$

/* sucht in der Matrix M in der Spalte s, angefangen in der Zeile zmin die Zeile z, in der der Koeffizient ungleich ist. Ist die Suche nicht spätestens in Zeile zmax
erfolgreich, so wird z = 0 zurückgegeben */

zeile_mit_Koeff_ungleich_null(M,s,zmin,zmax):=block([i],
  i:zmin,
  while (M[i,s] = 0) and (i < zmax) do
    i:i+1,
  if M[i,s] # 0 then z:i else z:0,
  z
)$

/* es wird die Gauss-Jordan-Form der Matrix erstellt, wobei die letzte Spalte
als rechte Seite eines zugeordneten Gleichungssystems aufgefasst wird */
gj(M):=block([r,c,zeile,spalte,znn,z],
  [r,c]:matrix_size(M),
  zeile:1,spalte:1,
  while (zeile <= r) and (spalte <= c-1) do
    (znn:zeile_mit_Koeff_ungleich_null(M,spalte,zeile,r),
     if znn # 0
     then
       ( M:rowswap(M,zeile,z),
         M:rowm(M,zeile,1/M[zeile,spalte]),
         for z:1 thru zeile-1 do M:rowop(M,z,zeile,M[z,spalte]),
         for z:zeile+1 thru r do M:rowop(M,z,zeile,M[z,spalte]),
         M:ratsimp(M),
         spalte:spalte+1,
         zeile:zeile+1
       )
     else
       spalte:spalte+1
    ),
  m:ratsimp(M),
  M
)$