```-rw-r--r-- 5249 mceliece-sage-20221023/echelon.sage rawdef echelon(X):
r'''
Return X reduced to row-echelon form.

INPUT:

"X" - a matrix
'''
return copy(X.rref())

def echelon_positions(R):
r'''
If R is in reduced row-echelon form,
return the list c such that
len(c) is the rank of R and, for each j in {0,1,...,r-1},
c[j] is the position of the leading 1 in row j of R.

If R is not in reduced row-echelon form,
raise an exception.

INPUT:

"R" - a matrix over a field
'''

r = R.rank()
c = [min(b for b,Rab in enumerate(R[a]) if Rab != 0) for a in range(r)]

# "c_0 < c_1 < ... < c_{r-1}"
assert c == sorted(set(c))

for a in range(r):
# "row 0 of R begins with a 1 in column c_0"
# "row 1 of R begins with a 1 in column c_1"
# etc.
assert R[a][:c[a]] == 0
assert R[a,c[a]] == 1

for a in range(r):
for b in range(R.nrows()):
if b != a:
# "this is the only nonzero entry in column c_0"
# "the only nonzero entry in column c_1"
# etc.
assert R[b,c[a]] == 0

for a in range(r,R.nrows()):
# "all subsequent rows of R are 0"
assert R[a] == 0

return c

def is_systematic(R):
r'''
Return True if R is in systematic form,
else False.

INPUT:

"R" - a matrix in reduced row-echelon form over a field
'''
c = echelon_positions(R)
r = len(c)

# "R has exactly zero rows, i.e., there are no zero rows"
if R.nrows() != r: return False

# "c_i = i for 0 <= i < r"
okc = all(c[i] == i for i in range(r))

# "equivalent to simply saying c_{r-1} == r-1,
#  except in the degenerate case r = 0"
if r > 0: assert okc == (c[r-1] == r-1)

if not okc: return False

# "in other words, R has the form (I_r | T),
#  where I is an rxr identity matrix."
for a in range(r):
for b in range(r):
assert R[a,b] == (a == b)

return True

def is_semi_systematic(R,mu,nu):
r'''
Return True if R is in systematic form,
else False.

INPUT:

"R" - a matrix in reduced row-echelon form over a field

"mu" - a nonnegative integer

"nu" - an integer with nu >= mu
'''
c = echelon_positions(R)
r = len(c)
assert mu >= 0
assert nu >= mu

# "r >= mu"
if r < mu: return False

# "there are at least r-mu+nu columns"
if R.ncols() < r-mu+nu: return False

# "R has r rows"
if R.nrows() != r: return False

# "c_i = i for 0 <= i < r-mu"
# c_i <= i-mu+nu for 0 <= i < r"
okc = all(c[i] == i for i in range(r-mu)) and all(c[i] <= i-mu+nu for i in range(r))

# "equivalent to simply c_{r-mu-1} = r-mu-1 and c_{r-1} <= r-mu+nu-1
#  except in the degenerate case r = mu"
if r > mu: assert okc == ((c[r-mu-1] == r-mu-1) and (c[r-1] <= r-mu+nu-1))

return okc

# ----- miscellaneous tests

def test_manual():
print('echelon misc')
sys.stdout.flush()

X = matrix(GF(2),[[1,1,0],[0,0,1]])
assert X == echelon(X)
assert echelon_positions(X) == [0,2]
assert echelon(matrix(GF(2),[[0,0,1],[1,1,1]])) == X
assert echelon(matrix(GF(2),[[1,1,1],[1,1,0]])) == X

X = matrix(GF(2),[[1,0,1],[0,1,1]])
assert echelon(X) == X
assert echelon_positions(X) == [0,1]
assert is_systematic(X)

X = matrix(GF(2),[[1,1,0,1],[0,0,1,1]])
assert echelon(X) == X
assert echelon_positions(X) == [0,2]
assert not is_semi_systematic(X,0,0)
assert not is_semi_systematic(X,1,1)
assert not is_semi_systematic(X,0,1)
assert is_semi_systematic(X,1,2)
assert not is_semi_systematic(X,0,2)
assert is_semi_systematic(X,1,3)

def test_smallrandom():
for q in range(100):
q = ZZ(q)
if not q.is_prime(): continue
k = GF(q)

nummatrices = 0
numsystematic = 0
numsemisystematic = 0
for rows in range(10):
for cols in range(10):
nummatrices += 1
X = matrix(k,[[randrange(q) for b in range(cols)] for a in range(rows)])
R = echelon(X)
assert R.row_space() == X.row_space()
assert R.rank() == X.rank()
c = echelon_positions(R)
assert len(c) == R.rank()
if len(c) == R.nrows() and (len(c) == 0 or c[-1] == len(c)-1):
numsystematic += 1
assert is_systematic(R)
assert matrix([Ra[:len(c)] for Ra in R]) == identity_matrix(len(c))
else:
assert not is_systematic(R)
assert matrix([Ra[:len(c)] for Ra in R]) != identity_matrix(len(c))
if len(c) == R.nrows():
if len(c) == 0:
assert is_semi_systematic(R,0,0)
else:
assert is_semi_systematic(R,len(c),c[-1]+1)
print('echelon field %d nummatrices %d numsystematic %d' % (q,nummatrices,numsystematic))
sys.stdout.flush()

for mu in range(5):
for nu in range(mu,mu+5):
for rows in range(mu,mu+5):
if rows == 0: continue
for cols in range(rows-mu+nu,rows-mu+nu+5):
X = matrix(k,[[a==b if b<rows else randrange(q) for b in range(cols)] for a in range(rows)])
perm = list(range(nu))
shuffle(perm)
perm = [b for b in range(rows-mu)]+[rows-mu+b for b in perm]+[b for b in range(rows-mu+nu,cols)]
assert len(perm) == cols
X = X.matrix_from_columns(perm)
assert is_semi_systematic(echelon(X),mu,nu)

if __name__ == '__main__':
test_manual()
test_smallrandom()
```