'''Functions returning normal forms of matrices''' from collections import defaultdict from .domainmatrix import DomainMatrix from .exceptions import DMDomainError, DMShapeError from sympy.ntheory.modular import symmetric_residue from sympy.polys.domains import QQ, ZZ # TODO (future work): # There are faster algorithms for Smith and Hermite normal forms, which # we should implement. See e.g. the Kannan-Bachem algorithm: # def smith_normal_form(m): ''' Return the Smith Normal Form of a matrix `m` over the ring `domain`. This will only work if the ring is a principal ideal domain. Examples ======== >>> from sympy import ZZ >>> from sympy.polys.matrices import DomainMatrix >>> from sympy.polys.matrices.normalforms import smith_normal_form >>> m = DomainMatrix([[ZZ(12), ZZ(6), ZZ(4)], ... [ZZ(3), ZZ(9), ZZ(6)], ... [ZZ(2), ZZ(16), ZZ(14)]], (3, 3), ZZ) >>> print(smith_normal_form(m).to_Matrix()) Matrix([[1, 0, 0], [0, 10, 0], [0, 0, -30]]) ''' invs = invariant_factors(m) smf = DomainMatrix.diag(invs, m.domain, m.shape) return smf def add_columns(m, i, j, a, b, c, d): # replace m[:, i] by a*m[:, i] + b*m[:, j] # and m[:, j] by c*m[:, i] + d*m[:, j] for k in range(len(m)): e = m[k][i] m[k][i] = a*e + b*m[k][j] m[k][j] = c*e + d*m[k][j] def invariant_factors(m): ''' Return the tuple of abelian invariants for a matrix `m` (as in the Smith-Normal form) References ========== [1] https://en.wikipedia.org/wiki/Smith_normal_form#Algorithm [2] https://web.archive.org/web/20200331143852/https://sierra.nmsu.edu/morandi/notes/SmithNormalForm.pdf ''' domain = m.domain if not domain.is_PID: msg = "The matrix entries must be over a principal ideal domain" raise ValueError(msg) if 0 in m.shape: return () rows, cols = shape = m.shape m = list(m.to_dense().rep.to_ddm()) def add_rows(m, i, j, a, b, c, d): # replace m[i, :] by a*m[i, :] + b*m[j, :] # and m[j, :] by c*m[i, :] + d*m[j, :] for k in range(cols): e = m[i][k] m[i][k] = a*e + b*m[j][k] m[j][k] = c*e + d*m[j][k] def clear_column(m): # make m[1:, 0] zero by row and column operations if m[0][0] == 0: return m # pragma: nocover pivot = m[0][0] for j in range(1, rows): if m[j][0] == 0: continue d, r = domain.div(m[j][0], pivot) if r == 0: add_rows(m, 0, j, 1, 0, -d, 1) else: a, b, g = domain.gcdex(pivot, m[j][0]) d_0 = domain.div(m[j][0], g)[0] d_j = domain.div(pivot, g)[0] add_rows(m, 0, j, a, b, d_0, -d_j) pivot = g return m def clear_row(m): # make m[0, 1:] zero by row and column operations if m[0][0] == 0: return m # pragma: nocover pivot = m[0][0] for j in range(1, cols): if m[0][j] == 0: continue d, r = domain.div(m[0][j], pivot) if r == 0: add_columns(m, 0, j, 1, 0, -d, 1) else: a, b, g = domain.gcdex(pivot, m[0][j]) d_0 = domain.div(m[0][j], g)[0] d_j = domain.div(pivot, g)[0] add_columns(m, 0, j, a, b, d_0, -d_j) pivot = g return m # permute the rows and columns until m[0,0] is non-zero if possible ind = [i for i in range(rows) if m[i][0] != 0] if ind and ind[0] != 0: m[0], m[ind[0]] = m[ind[0]], m[0] else: ind = [j for j in range(cols) if m[0][j] != 0] if ind and ind[0] != 0: for row in m: row[0], row[ind[0]] = row[ind[0]], row[0] # make the first row and column except m[0,0] zero while (any(m[0][i] != 0 for i in range(1,cols)) or any(m[i][0] != 0 for i in range(1,rows))): m = clear_column(m) m = clear_row(m) if 1 in shape: invs = () else: lower_right = DomainMatrix([r[1:] for r in m[1:]], (rows-1, cols-1), domain) invs = invariant_factors(lower_right) if m[0][0]: result = [m[0][0]] result.extend(invs) # in case m[0] doesn't divide the invariants of the rest of the matrix for i in range(len(result)-1): if result[i] and domain.div(result[i+1], result[i])[1] != 0: g = domain.gcd(result[i+1], result[i]) result[i+1] = domain.div(result[i], g)[0]*result[i+1] result[i] = g else: break else: result = invs + (m[0][0],) return tuple(result) def _gcdex(a, b): r""" This supports the functions that compute Hermite Normal Form. Explanation =========== Let x, y be the coefficients returned by the extended Euclidean Algorithm, so that x*a + y*b = g. In the algorithms for computing HNF, it is critical that x, y not only satisfy the condition of being small in magnitude -- namely that |x| <= |b|/g, |y| <- |a|/g -- but also that y == 0 when a | b. """ x, y, g = ZZ.gcdex(a, b) if a != 0 and b % a == 0: y = 0 x = -1 if a < 0 else 1 return x, y, g def _hermite_normal_form(A): r""" Compute the Hermite Normal Form of DomainMatrix *A* over :ref:`ZZ`. Parameters ========== A : :py:class:`~.DomainMatrix` over domain :ref:`ZZ`. Returns ======= :py:class:`~.DomainMatrix` The HNF of matrix *A*. Raises ====== DMDomainError If the domain of the matrix is not :ref:`ZZ`. References ========== .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.* (See Algorithm 2.4.5.) """ if not A.domain.is_ZZ: raise DMDomainError('Matrix must be over domain ZZ.') # We work one row at a time, starting from the bottom row, and working our # way up. m, n = A.shape A = A.to_dense().rep.to_ddm().copy() # Our goal is to put pivot entries in the rightmost columns. # Invariant: Before processing each row, k should be the index of the # leftmost column in which we have so far put a pivot. k = n for i in range(m - 1, -1, -1): if k == 0: # This case can arise when n < m and we've already found n pivots. # We don't need to consider any more rows, because this is already # the maximum possible number of pivots. break k -= 1 # k now points to the column in which we want to put a pivot. # We want zeros in all entries to the left of the pivot column. for j in range(k - 1, -1, -1): if A[i][j] != 0: # Replace cols j, k by lin combs of these cols such that, in row i, # col j has 0, while col k has the gcd of their row i entries. Note # that this ensures a nonzero entry in col k. u, v, d = _gcdex(A[i][k], A[i][j]) r, s = A[i][k] // d, A[i][j] // d add_columns(A, k, j, u, v, -s, r) b = A[i][k] # Do not want the pivot entry to be negative. if b < 0: add_columns(A, k, k, -1, 0, -1, 0) b = -b # The pivot entry will be 0 iff the row was 0 from the pivot col all the # way to the left. In this case, we are still working on the same pivot # col for the next row. Therefore: if b == 0: k += 1 # If the pivot entry is nonzero, then we want to reduce all entries to its # right in the sense of the division algorithm, i.e. make them all remainders # w.r.t. the pivot as divisor. else: for j in range(k + 1, n): q = A[i][j] // b add_columns(A, j, k, 1, -q, 0, 1) # Finally, the HNF consists of those columns of A in which we succeeded in making # a nonzero pivot. return DomainMatrix.from_rep(A.to_dfm_or_ddm())[:, k:] def _hermite_normal_form_modulo_D(A, D): r""" Perform the mod *D* Hermite Normal Form reduction algorithm on :py:class:`~.DomainMatrix` *A*. Explanation =========== If *A* is an $m \times n$ matrix of rank $m$, having Hermite Normal Form $W$, and if *D* is any positive integer known in advance to be a multiple of $\det(W)$, then the HNF of *A* can be computed by an algorithm that works mod *D* in order to prevent coefficient explosion. Parameters ========== A : :py:class:`~.DomainMatrix` over :ref:`ZZ` $m \times n$ matrix, having rank $m$. D : :ref:`ZZ` Positive integer, known to be a multiple of the determinant of the HNF of *A*. Returns ======= :py:class:`~.DomainMatrix` The HNF of matrix *A*. Raises ====== DMDomainError If the domain of the matrix is not :ref:`ZZ`, or if *D* is given but is not in :ref:`ZZ`. DMShapeError If the matrix has more rows than columns. References ========== .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.* (See Algorithm 2.4.8.) """ if not A.domain.is_ZZ: raise DMDomainError('Matrix must be over domain ZZ.') if not ZZ.of_type(D) or D < 1: raise DMDomainError('Modulus D must be positive element of domain ZZ.') def add_columns_mod_R(m, R, i, j, a, b, c, d): # replace m[:, i] by (a*m[:, i] + b*m[:, j]) % R # and m[:, j] by (c*m[:, i] + d*m[:, j]) % R for k in range(len(m)): e = m[k][i] m[k][i] = symmetric_residue((a * e + b * m[k][j]) % R, R) m[k][j] = symmetric_residue((c * e + d * m[k][j]) % R, R) W = defaultdict(dict) m, n = A.shape if n < m: raise DMShapeError('Matrix must have at least as many columns as rows.') A = A.to_dense().rep.to_ddm().copy() k = n R = D for i in range(m - 1, -1, -1): k -= 1 for j in range(k - 1, -1, -1): if A[i][j] != 0: u, v, d = _gcdex(A[i][k], A[i][j]) r, s = A[i][k] // d, A[i][j] // d add_columns_mod_R(A, R, k, j, u, v, -s, r) b = A[i][k] if b == 0: A[i][k] = b = R u, v, d = _gcdex(b, R) for ii in range(m): W[ii][i] = u*A[ii][k] % R if W[i][i] == 0: W[i][i] = R for j in range(i + 1, m): q = W[i][j] // W[i][i] add_columns(W, j, i, 1, -q, 0, 1) R //= d return DomainMatrix(W, (m, m), ZZ).to_dense() def hermite_normal_form(A, *, D=None, check_rank=False): r""" Compute the Hermite Normal Form of :py:class:`~.DomainMatrix` *A* over :ref:`ZZ`. Examples ======== >>> from sympy import ZZ >>> from sympy.polys.matrices import DomainMatrix >>> from sympy.polys.matrices.normalforms import hermite_normal_form >>> m = DomainMatrix([[ZZ(12), ZZ(6), ZZ(4)], ... [ZZ(3), ZZ(9), ZZ(6)], ... [ZZ(2), ZZ(16), ZZ(14)]], (3, 3), ZZ) >>> print(hermite_normal_form(m).to_Matrix()) Matrix([[10, 0, 2], [0, 15, 3], [0, 0, 2]]) Parameters ========== A : $m \times n$ ``DomainMatrix`` over :ref:`ZZ`. D : :ref:`ZZ`, optional Let $W$ be the HNF of *A*. If known in advance, a positive integer *D* being any multiple of $\det(W)$ may be provided. In this case, if *A* also has rank $m$, then we may use an alternative algorithm that works mod *D* in order to prevent coefficient explosion. check_rank : boolean, optional (default=False) The basic assumption is that, if you pass a value for *D*, then you already believe that *A* has rank $m$, so we do not waste time checking it for you. If you do want this to be checked (and the ordinary, non-modulo *D* algorithm to be used if the check fails), then set *check_rank* to ``True``. Returns ======= :py:class:`~.DomainMatrix` The HNF of matrix *A*. Raises ====== DMDomainError If the domain of the matrix is not :ref:`ZZ`, or if *D* is given but is not in :ref:`ZZ`. DMShapeError If the mod *D* algorithm is used but the matrix has more rows than columns. References ========== .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.* (See Algorithms 2.4.5 and 2.4.8.) """ if not A.domain.is_ZZ: raise DMDomainError('Matrix must be over domain ZZ.') if D is not None and (not check_rank or A.convert_to(QQ).rank() == A.shape[0]): return _hermite_normal_form_modulo_D(A, D) else: return _hermite_normal_form(A)