Coverage for /usr/lib/python3/dist-packages/sympy/polys/matrices/dense.py: 10%
204 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
3Module for the ddm_* routines for operating on a matrix in list of lists
4matrix representation.
6These routines are used internally by the DDM class which also provides a
7friendlier interface for them. The idea here is to implement core matrix
8routines in a way that can be applied to any simple list representation
9without the need to use any particular matrix class. For example we can
10compute the RREF of a matrix like:
12 >>> from sympy.polys.matrices.dense import ddm_irref
13 >>> M = [[1, 2, 3], [4, 5, 6]]
14 >>> pivots = ddm_irref(M)
15 >>> M
16 [[1.0, 0.0, -1.0], [0, 1.0, 2.0]]
18These are lower-level routines that work mostly in place.The routines at this
19level should not need to know what the domain of the elements is but should
20ideally document what operations they will use and what functions they need to
21be provided with.
23The next-level up is the DDM class which uses these routines but wraps them up
24with an interface that handles copying etc and keeps track of the Domain of
25the elements of the matrix:
27 >>> from sympy.polys.domains import QQ
28 >>> from sympy.polys.matrices.ddm import DDM
29 >>> M = DDM([[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)]], (2, 3), QQ)
30 >>> M
31 [[1, 2, 3], [4, 5, 6]]
32 >>> Mrref, pivots = M.rref()
33 >>> Mrref
34 [[1, 0, -1], [0, 1, 2]]
36"""
37from __future__ import annotations
38from operator import mul
39from .exceptions import (
40 DMShapeError,
41 DMNonInvertibleMatrixError,
42 DMNonSquareMatrixError,
43)
44from typing import Sequence, TypeVar
45from sympy.polys.matrices._typing import RingElement
48T = TypeVar('T')
49R = TypeVar('R', bound=RingElement)
52def ddm_transpose(matrix: Sequence[Sequence[T]]) -> list[list[T]]:
53 """matrix transpose"""
54 return list(map(list, zip(*matrix)))
57def ddm_iadd(a: list[list[R]], b: Sequence[Sequence[R]]) -> None:
58 """a += b"""
59 for ai, bi in zip(a, b):
60 for j, bij in enumerate(bi):
61 ai[j] += bij
64def ddm_isub(a: list[list[R]], b: Sequence[Sequence[R]]) -> None:
65 """a -= b"""
66 for ai, bi in zip(a, b):
67 for j, bij in enumerate(bi):
68 ai[j] -= bij
71def ddm_ineg(a: list[list[R]]) -> None:
72 """a <-- -a"""
73 for ai in a:
74 for j, aij in enumerate(ai):
75 ai[j] = -aij
78def ddm_imul(a: list[list[R]], b: R) -> None:
79 for ai in a:
80 for j, aij in enumerate(ai):
81 ai[j] = aij * b
84def ddm_irmul(a: list[list[R]], b: R) -> None:
85 for ai in a:
86 for j, aij in enumerate(ai):
87 ai[j] = b * aij
90def ddm_imatmul(
91 a: list[list[R]], b: Sequence[Sequence[R]], c: Sequence[Sequence[R]]
92) -> None:
93 """a += b @ c"""
94 cT = list(zip(*c))
96 for bi, ai in zip(b, a):
97 for j, cTj in enumerate(cT):
98 ai[j] = sum(map(mul, bi, cTj), ai[j])
101def ddm_irref(a, _partial_pivot=False):
102 """a <-- rref(a)"""
103 # a is (m x n)
104 m = len(a)
105 if not m:
106 return []
107 n = len(a[0])
109 i = 0
110 pivots = []
112 for j in range(n):
113 # Proper pivoting should be used for all domains for performance
114 # reasons but it is only strictly needed for RR and CC (and possibly
115 # other domains like RR(x)). This path is used by DDM.rref() if the
116 # domain is RR or CC. It uses partial (row) pivoting based on the
117 # absolute value of the pivot candidates.
118 if _partial_pivot:
119 ip = max(range(i, m), key=lambda ip: abs(a[ip][j]))
120 a[i], a[ip] = a[ip], a[i]
122 # pivot
123 aij = a[i][j]
125 # zero-pivot
126 if not aij:
127 for ip in range(i+1, m):
128 aij = a[ip][j]
129 # row-swap
130 if aij:
131 a[i], a[ip] = a[ip], a[i]
132 break
133 else:
134 # next column
135 continue
137 # normalise row
138 ai = a[i]
139 aijinv = aij**-1
140 for l in range(j, n):
141 ai[l] *= aijinv # ai[j] = one
143 # eliminate above and below to the right
144 for k, ak in enumerate(a):
145 if k == i or not ak[j]:
146 continue
147 akj = ak[j]
148 ak[j] -= akj # ak[j] = zero
149 for l in range(j+1, n):
150 ak[l] -= akj * ai[l]
152 # next row
153 pivots.append(j)
154 i += 1
156 # no more rows?
157 if i >= m:
158 break
160 return pivots
163def ddm_idet(a, K):
164 """a <-- echelon(a); return det"""
165 # Bareiss algorithm
166 # https://www.math.usm.edu/perry/Research/Thesis_DRL.pdf
168 # a is (m x n)
169 m = len(a)
170 if not m:
171 return K.one
172 n = len(a[0])
174 exquo = K.exquo
175 # uf keeps track of the sign change from row swaps
176 uf = K.one
178 for k in range(n-1):
179 if not a[k][k]:
180 for i in range(k+1, n):
181 if a[i][k]:
182 a[k], a[i] = a[i], a[k]
183 uf = -uf
184 break
185 else:
186 return K.zero
188 akkm1 = a[k-1][k-1] if k else K.one
190 for i in range(k+1, n):
191 for j in range(k+1, n):
192 a[i][j] = exquo(a[i][j]*a[k][k] - a[i][k]*a[k][j], akkm1)
194 return uf * a[-1][-1]
197def ddm_iinv(ainv, a, K):
198 if not K.is_Field:
199 raise ValueError('Not a field')
201 # a is (m x n)
202 m = len(a)
203 if not m:
204 return
205 n = len(a[0])
206 if m != n:
207 raise DMNonSquareMatrixError
209 eye = [[K.one if i==j else K.zero for j in range(n)] for i in range(n)]
210 Aaug = [row + eyerow for row, eyerow in zip(a, eye)]
211 pivots = ddm_irref(Aaug)
212 if pivots != list(range(n)):
213 raise DMNonInvertibleMatrixError('Matrix det == 0; not invertible.')
214 ainv[:] = [row[n:] for row in Aaug]
217def ddm_ilu_split(L, U, K):
218 """L, U <-- LU(U)"""
219 m = len(U)
220 if not m:
221 return []
222 n = len(U[0])
224 swaps = ddm_ilu(U)
226 zeros = [K.zero] * min(m, n)
227 for i in range(1, m):
228 j = min(i, n)
229 L[i][:j] = U[i][:j]
230 U[i][:j] = zeros[:j]
232 return swaps
235def ddm_ilu(a):
236 """a <-- LU(a)"""
237 m = len(a)
238 if not m:
239 return []
240 n = len(a[0])
242 swaps = []
244 for i in range(min(m, n)):
245 if not a[i][i]:
246 for ip in range(i+1, m):
247 if a[ip][i]:
248 swaps.append((i, ip))
249 a[i], a[ip] = a[ip], a[i]
250 break
251 else:
252 # M = Matrix([[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 1, 2]])
253 continue
254 for j in range(i+1, m):
255 l_ji = a[j][i] / a[i][i]
256 a[j][i] = l_ji
257 for k in range(i+1, n):
258 a[j][k] -= l_ji * a[i][k]
260 return swaps
263def ddm_ilu_solve(x, L, U, swaps, b):
264 """x <-- solve(L*U*x = swaps(b))"""
265 m = len(U)
266 if not m:
267 return
268 n = len(U[0])
270 m2 = len(b)
271 if not m2:
272 raise DMShapeError("Shape mismtch")
273 o = len(b[0])
275 if m != m2:
276 raise DMShapeError("Shape mismtch")
277 if m < n:
278 raise NotImplementedError("Underdetermined")
280 if swaps:
281 b = [row[:] for row in b]
282 for i1, i2 in swaps:
283 b[i1], b[i2] = b[i2], b[i1]
285 # solve Ly = b
286 y = [[None] * o for _ in range(m)]
287 for k in range(o):
288 for i in range(m):
289 rhs = b[i][k]
290 for j in range(i):
291 rhs -= L[i][j] * y[j][k]
292 y[i][k] = rhs
294 if m > n:
295 for i in range(n, m):
296 for j in range(o):
297 if y[i][j]:
298 raise DMNonInvertibleMatrixError
300 # Solve Ux = y
301 for k in range(o):
302 for i in reversed(range(n)):
303 if not U[i][i]:
304 raise DMNonInvertibleMatrixError
305 rhs = y[i][k]
306 for j in range(i+1, n):
307 rhs -= U[i][j] * x[j][k]
308 x[i][k] = rhs / U[i][i]
311def ddm_berk(M, K):
312 m = len(M)
313 if not m:
314 return [[K.one]]
315 n = len(M[0])
317 if m != n:
318 raise DMShapeError("Not square")
320 if n == 1:
321 return [[K.one], [-M[0][0]]]
323 a = M[0][0]
324 R = [M[0][1:]]
325 C = [[row[0]] for row in M[1:]]
326 A = [row[1:] for row in M[1:]]
328 q = ddm_berk(A, K)
330 T = [[K.zero] * n for _ in range(n+1)]
331 for i in range(n):
332 T[i][i] = K.one
333 T[i+1][i] = -a
334 for i in range(2, n+1):
335 if i == 2:
336 AnC = C
337 else:
338 C = AnC
339 AnC = [[K.zero] for row in C]
340 ddm_imatmul(AnC, A, C)
341 RAnC = [[K.zero]]
342 ddm_imatmul(RAnC, R, AnC)
343 for j in range(0, n+1-i):
344 T[i+j][j] = -RAnC[0][0]
346 qout = [[K.zero] for _ in range(n+1)]
347 ddm_imatmul(qout, T, q)
348 return qout