Coverage for /usr/lib/python3/dist-packages/sympy/polys/matrices/normalforms.py: 8%
157 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1'''Functions returning normal forms of matrices'''
3from collections import defaultdict
5from .domainmatrix import DomainMatrix
6from .exceptions import DMDomainError, DMShapeError
7from sympy.ntheory.modular import symmetric_residue
8from sympy.polys.domains import QQ, ZZ
11# TODO (future work):
12# There are faster algorithms for Smith and Hermite normal forms, which
13# we should implement. See e.g. the Kannan-Bachem algorithm:
14# <https://www.researchgate.net/publication/220617516_Polynomial_Algorithms_for_Computing_the_Smith_and_Hermite_Normal_Forms_of_an_Integer_Matrix>
17def smith_normal_form(m):
18 '''
19 Return the Smith Normal Form of a matrix `m` over the ring `domain`.
20 This will only work if the ring is a principal ideal domain.
22 Examples
23 ========
25 >>> from sympy import ZZ
26 >>> from sympy.polys.matrices import DomainMatrix
27 >>> from sympy.polys.matrices.normalforms import smith_normal_form
28 >>> m = DomainMatrix([[ZZ(12), ZZ(6), ZZ(4)],
29 ... [ZZ(3), ZZ(9), ZZ(6)],
30 ... [ZZ(2), ZZ(16), ZZ(14)]], (3, 3), ZZ)
31 >>> print(smith_normal_form(m).to_Matrix())
32 Matrix([[1, 0, 0], [0, 10, 0], [0, 0, -30]])
34 '''
35 invs = invariant_factors(m)
36 smf = DomainMatrix.diag(invs, m.domain, m.shape)
37 return smf
40def add_columns(m, i, j, a, b, c, d):
41 # replace m[:, i] by a*m[:, i] + b*m[:, j]
42 # and m[:, j] by c*m[:, i] + d*m[:, j]
43 for k in range(len(m)):
44 e = m[k][i]
45 m[k][i] = a*e + b*m[k][j]
46 m[k][j] = c*e + d*m[k][j]
49def invariant_factors(m):
50 '''
51 Return the tuple of abelian invariants for a matrix `m`
52 (as in the Smith-Normal form)
54 References
55 ==========
57 [1] https://en.wikipedia.org/wiki/Smith_normal_form#Algorithm
58 [2] https://web.archive.org/web/20200331143852/https://sierra.nmsu.edu/morandi/notes/SmithNormalForm.pdf
60 '''
61 domain = m.domain
62 if not domain.is_PID:
63 msg = "The matrix entries must be over a principal ideal domain"
64 raise ValueError(msg)
66 if 0 in m.shape:
67 return ()
69 rows, cols = shape = m.shape
70 m = list(m.to_dense().rep)
72 def add_rows(m, i, j, a, b, c, d):
73 # replace m[i, :] by a*m[i, :] + b*m[j, :]
74 # and m[j, :] by c*m[i, :] + d*m[j, :]
75 for k in range(cols):
76 e = m[i][k]
77 m[i][k] = a*e + b*m[j][k]
78 m[j][k] = c*e + d*m[j][k]
80 def clear_column(m):
81 # make m[1:, 0] zero by row and column operations
82 if m[0][0] == 0:
83 return m # pragma: nocover
84 pivot = m[0][0]
85 for j in range(1, rows):
86 if m[j][0] == 0:
87 continue
88 d, r = domain.div(m[j][0], pivot)
89 if r == 0:
90 add_rows(m, 0, j, 1, 0, -d, 1)
91 else:
92 a, b, g = domain.gcdex(pivot, m[j][0])
93 d_0 = domain.div(m[j][0], g)[0]
94 d_j = domain.div(pivot, g)[0]
95 add_rows(m, 0, j, a, b, d_0, -d_j)
96 pivot = g
97 return m
99 def clear_row(m):
100 # make m[0, 1:] zero by row and column operations
101 if m[0][0] == 0:
102 return m # pragma: nocover
103 pivot = m[0][0]
104 for j in range(1, cols):
105 if m[0][j] == 0:
106 continue
107 d, r = domain.div(m[0][j], pivot)
108 if r == 0:
109 add_columns(m, 0, j, 1, 0, -d, 1)
110 else:
111 a, b, g = domain.gcdex(pivot, m[0][j])
112 d_0 = domain.div(m[0][j], g)[0]
113 d_j = domain.div(pivot, g)[0]
114 add_columns(m, 0, j, a, b, d_0, -d_j)
115 pivot = g
116 return m
118 # permute the rows and columns until m[0,0] is non-zero if possible
119 ind = [i for i in range(rows) if m[i][0] != 0]
120 if ind and ind[0] != 0:
121 m[0], m[ind[0]] = m[ind[0]], m[0]
122 else:
123 ind = [j for j in range(cols) if m[0][j] != 0]
124 if ind and ind[0] != 0:
125 for row in m:
126 row[0], row[ind[0]] = row[ind[0]], row[0]
128 # make the first row and column except m[0,0] zero
129 while (any(m[0][i] != 0 for i in range(1,cols)) or
130 any(m[i][0] != 0 for i in range(1,rows))):
131 m = clear_column(m)
132 m = clear_row(m)
134 if 1 in shape:
135 invs = ()
136 else:
137 lower_right = DomainMatrix([r[1:] for r in m[1:]], (rows-1, cols-1), domain)
138 invs = invariant_factors(lower_right)
140 if m[0][0]:
141 result = [m[0][0]]
142 result.extend(invs)
143 # in case m[0] doesn't divide the invariants of the rest of the matrix
144 for i in range(len(result)-1):
145 if result[i] and domain.div(result[i+1], result[i])[1] != 0:
146 g = domain.gcd(result[i+1], result[i])
147 result[i+1] = domain.div(result[i], g)[0]*result[i+1]
148 result[i] = g
149 else:
150 break
151 else:
152 result = invs + (m[0][0],)
153 return tuple(result)
156def _gcdex(a, b):
157 r"""
158 This supports the functions that compute Hermite Normal Form.
160 Explanation
161 ===========
163 Let x, y be the coefficients returned by the extended Euclidean
164 Algorithm, so that x*a + y*b = g. In the algorithms for computing HNF,
165 it is critical that x, y not only satisfy the condition of being small
166 in magnitude -- namely that |x| <= |b|/g, |y| <- |a|/g -- but also that
167 y == 0 when a | b.
169 """
170 x, y, g = ZZ.gcdex(a, b)
171 if a != 0 and b % a == 0:
172 y = 0
173 x = -1 if a < 0 else 1
174 return x, y, g
177def _hermite_normal_form(A):
178 r"""
179 Compute the Hermite Normal Form of DomainMatrix *A* over :ref:`ZZ`.
181 Parameters
182 ==========
184 A : :py:class:`~.DomainMatrix` over domain :ref:`ZZ`.
186 Returns
187 =======
189 :py:class:`~.DomainMatrix`
190 The HNF of matrix *A*.
192 Raises
193 ======
195 DMDomainError
196 If the domain of the matrix is not :ref:`ZZ`.
198 References
199 ==========
201 .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
202 (See Algorithm 2.4.5.)
204 """
205 if not A.domain.is_ZZ:
206 raise DMDomainError('Matrix must be over domain ZZ.')
207 # We work one row at a time, starting from the bottom row, and working our
208 # way up.
209 m, n = A.shape
210 A = A.to_dense().rep.copy()
211 # Our goal is to put pivot entries in the rightmost columns.
212 # Invariant: Before processing each row, k should be the index of the
213 # leftmost column in which we have so far put a pivot.
214 k = n
215 for i in range(m - 1, -1, -1):
216 if k == 0:
217 # This case can arise when n < m and we've already found n pivots.
218 # We don't need to consider any more rows, because this is already
219 # the maximum possible number of pivots.
220 break
221 k -= 1
222 # k now points to the column in which we want to put a pivot.
223 # We want zeros in all entries to the left of the pivot column.
224 for j in range(k - 1, -1, -1):
225 if A[i][j] != 0:
226 # Replace cols j, k by lin combs of these cols such that, in row i,
227 # col j has 0, while col k has the gcd of their row i entries. Note
228 # that this ensures a nonzero entry in col k.
229 u, v, d = _gcdex(A[i][k], A[i][j])
230 r, s = A[i][k] // d, A[i][j] // d
231 add_columns(A, k, j, u, v, -s, r)
232 b = A[i][k]
233 # Do not want the pivot entry to be negative.
234 if b < 0:
235 add_columns(A, k, k, -1, 0, -1, 0)
236 b = -b
237 # The pivot entry will be 0 iff the row was 0 from the pivot col all the
238 # way to the left. In this case, we are still working on the same pivot
239 # col for the next row. Therefore:
240 if b == 0:
241 k += 1
242 # If the pivot entry is nonzero, then we want to reduce all entries to its
243 # right in the sense of the division algorithm, i.e. make them all remainders
244 # w.r.t. the pivot as divisor.
245 else:
246 for j in range(k + 1, n):
247 q = A[i][j] // b
248 add_columns(A, j, k, 1, -q, 0, 1)
249 # Finally, the HNF consists of those columns of A in which we succeeded in making
250 # a nonzero pivot.
251 return DomainMatrix.from_rep(A)[:, k:]
254def _hermite_normal_form_modulo_D(A, D):
255 r"""
256 Perform the mod *D* Hermite Normal Form reduction algorithm on
257 :py:class:`~.DomainMatrix` *A*.
259 Explanation
260 ===========
262 If *A* is an $m \times n$ matrix of rank $m$, having Hermite Normal Form
263 $W$, and if *D* is any positive integer known in advance to be a multiple
264 of $\det(W)$, then the HNF of *A* can be computed by an algorithm that
265 works mod *D* in order to prevent coefficient explosion.
267 Parameters
268 ==========
270 A : :py:class:`~.DomainMatrix` over :ref:`ZZ`
271 $m \times n$ matrix, having rank $m$.
272 D : :ref:`ZZ`
273 Positive integer, known to be a multiple of the determinant of the
274 HNF of *A*.
276 Returns
277 =======
279 :py:class:`~.DomainMatrix`
280 The HNF of matrix *A*.
282 Raises
283 ======
285 DMDomainError
286 If the domain of the matrix is not :ref:`ZZ`, or
287 if *D* is given but is not in :ref:`ZZ`.
289 DMShapeError
290 If the matrix has more rows than columns.
292 References
293 ==========
295 .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
296 (See Algorithm 2.4.8.)
298 """
299 if not A.domain.is_ZZ:
300 raise DMDomainError('Matrix must be over domain ZZ.')
301 if not ZZ.of_type(D) or D < 1:
302 raise DMDomainError('Modulus D must be positive element of domain ZZ.')
304 def add_columns_mod_R(m, R, i, j, a, b, c, d):
305 # replace m[:, i] by (a*m[:, i] + b*m[:, j]) % R
306 # and m[:, j] by (c*m[:, i] + d*m[:, j]) % R
307 for k in range(len(m)):
308 e = m[k][i]
309 m[k][i] = symmetric_residue((a * e + b * m[k][j]) % R, R)
310 m[k][j] = symmetric_residue((c * e + d * m[k][j]) % R, R)
312 W = defaultdict(dict)
314 m, n = A.shape
315 if n < m:
316 raise DMShapeError('Matrix must have at least as many columns as rows.')
317 A = A.to_dense().rep.copy()
318 k = n
319 R = D
320 for i in range(m - 1, -1, -1):
321 k -= 1
322 for j in range(k - 1, -1, -1):
323 if A[i][j] != 0:
324 u, v, d = _gcdex(A[i][k], A[i][j])
325 r, s = A[i][k] // d, A[i][j] // d
326 add_columns_mod_R(A, R, k, j, u, v, -s, r)
327 b = A[i][k]
328 if b == 0:
329 A[i][k] = b = R
330 u, v, d = _gcdex(b, R)
331 for ii in range(m):
332 W[ii][i] = u*A[ii][k] % R
333 if W[i][i] == 0:
334 W[i][i] = R
335 for j in range(i + 1, m):
336 q = W[i][j] // W[i][i]
337 add_columns(W, j, i, 1, -q, 0, 1)
338 R //= d
339 return DomainMatrix(W, (m, m), ZZ).to_dense()
342def hermite_normal_form(A, *, D=None, check_rank=False):
343 r"""
344 Compute the Hermite Normal Form of :py:class:`~.DomainMatrix` *A* over
345 :ref:`ZZ`.
347 Examples
348 ========
350 >>> from sympy import ZZ
351 >>> from sympy.polys.matrices import DomainMatrix
352 >>> from sympy.polys.matrices.normalforms import hermite_normal_form
353 >>> m = DomainMatrix([[ZZ(12), ZZ(6), ZZ(4)],
354 ... [ZZ(3), ZZ(9), ZZ(6)],
355 ... [ZZ(2), ZZ(16), ZZ(14)]], (3, 3), ZZ)
356 >>> print(hermite_normal_form(m).to_Matrix())
357 Matrix([[10, 0, 2], [0, 15, 3], [0, 0, 2]])
359 Parameters
360 ==========
362 A : $m \times n$ ``DomainMatrix`` over :ref:`ZZ`.
364 D : :ref:`ZZ`, optional
365 Let $W$ be the HNF of *A*. If known in advance, a positive integer *D*
366 being any multiple of $\det(W)$ may be provided. In this case, if *A*
367 also has rank $m$, then we may use an alternative algorithm that works
368 mod *D* in order to prevent coefficient explosion.
370 check_rank : boolean, optional (default=False)
371 The basic assumption is that, if you pass a value for *D*, then
372 you already believe that *A* has rank $m$, so we do not waste time
373 checking it for you. If you do want this to be checked (and the
374 ordinary, non-modulo *D* algorithm to be used if the check fails), then
375 set *check_rank* to ``True``.
377 Returns
378 =======
380 :py:class:`~.DomainMatrix`
381 The HNF of matrix *A*.
383 Raises
384 ======
386 DMDomainError
387 If the domain of the matrix is not :ref:`ZZ`, or
388 if *D* is given but is not in :ref:`ZZ`.
390 DMShapeError
391 If the mod *D* algorithm is used but the matrix has more rows than
392 columns.
394 References
395 ==========
397 .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
398 (See Algorithms 2.4.5 and 2.4.8.)
400 """
401 if not A.domain.is_ZZ:
402 raise DMDomainError('Matrix must be over domain ZZ.')
403 if D is not None and (not check_rank or A.convert_to(QQ).rank() == A.shape[0]):
404 return _hermite_normal_form_modulo_D(A, D)
405 else:
406 return _hermite_normal_form(A)