Coverage for /usr/lib/python3/dist-packages/sympy/matrices/reductions.py: 9%
99 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from types import FunctionType
3from .utilities import _get_intermediate_simp, _iszero, _dotprodsimp, _simplify
4from .determinant import _find_reasonable_pivot
7def _row_reduce_list(mat, rows, cols, one, iszerofunc, simpfunc,
8 normalize_last=True, normalize=True, zero_above=True):
9 """Row reduce a flat list representation of a matrix and return a tuple
10 (rref_matrix, pivot_cols, swaps) where ``rref_matrix`` is a flat list,
11 ``pivot_cols`` are the pivot columns and ``swaps`` are any row swaps that
12 were used in the process of row reduction.
14 Parameters
15 ==========
17 mat : list
18 list of matrix elements, must be ``rows`` * ``cols`` in length
20 rows, cols : integer
21 number of rows and columns in flat list representation
23 one : SymPy object
24 represents the value one, from ``Matrix.one``
26 iszerofunc : determines if an entry can be used as a pivot
28 simpfunc : used to simplify elements and test if they are
29 zero if ``iszerofunc`` returns `None`
31 normalize_last : indicates where all row reduction should
32 happen in a fraction-free manner and then the rows are
33 normalized (so that the pivots are 1), or whether
34 rows should be normalized along the way (like the naive
35 row reduction algorithm)
37 normalize : whether pivot rows should be normalized so that
38 the pivot value is 1
40 zero_above : whether entries above the pivot should be zeroed.
41 If ``zero_above=False``, an echelon matrix will be returned.
42 """
44 def get_col(i):
45 return mat[i::cols]
47 def row_swap(i, j):
48 mat[i*cols:(i + 1)*cols], mat[j*cols:(j + 1)*cols] = \
49 mat[j*cols:(j + 1)*cols], mat[i*cols:(i + 1)*cols]
51 def cross_cancel(a, i, b, j):
52 """Does the row op row[i] = a*row[i] - b*row[j]"""
53 q = (j - i)*cols
54 for p in range(i*cols, (i + 1)*cols):
55 mat[p] = isimp(a*mat[p] - b*mat[p + q])
57 isimp = _get_intermediate_simp(_dotprodsimp)
58 piv_row, piv_col = 0, 0
59 pivot_cols = []
60 swaps = []
62 # use a fraction free method to zero above and below each pivot
63 while piv_col < cols and piv_row < rows:
64 pivot_offset, pivot_val, \
65 assumed_nonzero, newly_determined = _find_reasonable_pivot(
66 get_col(piv_col)[piv_row:], iszerofunc, simpfunc)
68 # _find_reasonable_pivot may have simplified some things
69 # in the process. Let's not let them go to waste
70 for (offset, val) in newly_determined:
71 offset += piv_row
72 mat[offset*cols + piv_col] = val
74 if pivot_offset is None:
75 piv_col += 1
76 continue
78 pivot_cols.append(piv_col)
79 if pivot_offset != 0:
80 row_swap(piv_row, pivot_offset + piv_row)
81 swaps.append((piv_row, pivot_offset + piv_row))
83 # if we aren't normalizing last, we normalize
84 # before we zero the other rows
85 if normalize_last is False:
86 i, j = piv_row, piv_col
87 mat[i*cols + j] = one
88 for p in range(i*cols + j + 1, (i + 1)*cols):
89 mat[p] = isimp(mat[p] / pivot_val)
90 # after normalizing, the pivot value is 1
91 pivot_val = one
93 # zero above and below the pivot
94 for row in range(rows):
95 # don't zero our current row
96 if row == piv_row:
97 continue
98 # don't zero above the pivot unless we're told.
99 if zero_above is False and row < piv_row:
100 continue
101 # if we're already a zero, don't do anything
102 val = mat[row*cols + piv_col]
103 if iszerofunc(val):
104 continue
106 cross_cancel(pivot_val, row, val, piv_row)
107 piv_row += 1
109 # normalize each row
110 if normalize_last is True and normalize is True:
111 for piv_i, piv_j in enumerate(pivot_cols):
112 pivot_val = mat[piv_i*cols + piv_j]
113 mat[piv_i*cols + piv_j] = one
114 for p in range(piv_i*cols + piv_j + 1, (piv_i + 1)*cols):
115 mat[p] = isimp(mat[p] / pivot_val)
117 return mat, tuple(pivot_cols), tuple(swaps)
120# This functions is a candidate for caching if it gets implemented for matrices.
121def _row_reduce(M, iszerofunc, simpfunc, normalize_last=True,
122 normalize=True, zero_above=True):
124 mat, pivot_cols, swaps = _row_reduce_list(list(M), M.rows, M.cols, M.one,
125 iszerofunc, simpfunc, normalize_last=normalize_last,
126 normalize=normalize, zero_above=zero_above)
128 return M._new(M.rows, M.cols, mat), pivot_cols, swaps
131def _is_echelon(M, iszerofunc=_iszero):
132 """Returns `True` if the matrix is in echelon form. That is, all rows of
133 zeros are at the bottom, and below each leading non-zero in a row are
134 exclusively zeros."""
136 if M.rows <= 0 or M.cols <= 0:
137 return True
139 zeros_below = all(iszerofunc(t) for t in M[1:, 0])
141 if iszerofunc(M[0, 0]):
142 return zeros_below and _is_echelon(M[:, 1:], iszerofunc)
144 return zeros_below and _is_echelon(M[1:, 1:], iszerofunc)
147def _echelon_form(M, iszerofunc=_iszero, simplify=False, with_pivots=False):
148 """Returns a matrix row-equivalent to ``M`` that is in echelon form. Note
149 that echelon form of a matrix is *not* unique, however, properties like the
150 row space and the null space are preserved.
152 Examples
153 ========
155 >>> from sympy import Matrix
156 >>> M = Matrix([[1, 2], [3, 4]])
157 >>> M.echelon_form()
158 Matrix([
159 [1, 2],
160 [0, -2]])
161 """
163 simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
165 mat, pivots, _ = _row_reduce(M, iszerofunc, simpfunc,
166 normalize_last=True, normalize=False, zero_above=False)
168 if with_pivots:
169 return mat, pivots
171 return mat
174# This functions is a candidate for caching if it gets implemented for matrices.
175def _rank(M, iszerofunc=_iszero, simplify=False):
176 """Returns the rank of a matrix.
178 Examples
179 ========
181 >>> from sympy import Matrix
182 >>> from sympy.abc import x
183 >>> m = Matrix([[1, 2], [x, 1 - 1/x]])
184 >>> m.rank()
185 2
186 >>> n = Matrix(3, 3, range(1, 10))
187 >>> n.rank()
188 2
189 """
191 def _permute_complexity_right(M, iszerofunc):
192 """Permute columns with complicated elements as
193 far right as they can go. Since the ``sympy`` row reduction
194 algorithms start on the left, having complexity right-shifted
195 speeds things up.
197 Returns a tuple (mat, perm) where perm is a permutation
198 of the columns to perform to shift the complex columns right, and mat
199 is the permuted matrix."""
201 def complexity(i):
202 # the complexity of a column will be judged by how many
203 # element's zero-ness cannot be determined
204 return sum(1 if iszerofunc(e) is None else 0 for e in M[:, i])
206 complex = [(complexity(i), i) for i in range(M.cols)]
207 perm = [j for (i, j) in sorted(complex)]
209 return (M.permute(perm, orientation='cols'), perm)
211 simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
213 # for small matrices, we compute the rank explicitly
214 # if is_zero on elements doesn't answer the question
215 # for small matrices, we fall back to the full routine.
216 if M.rows <= 0 or M.cols <= 0:
217 return 0
219 if M.rows <= 1 or M.cols <= 1:
220 zeros = [iszerofunc(x) for x in M]
222 if False in zeros:
223 return 1
225 if M.rows == 2 and M.cols == 2:
226 zeros = [iszerofunc(x) for x in M]
228 if False not in zeros and None not in zeros:
229 return 0
231 d = M.det()
233 if iszerofunc(d) and False in zeros:
234 return 1
235 if iszerofunc(d) is False:
236 return 2
238 mat, _ = _permute_complexity_right(M, iszerofunc=iszerofunc)
239 _, pivots, _ = _row_reduce(mat, iszerofunc, simpfunc, normalize_last=True,
240 normalize=False, zero_above=False)
242 return len(pivots)
245def _rref(M, iszerofunc=_iszero, simplify=False, pivots=True,
246 normalize_last=True):
247 """Return reduced row-echelon form of matrix and indices of pivot vars.
249 Parameters
250 ==========
252 iszerofunc : Function
253 A function used for detecting whether an element can
254 act as a pivot. ``lambda x: x.is_zero`` is used by default.
256 simplify : Function
257 A function used to simplify elements when looking for a pivot.
258 By default SymPy's ``simplify`` is used.
260 pivots : True or False
261 If ``True``, a tuple containing the row-reduced matrix and a tuple
262 of pivot columns is returned. If ``False`` just the row-reduced
263 matrix is returned.
265 normalize_last : True or False
266 If ``True``, no pivots are normalized to `1` until after all
267 entries above and below each pivot are zeroed. This means the row
268 reduction algorithm is fraction free until the very last step.
269 If ``False``, the naive row reduction procedure is used where
270 each pivot is normalized to be `1` before row operations are
271 used to zero above and below the pivot.
273 Examples
274 ========
276 >>> from sympy import Matrix
277 >>> from sympy.abc import x
278 >>> m = Matrix([[1, 2], [x, 1 - 1/x]])
279 >>> m.rref()
280 (Matrix([
281 [1, 0],
282 [0, 1]]), (0, 1))
283 >>> rref_matrix, rref_pivots = m.rref()
284 >>> rref_matrix
285 Matrix([
286 [1, 0],
287 [0, 1]])
288 >>> rref_pivots
289 (0, 1)
291 ``iszerofunc`` can correct rounding errors in matrices with float
292 values. In the following example, calling ``rref()`` leads to
293 floating point errors, incorrectly row reducing the matrix.
294 ``iszerofunc= lambda x: abs(x)<1e-9`` sets sufficiently small numbers
295 to zero, avoiding this error.
297 >>> m = Matrix([[0.9, -0.1, -0.2, 0], [-0.8, 0.9, -0.4, 0], [-0.1, -0.8, 0.6, 0]])
298 >>> m.rref()
299 (Matrix([
300 [1, 0, 0, 0],
301 [0, 1, 0, 0],
302 [0, 0, 1, 0]]), (0, 1, 2))
303 >>> m.rref(iszerofunc=lambda x:abs(x)<1e-9)
304 (Matrix([
305 [1, 0, -0.301369863013699, 0],
306 [0, 1, -0.712328767123288, 0],
307 [0, 0, 0, 0]]), (0, 1))
309 Notes
310 =====
312 The default value of ``normalize_last=True`` can provide significant
313 speedup to row reduction, especially on matrices with symbols. However,
314 if you depend on the form row reduction algorithm leaves entries
315 of the matrix, set ``noramlize_last=False``
316 """
318 simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
320 mat, pivot_cols, _ = _row_reduce(M, iszerofunc, simpfunc,
321 normalize_last, normalize=True, zero_above=True)
323 if pivots:
324 mat = (mat, pivot_cols)
326 return mat