Coverage for /usr/lib/python3/dist-packages/sympy/matrices/determinant.py: 11%
239 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from types import FunctionType
3from sympy.core.numbers import Float, Integer
4from sympy.core.singleton import S
5from sympy.core.symbol import uniquely_named_symbol
6from sympy.core.mul import Mul
7from sympy.polys import PurePoly, cancel
8from sympy.functions.combinatorial.numbers import nC
9from sympy.polys.matrices.domainmatrix import DomainMatrix
11from .common import NonSquareMatrixError
12from .utilities import (
13 _get_intermediate_simp, _get_intermediate_simp_bool,
14 _iszero, _is_zero_after_expand_mul, _dotprodsimp, _simplify)
17def _find_reasonable_pivot(col, iszerofunc=_iszero, simpfunc=_simplify):
18 """ Find the lowest index of an item in ``col`` that is
19 suitable for a pivot. If ``col`` consists only of
20 Floats, the pivot with the largest norm is returned.
21 Otherwise, the first element where ``iszerofunc`` returns
22 False is used. If ``iszerofunc`` does not return false,
23 items are simplified and retested until a suitable
24 pivot is found.
26 Returns a 4-tuple
27 (pivot_offset, pivot_val, assumed_nonzero, newly_determined)
28 where pivot_offset is the index of the pivot, pivot_val is
29 the (possibly simplified) value of the pivot, assumed_nonzero
30 is True if an assumption that the pivot was non-zero
31 was made without being proved, and newly_determined are
32 elements that were simplified during the process of pivot
33 finding."""
35 newly_determined = []
36 col = list(col)
37 # a column that contains a mix of floats and integers
38 # but at least one float is considered a numerical
39 # column, and so we do partial pivoting
40 if all(isinstance(x, (Float, Integer)) for x in col) and any(
41 isinstance(x, Float) for x in col):
42 col_abs = [abs(x) for x in col]
43 max_value = max(col_abs)
44 if iszerofunc(max_value):
45 # just because iszerofunc returned True, doesn't
46 # mean the value is numerically zero. Make sure
47 # to replace all entries with numerical zeros
48 if max_value != 0:
49 newly_determined = [(i, 0) for i, x in enumerate(col) if x != 0]
50 return (None, None, False, newly_determined)
51 index = col_abs.index(max_value)
52 return (index, col[index], False, newly_determined)
54 # PASS 1 (iszerofunc directly)
55 possible_zeros = []
56 for i, x in enumerate(col):
57 is_zero = iszerofunc(x)
58 # is someone wrote a custom iszerofunc, it may return
59 # BooleanFalse or BooleanTrue instead of True or False,
60 # so use == for comparison instead of `is`
61 if is_zero == False:
62 # we found something that is definitely not zero
63 return (i, x, False, newly_determined)
64 possible_zeros.append(is_zero)
66 # by this point, we've found no certain non-zeros
67 if all(possible_zeros):
68 # if everything is definitely zero, we have
69 # no pivot
70 return (None, None, False, newly_determined)
72 # PASS 2 (iszerofunc after simplify)
73 # we haven't found any for-sure non-zeros, so
74 # go through the elements iszerofunc couldn't
75 # make a determination about and opportunistically
76 # simplify to see if we find something
77 for i, x in enumerate(col):
78 if possible_zeros[i] is not None:
79 continue
80 simped = simpfunc(x)
81 is_zero = iszerofunc(simped)
82 if is_zero in (True, False):
83 newly_determined.append((i, simped))
84 if is_zero == False:
85 return (i, simped, False, newly_determined)
86 possible_zeros[i] = is_zero
88 # after simplifying, some things that were recognized
89 # as zeros might be zeros
90 if all(possible_zeros):
91 # if everything is definitely zero, we have
92 # no pivot
93 return (None, None, False, newly_determined)
95 # PASS 3 (.equals(0))
96 # some expressions fail to simplify to zero, but
97 # ``.equals(0)`` evaluates to True. As a last-ditch
98 # attempt, apply ``.equals`` to these expressions
99 for i, x in enumerate(col):
100 if possible_zeros[i] is not None:
101 continue
102 if x.equals(S.Zero):
103 # ``.iszero`` may return False with
104 # an implicit assumption (e.g., ``x.equals(0)``
105 # when ``x`` is a symbol), so only treat it
106 # as proved when ``.equals(0)`` returns True
107 possible_zeros[i] = True
108 newly_determined.append((i, S.Zero))
110 if all(possible_zeros):
111 return (None, None, False, newly_determined)
113 # at this point there is nothing that could definitely
114 # be a pivot. To maintain compatibility with existing
115 # behavior, we'll assume that an illdetermined thing is
116 # non-zero. We should probably raise a warning in this case
117 i = possible_zeros.index(None)
118 return (i, col[i], True, newly_determined)
121def _find_reasonable_pivot_naive(col, iszerofunc=_iszero, simpfunc=None):
122 """
123 Helper that computes the pivot value and location from a
124 sequence of contiguous matrix column elements. As a side effect
125 of the pivot search, this function may simplify some of the elements
126 of the input column. A list of these simplified entries and their
127 indices are also returned.
128 This function mimics the behavior of _find_reasonable_pivot(),
129 but does less work trying to determine if an indeterminate candidate
130 pivot simplifies to zero. This more naive approach can be much faster,
131 with the trade-off that it may erroneously return a pivot that is zero.
133 ``col`` is a sequence of contiguous column entries to be searched for
134 a suitable pivot.
135 ``iszerofunc`` is a callable that returns a Boolean that indicates
136 if its input is zero, or None if no such determination can be made.
137 ``simpfunc`` is a callable that simplifies its input. It must return
138 its input if it does not simplify its input. Passing in
139 ``simpfunc=None`` indicates that the pivot search should not attempt
140 to simplify any candidate pivots.
142 Returns a 4-tuple:
143 (pivot_offset, pivot_val, assumed_nonzero, newly_determined)
144 ``pivot_offset`` is the sequence index of the pivot.
145 ``pivot_val`` is the value of the pivot.
146 pivot_val and col[pivot_index] are equivalent, but will be different
147 when col[pivot_index] was simplified during the pivot search.
148 ``assumed_nonzero`` is a boolean indicating if the pivot cannot be
149 guaranteed to be zero. If assumed_nonzero is true, then the pivot
150 may or may not be non-zero. If assumed_nonzero is false, then
151 the pivot is non-zero.
152 ``newly_determined`` is a list of index-value pairs of pivot candidates
153 that were simplified during the pivot search.
154 """
156 # indeterminates holds the index-value pairs of each pivot candidate
157 # that is neither zero or non-zero, as determined by iszerofunc().
158 # If iszerofunc() indicates that a candidate pivot is guaranteed
159 # non-zero, or that every candidate pivot is zero then the contents
160 # of indeterminates are unused.
161 # Otherwise, the only viable candidate pivots are symbolic.
162 # In this case, indeterminates will have at least one entry,
163 # and all but the first entry are ignored when simpfunc is None.
164 indeterminates = []
165 for i, col_val in enumerate(col):
166 col_val_is_zero = iszerofunc(col_val)
167 if col_val_is_zero == False:
168 # This pivot candidate is non-zero.
169 return i, col_val, False, []
170 elif col_val_is_zero is None:
171 # The candidate pivot's comparison with zero
172 # is indeterminate.
173 indeterminates.append((i, col_val))
175 if len(indeterminates) == 0:
176 # All candidate pivots are guaranteed to be zero, i.e. there is
177 # no pivot.
178 return None, None, False, []
180 if simpfunc is None:
181 # Caller did not pass in a simplification function that might
182 # determine if an indeterminate pivot candidate is guaranteed
183 # to be nonzero, so assume the first indeterminate candidate
184 # is non-zero.
185 return indeterminates[0][0], indeterminates[0][1], True, []
187 # newly_determined holds index-value pairs of candidate pivots
188 # that were simplified during the search for a non-zero pivot.
189 newly_determined = []
190 for i, col_val in indeterminates:
191 tmp_col_val = simpfunc(col_val)
192 if id(col_val) != id(tmp_col_val):
193 # simpfunc() simplified this candidate pivot.
194 newly_determined.append((i, tmp_col_val))
195 if iszerofunc(tmp_col_val) == False:
196 # Candidate pivot simplified to a guaranteed non-zero value.
197 return i, tmp_col_val, False, newly_determined
199 return indeterminates[0][0], indeterminates[0][1], True, newly_determined
202# This functions is a candidate for caching if it gets implemented for matrices.
203def _berkowitz_toeplitz_matrix(M):
204 """Return (A,T) where T the Toeplitz matrix used in the Berkowitz algorithm
205 corresponding to ``M`` and A is the first principal submatrix.
206 """
208 # the 0 x 0 case is trivial
209 if M.rows == 0 and M.cols == 0:
210 return M._new(1,1, [M.one])
212 #
213 # Partition M = [ a_11 R ]
214 # [ C A ]
215 #
217 a, R = M[0,0], M[0, 1:]
218 C, A = M[1:, 0], M[1:,1:]
220 #
221 # The Toeplitz matrix looks like
222 #
223 # [ 1 ]
224 # [ -a 1 ]
225 # [ -RC -a 1 ]
226 # [ -RAC -RC -a 1 ]
227 # [ -RA**2C -RAC -RC -a 1 ]
228 # etc.
230 # Compute the diagonal entries.
231 # Because multiplying matrix times vector is so much
232 # more efficient than matrix times matrix, recursively
233 # compute -R * A**n * C.
234 diags = [C]
235 for i in range(M.rows - 2):
236 diags.append(A.multiply(diags[i], dotprodsimp=None))
237 diags = [(-R).multiply(d, dotprodsimp=None)[0, 0] for d in diags]
238 diags = [M.one, -a] + diags
240 def entry(i,j):
241 if j > i:
242 return M.zero
243 return diags[i - j]
245 toeplitz = M._new(M.cols + 1, M.rows, entry)
246 return (A, toeplitz)
249# This functions is a candidate for caching if it gets implemented for matrices.
250def _berkowitz_vector(M):
251 """ Run the Berkowitz algorithm and return a vector whose entries
252 are the coefficients of the characteristic polynomial of ``M``.
254 Given N x N matrix, efficiently compute
255 coefficients of characteristic polynomials of ``M``
256 without division in the ground domain.
258 This method is particularly useful for computing determinant,
259 principal minors and characteristic polynomial when ``M``
260 has complicated coefficients e.g. polynomials. Semi-direct
261 usage of this algorithm is also important in computing
262 efficiently sub-resultant PRS.
264 Assuming that M is a square matrix of dimension N x N and
265 I is N x N identity matrix, then the Berkowitz vector is
266 an N x 1 vector whose entries are coefficients of the
267 polynomial
269 charpoly(M) = det(t*I - M)
271 As a consequence, all polynomials generated by Berkowitz
272 algorithm are monic.
274 For more information on the implemented algorithm refer to:
276 [1] S.J. Berkowitz, On computing the determinant in small
277 parallel time using a small number of processors, ACM,
278 Information Processing Letters 18, 1984, pp. 147-150
280 [2] M. Keber, Division-Free computation of sub-resultants
281 using Bezout matrices, Tech. Report MPI-I-2006-1-006,
282 Saarbrucken, 2006
283 """
285 # handle the trivial cases
286 if M.rows == 0 and M.cols == 0:
287 return M._new(1, 1, [M.one])
288 elif M.rows == 1 and M.cols == 1:
289 return M._new(2, 1, [M.one, -M[0,0]])
291 submat, toeplitz = _berkowitz_toeplitz_matrix(M)
293 return toeplitz.multiply(_berkowitz_vector(submat), dotprodsimp=None)
296def _adjugate(M, method="berkowitz"):
297 """Returns the adjugate, or classical adjoint, of
298 a matrix. That is, the transpose of the matrix of cofactors.
300 https://en.wikipedia.org/wiki/Adjugate
302 Parameters
303 ==========
305 method : string, optional
306 Method to use to find the cofactors, can be "bareiss", "berkowitz" or
307 "lu".
309 Examples
310 ========
312 >>> from sympy import Matrix
313 >>> M = Matrix([[1, 2], [3, 4]])
314 >>> M.adjugate()
315 Matrix([
316 [ 4, -2],
317 [-3, 1]])
319 See Also
320 ========
322 cofactor_matrix
323 sympy.matrices.common.MatrixCommon.transpose
324 """
326 return M.cofactor_matrix(method=method).transpose()
329# This functions is a candidate for caching if it gets implemented for matrices.
330def _charpoly(M, x='lambda', simplify=_simplify):
331 """Computes characteristic polynomial det(x*I - M) where I is
332 the identity matrix.
334 A PurePoly is returned, so using different variables for ``x`` does
335 not affect the comparison or the polynomials:
337 Parameters
338 ==========
340 x : string, optional
341 Name for the "lambda" variable, defaults to "lambda".
343 simplify : function, optional
344 Simplification function to use on the characteristic polynomial
345 calculated. Defaults to ``simplify``.
347 Examples
348 ========
350 >>> from sympy import Matrix
351 >>> from sympy.abc import x, y
352 >>> M = Matrix([[1, 3], [2, 0]])
353 >>> M.charpoly()
354 PurePoly(lambda**2 - lambda - 6, lambda, domain='ZZ')
355 >>> M.charpoly(x) == M.charpoly(y)
356 True
357 >>> M.charpoly(x) == M.charpoly(y)
358 True
360 Specifying ``x`` is optional; a symbol named ``lambda`` is used by
361 default (which looks good when pretty-printed in unicode):
363 >>> M.charpoly().as_expr()
364 lambda**2 - lambda - 6
366 And if ``x`` clashes with an existing symbol, underscores will
367 be prepended to the name to make it unique:
369 >>> M = Matrix([[1, 2], [x, 0]])
370 >>> M.charpoly(x).as_expr()
371 _x**2 - _x - 2*x
373 Whether you pass a symbol or not, the generator can be obtained
374 with the gen attribute since it may not be the same as the symbol
375 that was passed:
377 >>> M.charpoly(x).gen
378 _x
379 >>> M.charpoly(x).gen == x
380 False
382 Notes
383 =====
385 The Samuelson-Berkowitz algorithm is used to compute
386 the characteristic polynomial efficiently and without any
387 division operations. Thus the characteristic polynomial over any
388 commutative ring without zero divisors can be computed.
390 If the determinant det(x*I - M) can be found out easily as
391 in the case of an upper or a lower triangular matrix, then
392 instead of Samuelson-Berkowitz algorithm, eigenvalues are computed
393 and the characteristic polynomial with their help.
395 See Also
396 ========
398 det
399 """
401 if not M.is_square:
402 raise NonSquareMatrixError()
403 if M.is_lower or M.is_upper:
404 diagonal_elements = M.diagonal()
405 x = uniquely_named_symbol(x, diagonal_elements, modify=lambda s: '_' + s)
406 m = 1
407 for i in diagonal_elements:
408 m = m * (x - simplify(i))
409 return PurePoly(m, x)
411 berk_vector = _berkowitz_vector(M)
412 x = uniquely_named_symbol(x, berk_vector, modify=lambda s: '_' + s)
414 return PurePoly([simplify(a) for a in berk_vector], x)
417def _cofactor(M, i, j, method="berkowitz"):
418 """Calculate the cofactor of an element.
420 Parameters
421 ==========
423 method : string, optional
424 Method to use to find the cofactors, can be "bareiss", "berkowitz" or
425 "lu".
427 Examples
428 ========
430 >>> from sympy import Matrix
431 >>> M = Matrix([[1, 2], [3, 4]])
432 >>> M.cofactor(0, 1)
433 -3
435 See Also
436 ========
438 cofactor_matrix
439 minor
440 minor_submatrix
441 """
443 if not M.is_square or M.rows < 1:
444 raise NonSquareMatrixError()
446 return S.NegativeOne**((i + j) % 2) * M.minor(i, j, method)
449def _cofactor_matrix(M, method="berkowitz"):
450 """Return a matrix containing the cofactor of each element.
452 Parameters
453 ==========
455 method : string, optional
456 Method to use to find the cofactors, can be "bareiss", "berkowitz" or
457 "lu".
459 Examples
460 ========
462 >>> from sympy import Matrix
463 >>> M = Matrix([[1, 2], [3, 4]])
464 >>> M.cofactor_matrix()
465 Matrix([
466 [ 4, -3],
467 [-2, 1]])
469 See Also
470 ========
472 cofactor
473 minor
474 minor_submatrix
475 """
477 if not M.is_square or M.rows < 1:
478 raise NonSquareMatrixError()
480 return M._new(M.rows, M.cols,
481 lambda i, j: M.cofactor(i, j, method))
483def _per(M):
484 """Returns the permanent of a matrix. Unlike determinant,
485 permanent is defined for both square and non-square matrices.
487 For an m x n matrix, with m less than or equal to n,
488 it is given as the sum over the permutations s of size
489 less than or equal to m on [1, 2, . . . n] of the product
490 from i = 1 to m of M[i, s[i]]. Taking the transpose will
491 not affect the value of the permanent.
493 In the case of a square matrix, this is the same as the permutation
494 definition of the determinant, but it does not take the sign of the
495 permutation into account. Computing the permanent with this definition
496 is quite inefficient, so here the Ryser formula is used.
498 Examples
499 ========
501 >>> from sympy import Matrix
502 >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
503 >>> M.per()
504 450
505 >>> M = Matrix([1, 5, 7])
506 >>> M.per()
507 13
509 References
510 ==========
512 .. [1] Prof. Frank Ben's notes: https://math.berkeley.edu/~bernd/ban275.pdf
513 .. [2] Wikipedia article on Permanent: https://en.wikipedia.org/wiki/Permanent_%28mathematics%29
514 .. [3] https://reference.wolfram.com/language/ref/Permanent.html
515 .. [4] Permanent of a rectangular matrix : https://arxiv.org/pdf/0904.3251.pdf
516 """
517 import itertools
519 m, n = M.shape
520 if m > n:
521 M = M.T
522 m, n = n, m
523 s = list(range(n))
525 subsets = []
526 for i in range(1, m + 1):
527 subsets += list(map(list, itertools.combinations(s, i)))
529 perm = 0
530 for subset in subsets:
531 prod = 1
532 sub_len = len(subset)
533 for i in range(m):
534 prod *= sum([M[i, j] for j in subset])
535 perm += prod * S.NegativeOne**sub_len * nC(n - sub_len, m - sub_len)
536 perm *= S.NegativeOne**m
537 return perm.simplify()
539def _det_DOM(M):
540 DOM = DomainMatrix.from_Matrix(M, field=True, extension=True)
541 K = DOM.domain
542 return K.to_sympy(DOM.det())
544# This functions is a candidate for caching if it gets implemented for matrices.
545def _det(M, method="bareiss", iszerofunc=None):
546 """Computes the determinant of a matrix if ``M`` is a concrete matrix object
547 otherwise return an expressions ``Determinant(M)`` if ``M`` is a
548 ``MatrixSymbol`` or other expression.
550 Parameters
551 ==========
553 method : string, optional
554 Specifies the algorithm used for computing the matrix determinant.
556 If the matrix is at most 3x3, a hard-coded formula is used and the
557 specified method is ignored. Otherwise, it defaults to
558 ``'bareiss'``.
560 Also, if the matrix is an upper or a lower triangular matrix, determinant
561 is computed by simple multiplication of diagonal elements, and the
562 specified method is ignored.
564 If it is set to ``'domain-ge'``, then Gaussian elimination method will
565 be used via using DomainMatrix.
567 If it is set to ``'bareiss'``, Bareiss' fraction-free algorithm will
568 be used.
570 If it is set to ``'berkowitz'``, Berkowitz' algorithm will be used.
572 Otherwise, if it is set to ``'lu'``, LU decomposition will be used.
574 .. note::
575 For backward compatibility, legacy keys like "bareis" and
576 "det_lu" can still be used to indicate the corresponding
577 methods.
578 And the keys are also case-insensitive for now. However, it is
579 suggested to use the precise keys for specifying the method.
581 iszerofunc : FunctionType or None, optional
582 If it is set to ``None``, it will be defaulted to ``_iszero`` if the
583 method is set to ``'bareiss'``, and ``_is_zero_after_expand_mul`` if
584 the method is set to ``'lu'``.
586 It can also accept any user-specified zero testing function, if it
587 is formatted as a function which accepts a single symbolic argument
588 and returns ``True`` if it is tested as zero and ``False`` if it
589 tested as non-zero, and also ``None`` if it is undecidable.
591 Returns
592 =======
594 det : Basic
595 Result of determinant.
597 Raises
598 ======
600 ValueError
601 If unrecognized keys are given for ``method`` or ``iszerofunc``.
603 NonSquareMatrixError
604 If attempted to calculate determinant from a non-square matrix.
606 Examples
607 ========
609 >>> from sympy import Matrix, eye, det
610 >>> I3 = eye(3)
611 >>> det(I3)
612 1
613 >>> M = Matrix([[1, 2], [3, 4]])
614 >>> det(M)
615 -2
616 >>> det(M) == M.det()
617 True
618 >>> M.det(method="domain-ge")
619 -2
620 """
622 # sanitize `method`
623 method = method.lower()
625 if method == "bareis":
626 method = "bareiss"
627 elif method == "det_lu":
628 method = "lu"
630 if method not in ("bareiss", "berkowitz", "lu", "domain-ge"):
631 raise ValueError("Determinant method '%s' unrecognized" % method)
633 if iszerofunc is None:
634 if method == "bareiss":
635 iszerofunc = _is_zero_after_expand_mul
636 elif method == "lu":
637 iszerofunc = _iszero
639 elif not isinstance(iszerofunc, FunctionType):
640 raise ValueError("Zero testing method '%s' unrecognized" % iszerofunc)
642 n = M.rows
644 if n == M.cols: # square check is done in individual method functions
645 if n == 0:
646 return M.one
647 elif n == 1:
648 return M[0, 0]
649 elif n == 2:
650 m = M[0, 0] * M[1, 1] - M[0, 1] * M[1, 0]
651 return _get_intermediate_simp(_dotprodsimp)(m)
652 elif n == 3:
653 m = (M[0, 0] * M[1, 1] * M[2, 2]
654 + M[0, 1] * M[1, 2] * M[2, 0]
655 + M[0, 2] * M[1, 0] * M[2, 1]
656 - M[0, 2] * M[1, 1] * M[2, 0]
657 - M[0, 0] * M[1, 2] * M[2, 1]
658 - M[0, 1] * M[1, 0] * M[2, 2])
659 return _get_intermediate_simp(_dotprodsimp)(m)
661 dets = []
662 for b in M.strongly_connected_components():
663 if method == "domain-ge": # uses DomainMatrix to evaluate determinant
664 det = _det_DOM(M[b, b])
665 elif method == "bareiss":
666 det = M[b, b]._eval_det_bareiss(iszerofunc=iszerofunc)
667 elif method == "berkowitz":
668 det = M[b, b]._eval_det_berkowitz()
669 elif method == "lu":
670 det = M[b, b]._eval_det_lu(iszerofunc=iszerofunc)
671 dets.append(det)
672 return Mul(*dets)
675# This functions is a candidate for caching if it gets implemented for matrices.
676def _det_bareiss(M, iszerofunc=_is_zero_after_expand_mul):
677 """Compute matrix determinant using Bareiss' fraction-free
678 algorithm which is an extension of the well known Gaussian
679 elimination method. This approach is best suited for dense
680 symbolic matrices and will result in a determinant with
681 minimal number of fractions. It means that less term
682 rewriting is needed on resulting formulae.
684 Parameters
685 ==========
687 iszerofunc : function, optional
688 The function to use to determine zeros when doing an LU decomposition.
689 Defaults to ``lambda x: x.is_zero``.
691 TODO: Implement algorithm for sparse matrices (SFF),
692 http://www.eecis.udel.edu/~saunders/papers/sffge/it5.ps.
693 """
695 # Recursively implemented Bareiss' algorithm as per Deanna Richelle Leggett's
696 # thesis http://www.math.usm.edu/perry/Research/Thesis_DRL.pdf
697 def bareiss(mat, cumm=1):
698 if mat.rows == 0:
699 return mat.one
700 elif mat.rows == 1:
701 return mat[0, 0]
703 # find a pivot and extract the remaining matrix
704 # With the default iszerofunc, _find_reasonable_pivot slows down
705 # the computation by the factor of 2.5 in one test.
706 # Relevant issues: #10279 and #13877.
707 pivot_pos, pivot_val, _, _ = _find_reasonable_pivot(mat[:, 0], iszerofunc=iszerofunc)
708 if pivot_pos is None:
709 return mat.zero
711 # if we have a valid pivot, we'll do a "row swap", so keep the
712 # sign of the det
713 sign = (-1) ** (pivot_pos % 2)
715 # we want every row but the pivot row and every column
716 rows = [i for i in range(mat.rows) if i != pivot_pos]
717 cols = list(range(mat.cols))
718 tmp_mat = mat.extract(rows, cols)
720 def entry(i, j):
721 ret = (pivot_val*tmp_mat[i, j + 1] - mat[pivot_pos, j + 1]*tmp_mat[i, 0]) / cumm
722 if _get_intermediate_simp_bool(True):
723 return _dotprodsimp(ret)
724 elif not ret.is_Atom:
725 return cancel(ret)
726 return ret
728 return sign*bareiss(M._new(mat.rows - 1, mat.cols - 1, entry), pivot_val)
730 if not M.is_square:
731 raise NonSquareMatrixError()
733 if M.rows == 0:
734 return M.one
735 # sympy/matrices/tests/test_matrices.py contains a test that
736 # suggests that the determinant of a 0 x 0 matrix is one, by
737 # convention.
739 return bareiss(M)
742def _det_berkowitz(M):
743 """ Use the Berkowitz algorithm to compute the determinant."""
745 if not M.is_square:
746 raise NonSquareMatrixError()
748 if M.rows == 0:
749 return M.one
750 # sympy/matrices/tests/test_matrices.py contains a test that
751 # suggests that the determinant of a 0 x 0 matrix is one, by
752 # convention.
754 berk_vector = _berkowitz_vector(M)
755 return (-1)**(len(berk_vector) - 1) * berk_vector[-1]
758# This functions is a candidate for caching if it gets implemented for matrices.
759def _det_LU(M, iszerofunc=_iszero, simpfunc=None):
760 """ Computes the determinant of a matrix from its LU decomposition.
761 This function uses the LU decomposition computed by
762 LUDecomposition_Simple().
764 The keyword arguments iszerofunc and simpfunc are passed to
765 LUDecomposition_Simple().
766 iszerofunc is a callable that returns a boolean indicating if its
767 input is zero, or None if it cannot make the determination.
768 simpfunc is a callable that simplifies its input.
769 The default is simpfunc=None, which indicate that the pivot search
770 algorithm should not attempt to simplify any candidate pivots.
771 If simpfunc fails to simplify its input, then it must return its input
772 instead of a copy.
774 Parameters
775 ==========
777 iszerofunc : function, optional
778 The function to use to determine zeros when doing an LU decomposition.
779 Defaults to ``lambda x: x.is_zero``.
781 simpfunc : function, optional
782 The simplification function to use when looking for zeros for pivots.
783 """
785 if not M.is_square:
786 raise NonSquareMatrixError()
788 if M.rows == 0:
789 return M.one
790 # sympy/matrices/tests/test_matrices.py contains a test that
791 # suggests that the determinant of a 0 x 0 matrix is one, by
792 # convention.
794 lu, row_swaps = M.LUdecomposition_Simple(iszerofunc=iszerofunc,
795 simpfunc=simpfunc)
796 # P*A = L*U => det(A) = det(L)*det(U)/det(P) = det(P)*det(U).
797 # Lower triangular factor L encoded in lu has unit diagonal => det(L) = 1.
798 # P is a permutation matrix => det(P) in {-1, 1} => 1/det(P) = det(P).
799 # LUdecomposition_Simple() returns a list of row exchange index pairs, rather
800 # than a permutation matrix, but det(P) = (-1)**len(row_swaps).
802 # Avoid forming the potentially time consuming product of U's diagonal entries
803 # if the product is zero.
804 # Bottom right entry of U is 0 => det(A) = 0.
805 # It may be impossible to determine if this entry of U is zero when it is symbolic.
806 if iszerofunc(lu[lu.rows-1, lu.rows-1]):
807 return M.zero
809 # Compute det(P)
810 det = -M.one if len(row_swaps)%2 else M.one
812 # Compute det(U) by calculating the product of U's diagonal entries.
813 # The upper triangular portion of lu is the upper triangular portion of the
814 # U factor in the LU decomposition.
815 for k in range(lu.rows):
816 det *= lu[k, k]
818 # return det(P)*det(U)
819 return det
822def _minor(M, i, j, method="berkowitz"):
823 """Return the (i,j) minor of ``M``. That is,
824 return the determinant of the matrix obtained by deleting
825 the `i`th row and `j`th column from ``M``.
827 Parameters
828 ==========
830 i, j : int
831 The row and column to exclude to obtain the submatrix.
833 method : string, optional
834 Method to use to find the determinant of the submatrix, can be
835 "bareiss", "berkowitz" or "lu".
837 Examples
838 ========
840 >>> from sympy import Matrix
841 >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
842 >>> M.minor(1, 1)
843 -12
845 See Also
846 ========
848 minor_submatrix
849 cofactor
850 det
851 """
853 if not M.is_square:
854 raise NonSquareMatrixError()
856 return M.minor_submatrix(i, j).det(method=method)
859def _minor_submatrix(M, i, j):
860 """Return the submatrix obtained by removing the `i`th row
861 and `j`th column from ``M`` (works with Pythonic negative indices).
863 Parameters
864 ==========
866 i, j : int
867 The row and column to exclude to obtain the submatrix.
869 Examples
870 ========
872 >>> from sympy import Matrix
873 >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
874 >>> M.minor_submatrix(1, 1)
875 Matrix([
876 [1, 3],
877 [7, 9]])
879 See Also
880 ========
882 minor
883 cofactor
884 """
886 if i < 0:
887 i += M.rows
888 if j < 0:
889 j += M.cols
891 if not 0 <= i < M.rows or not 0 <= j < M.cols:
892 raise ValueError("`i` and `j` must satisfy 0 <= i < ``M.rows`` "
893 "(%d)" % M.rows + "and 0 <= j < ``M.cols`` (%d)." % M.cols)
895 rows = [a for a in range(M.rows) if a != i]
896 cols = [a for a in range(M.cols) if a != j]
898 return M.extract(rows, cols)