def GaussJordan(A, k=-1, d0=None):
  if d0 == None:
    d0 = A.base_ring()(1)
  n = A.nrows()
  m = A.ncols()
  I = MatrixSpace(A.base_ring(),A.nrows(),A.nrows())(1)
  for i in xrange(k+1,n):
    if any(A[i,j] for j in xrange(m)):
      break
  else:
      U,P,r,h,d = d0*I, I, 0, n-k, d0
      return (U,P,r,h,d)

  if m == 1:
    i = min([i for i in xrange(k+1,n) if A[i,0] != 0])
    P = copy(I)
    P.swap_rows(i,k+1)
    r,h,d = 1, n-i, (P*A)[k+1,0]
    U = d*I
    for j in range(n):
      U[j,k+1] = -(P*A)[j,0]
    U[k+1,k+1] = d0
  else:
     m1,m2 = m//2, m-m//2
     A1 = A.matrix_from_columns(range(m1))
     B = A.matrix_from_columns(range(m1,m))
     U1, P1, r1, h1, d1 = GaussJordan(A1, k, d0)
     A2 = d0**(-1)*U1*P1*B
     U2, P2, r2, h2, d = GaussJordan(A2, k+r1, d1)
     U = d1**(-1) * U2*(P2*(U1 - d1*I) + d1*I)
     P,r,h,d = P2*P1, r1+r2, min(h1,h2),d
  return U,P,r,h,d