From 33dd3950a3fe7afada7359a68c6a9537ddd195d6 Mon Sep 17 00:00:00 2001 From: Edgar Costa Date: Thu, 25 Mar 2021 12:42:03 -0400 Subject: [PATCH 1/6] comment --- src/sage/matrix/matrix_nmod_dense.pyx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/sage/matrix/matrix_nmod_dense.pyx b/src/sage/matrix/matrix_nmod_dense.pyx index 7f87ba7de87..bbed69174f7 100644 --- a/src/sage/matrix/matrix_nmod_dense.pyx +++ b/src/sage/matrix/matrix_nmod_dense.pyx @@ -216,6 +216,8 @@ cdef class Matrix_nmod_dense(Matrix_dense): self.cache('rank', rank) def echelon_form(self): + #FIXME: add a warning regarding the output potentially having more rows + # call howell_form and key='echelon_form' ans = self.fetch(key) if ans is not None: From c1961a223c074ea026d3ecff1ff076fa4a6bdd54 Mon Sep 17 00:00:00 2001 From: Edgar Costa Date: Sun, 28 Mar 2021 13:03:54 -0400 Subject: [PATCH 2/6] right kernel matrix --- src/sage/matrix/matrix_nmod_dense.pyx | 119 ++++++++++++++++++++------ 1 file changed, 93 insertions(+), 26 deletions(-) diff --git a/src/sage/matrix/matrix_nmod_dense.pyx b/src/sage/matrix/matrix_nmod_dense.pyx index ca58437ba99..035e20ed358 100644 --- a/src/sage/matrix/matrix_nmod_dense.pyx +++ b/src/sage/matrix/matrix_nmod_dense.pyx @@ -12,6 +12,7 @@ from cpython.sequence cimport * from cysignals.signals cimport sig_on, sig_str, sig_off +from sage.arith.all import gcd from sage.arith.power cimport generic_power from sage.arith.long cimport integer_check_long_py from sage.structure.sage_object cimport SageObject @@ -45,13 +46,11 @@ cdef class Matrix_nmod_dense(Matrix_dense): def __init__(self, parent, entries=None, copy=None, bint coerce=True): self._parent = parent # MatrixSpace over IntegerMod_int or IntegerMod_int64 ma = MatrixArgs_init(parent, entries) - cdef long z for t in ma.iter(coerce, True): #???? se = t z = se.entry nmod_mat_set_entry(self._matrix, se.i, se.j, z) - nmod_mat_print_pretty(self._matrix) def __dealloc__(self): nmod_mat_clear(self._matrix) @@ -293,6 +292,7 @@ cdef class Matrix_nmod_dense(Matrix_dense): # Extra + #FIXME returning wrong type def tranpose(self): cdef Matrix_nmod_dense M = self._new(self._nrows, self._ncols) sig_on() @@ -302,26 +302,27 @@ cdef class Matrix_nmod_dense(Matrix_dense): def echelonize(self): + #FIXME: add a warning regarding the output potentially having more rows than input """ Echelon form in place """ - if not self._parent._base.is_field(): - raise NotImplementedError("Only implemented over fields") self.check_mutability() self.clear_cache() - rank = nmod_mat_rref(self._matrix) - self.cache('rank', rank) + if self._parent._base.is_field(): + rank = nmod_mat_rref(self._matrix) + self.cache('rank', rank) + else: + self._howell_form() def echelon_form(self): - #FIXME: add a warning regarding the output potentially having more rows - # call howell_form and + #FIXME: add a warning regarding the output potentially having more rows than input key='echelon_form' ans = self.fetch(key) if ans is not None: return ans ans = self.__copy__() - self.cache(key, ans) ans.echelonize() + self.cache(key, ans) return ans @@ -331,13 +332,14 @@ cdef class Matrix_nmod_dense(Matrix_dense): if ans is not None: return ans cdef Matrix_nmod_dense E + cdef Py_ssize_t i, j, k # howell form has all the zero rows at the bottom - E = self.echelon_form() if self._parent._base.is_field() else self.howell_form() + E = self.echelon_form() p = [] zdp = [] k = 0 - for i from 0 <= i < E._nrows: - for j from k <= j < E._ncols: + for i in range(E._nrows): + for j in range(k, E._ncols): if nmod_mat_get_entry(E._matrix, i, j) != 0: # nonzero position if nmod_mat_get_entry(E._matrix, i, j) == 1: p.append(j) @@ -358,10 +360,17 @@ cdef class Matrix_nmod_dense(Matrix_dense): ans = self.fetch(key) if ans is not None: return ans - if self._parent._base.is_field(): - sig_on() - ans = nmod_mat_rank(self._matrix) - sig_off() + if not self._parent._base.is_field(): + p = self.pivots() + ans = len(p[0]) + else: + p = self.fetch('pivots') + if ans is not None: + ans = len(p[0]) + else: + sig_on() + ans = nmod_mat_rank(self._matrix) + sig_off() self.cache(key, ans) return ans @@ -412,7 +421,7 @@ cdef class Matrix_nmod_dense(Matrix_dense): """ In place strong echelon form of self """ - if self._nrows >= self._ncols: + if self._nrows < self._ncols: raise ValueError("self must must have at least as many rows as columns.") self.check_mutability() self.clear_cache() @@ -421,6 +430,7 @@ cdef class Matrix_nmod_dense(Matrix_dense): sig_off() def strong_echelon_form(self): + #FXIME add warning """ Strong echelon form of self """ @@ -428,7 +438,11 @@ cdef class Matrix_nmod_dense(Matrix_dense): ans = self.fetch(key) if ans is not None: return ans - ans = self.__copy__() + if self._nrows < self._ncols: + M = self._new(self._ncols - self._nrows, self._ncols) + ans = self.stack(M) + else: + ans = self.__copy__() ans._strong_echelon_form() self.cache(key, ans) return ans @@ -437,13 +451,17 @@ cdef class Matrix_nmod_dense(Matrix_dense): """ In place Howell form of self """ - if self._nrows >= self._ncols: - raise ValueError("self must must have at least as many rows as columns.") + if self._nrows < self._ncols: + raise ValueError("self must have at least as many rows as columns.") self.check_mutability() + cdef Matrix_nmod_dense M self.clear_cache() + sig_on() nmod_mat_howell_form(self._matrix) + sig_off() def howell_form(self): + #FXIME add warning """ Howell form of self """ @@ -451,7 +469,12 @@ cdef class Matrix_nmod_dense(Matrix_dense): ans = self.fetch(key) if ans is not None: return ans - ans = self.__copy__() + + if self._nrows < self._ncols: + M = self._new(self._ncols - self._nrows, self._ncols) + ans = self.stack(M) + else: + ans = self.__copy__() ans._howell_form() self.cache(key, ans) return ans @@ -505,21 +528,65 @@ cdef class Matrix_nmod_dense(Matrix_dense): return M def _right_kernel_matrix(self): - cdef Matrix_nmod_dense X, ans, echelon_form; + cdef Matrix_nmod_dense X, ans, E + cdef Py_ssize_t i, j, k, l + cdef long x, y, s + E = self.echelon_form() if self._parent._base.is_field(): # nmod_mut_nullspace will do this regardless # so we are better off to start in echelon form to have the rank - echelon_form = self.echelon_form() X = self._new(self.ncols, self._nrows - self.rank()) ans = self._new(self._nrows - self.rank(), self._ncols) sig_on() - nmod_mat_nullspace(X._matrix, echelon_form._matrix) # columns of X form a basis + nmod_mat_nullspace(X._matrix, E._matrix) # columns of X form a basis nmod_mat_transpose(ans._matrix, X._matrix) sig_off() return ans else: - #I'm here - strong_echelon_form = self.strong_echelon_form() + zero = self._parent._base.zero() + one = self._parent._base.one() + p, zdp = map(set, self._pivots()) + set_pivots = p.union(zdp) + pivot = sorted(list(set_pivots)) + basis = [] + k = 0 + for j in range(self._ncols): + if j in p: + k += 1 + continue + v = [zero] * self._ncols + i = k + if j in zdp: + k += 1 + v[j] = self._modulus.int64/nmod_mat_get_entry(E._matrix, i, j) + else: + v[j] = one + + # figure out the remaining coefficients of v + # note that v might need to be rescaled several times + for l in reversed(range(i)): + x = nmod_mat_get_entry(E._matrix, l, pivot[l]) + # solve + # v[pivot[l]] E[l, pivot[l]] + y*sum(E[l,m] * v[m] for m in range(pivot[l] + 1, j)) = 0 + s = sum(v[m]*nmod_mat_get_entry(E._matrix, l, m) for m in range(pivot[l] + 1, j + 1)) % self._modulus.int64 + if s % x != 0: # make sure we can work mod N/x + y = x//gcd(s, x) + s *= y # now s is divisible by x + for m in range(pivot[l] + 1, j + 1): + v[m] *= y + assert v[j] % self._modulus.int64 != 0 + # QUESTION: this is correct modulo N/x, does one need to consider the various lifts? + # FIXME, this feels wrong + v[pivot[l]] = self._parent._base((-s//x)) + basis.append(v) + # FIXME, this feels wrong + ans = self._new(len(basis), self._ncols) + ma = MatrixArgs_init(ans._parent, basis) + for t in ma.iter(False, True): #???? + se = t + x = se.entry + nmod_mat_set_entry(ans._matrix, se.i, se.j, x) + return ans # random matrix generation (David) From b9965be799b4612afbdec8f8b75e05a6037a173d Mon Sep 17 00:00:00 2001 From: Edgar Costa Date: Sun, 28 Mar 2021 14:27:03 -0400 Subject: [PATCH 3/6] remove comment --- src/sage/matrix/matrix_nmod_dense.pyx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/sage/matrix/matrix_nmod_dense.pyx b/src/sage/matrix/matrix_nmod_dense.pyx index b33156c2ac4..f83ad118451 100644 --- a/src/sage/matrix/matrix_nmod_dense.pyx +++ b/src/sage/matrix/matrix_nmod_dense.pyx @@ -492,7 +492,6 @@ cdef class Matrix_nmod_dense(Matrix_dense): # Extra - #FIXME returning wrong type def tranpose(self): cdef Matrix_nmod_dense M = self._new(self._nrows, self._ncols) sig_on() From 10452071cf97d5f4c85ca926c46e7dd1dfa0ca00 Mon Sep 17 00:00:00 2001 From: Edgar Costa Date: Mon, 29 Mar 2021 13:25:01 -0400 Subject: [PATCH 4/6] touch up right kernel --- src/sage/matrix/matrix_nmod_dense.pyx | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix_nmod_dense.pyx b/src/sage/matrix/matrix_nmod_dense.pyx index f83ad118451..22c5923d2f9 100644 --- a/src/sage/matrix/matrix_nmod_dense.pyx +++ b/src/sage/matrix/matrix_nmod_dense.pyx @@ -415,7 +415,9 @@ cdef class Matrix_nmod_dense(Matrix_dense): algorithm = "crt" if algorithm == "FLINT": f = R[var]() + sig_on() nmod_mat_charpoly(&f.x, self._matrix) + sig_off() return f if algorithm == "lift": return self.lift_centered().charpoly(var).change_ring(R) @@ -508,7 +510,9 @@ cdef class Matrix_nmod_dense(Matrix_dense): self.check_mutability() self.clear_cache() if self._parent._base.is_field(): + sig_on() rank = nmod_mat_rref(self._matrix) + sig_off() self.cache('rank', rank) else: self._howell_form() @@ -726,7 +730,7 @@ cdef class Matrix_nmod_dense(Matrix_dense): sig_off() return M - def _right_kernel_matrix(self): + def _right_kernel_matrix(self, zero_divisors_are_pivots=False): cdef Matrix_nmod_dense X, ans, E cdef Py_ssize_t i, j, k, l cdef long x, y, s @@ -757,6 +761,8 @@ cdef class Matrix_nmod_dense(Matrix_dense): i = k if j in zdp: k += 1 + if zero_divisors_are_pivots: + continue v[j] = self._modulus.int64/nmod_mat_get_entry(E._matrix, i, j) else: v[j] = one @@ -785,7 +791,10 @@ cdef class Matrix_nmod_dense(Matrix_dense): se = t x = se.entry nmod_mat_set_entry(ans._matrix, se.i, se.j, x) - return ans + if zero_divisors_are_pivots: + return ans + else: + return ans.howell_form() # random matrix generation (David) From fb619b0ab2ccdb99d85fd353c9a757bd40c16b71 Mon Sep 17 00:00:00 2001 From: Edgar Costa Date: Mon, 29 Mar 2021 17:00:05 -0400 Subject: [PATCH 5/6] work around echenolize --- src/sage/matrix/matrix_nmod_dense.pyx | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/sage/matrix/matrix_nmod_dense.pyx b/src/sage/matrix/matrix_nmod_dense.pyx index 22c5923d2f9..866d692533b 100644 --- a/src/sage/matrix/matrix_nmod_dense.pyx +++ b/src/sage/matrix/matrix_nmod_dense.pyx @@ -515,6 +515,10 @@ cdef class Matrix_nmod_dense(Matrix_dense): sig_off() self.cache('rank', rank) else: + if self._nrows < self._ncols: + M = self._new(self._ncols - self._nrows, self._ncols) + #FIXME: this doesn't work + self = M.stack(self) self._howell_form() def echelon_form(self): @@ -523,7 +527,11 @@ cdef class Matrix_nmod_dense(Matrix_dense): ans = self.fetch(key) if ans is not None: return ans - ans = self.__copy__() + if self._nrows < self._ncols: + M = self._new(self._ncols - self._nrows, self._ncols) + ans = M.stack(self) + else: + ans = self.__copy__() ans.echelonize() self.cache(key, ans) return ans From 45de1c83ffcf3f6e95192be8f78a547786c95479 Mon Sep 17 00:00:00 2001 From: Edgar Costa Date: Mon, 29 Mar 2021 17:01:31 -0400 Subject: [PATCH 6/6] remove stack --- src/sage/matrix/matrix_nmod_dense.pyx | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/sage/matrix/matrix_nmod_dense.pyx b/src/sage/matrix/matrix_nmod_dense.pyx index 866d692533b..d84db72a17e 100644 --- a/src/sage/matrix/matrix_nmod_dense.pyx +++ b/src/sage/matrix/matrix_nmod_dense.pyx @@ -515,10 +515,6 @@ cdef class Matrix_nmod_dense(Matrix_dense): sig_off() self.cache('rank', rank) else: - if self._nrows < self._ncols: - M = self._new(self._ncols - self._nrows, self._ncols) - #FIXME: this doesn't work - self = M.stack(self) self._howell_form() def echelon_form(self):