code: def gauss_jordan(A, x, b):
      
      n = len(x)
      m = len(b)
      assert len(A) == m and len(A[0]) == n
      S = []                        # put linear system in a single matrix S
      for i in range(m):
          S.append(A[i][:] + [b[i]])
      S.append(list(range(n)))      # indices in x
      k = diagonalize(S, n, m)
      if k < m:
          for i in range(k, m):
              if not is_zero(S[i][n]):
                  return GJ_ZERO_SOLUTIONS
      for j in range(k):
          x[S[m][j]] = S[j][n]
      if k < n:
          for j in range(k, n):
              x[S[m][j]] = 0
          return GJ_SEVERAL_SOLUTIONS
      return GJ_SINGLE_SOLUTION
docs: Linear equation system Ax=b by Gauss-Jordan
  
      :param A: n by m matrix
      :param x: table of size n
      :param b: table of size m
      :modifies: x will contain solution if any
      :returns int:
            0 if no solution,
            1 if solution unique,
            2 otherwise
      :complexity: :math:`O(n^2m)`