-rw-r--r-- 5249 mceliece-sage-20221023/echelon.sage raw
def 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()