2025-summer/team-10
10
Fork 0
Code Issues Pull requests Projects Releases Activity
team-10/venv/Lib/site-packages/sympy/polys/matrices/domainmatrix.py

3984 lines
113 KiB
Python
Raw Normal View History

Adding all project files
2025-08-02 02:00:33 +02:00
"""
Module for the DomainMatrix class.
A DomainMatrix represents a matrix with elements that are in a particular
Domain. Each DomainMatrix internally wraps a DDM which is used for the
lower-level operations. The idea is that the DomainMatrix class provides the
convenience routines for converting between Expr and the poly domains as well
as unifying matrices with different domains.
"""
from __future__ import annotations
from collections import Counter
from functools import reduce
from sympy.external.gmpy import GROUND_TYPES
from sympy.utilities.decorator import doctest_depends_on
from sympy.core.sympify import _sympify
from ..domains import Domain
from ..constructor import construct_domain
from .exceptions import (
DMFormatError,
DMBadInputError,
DMShapeError,
DMDomainError,
DMNotAField,
DMNonSquareMatrixError,
DMNonInvertibleMatrixError
)
from .domainscalar import DomainScalar
from sympy.polys.domains import ZZ, EXRAW, QQ
from sympy.polys.densearith import dup_mul
from sympy.polys.densebasic import dup_convert
from sympy.polys.densetools import (
dup_mul_ground,
dup_quo_ground,
dup_content,
dup_clear_denoms,
dup_primitive,
dup_transform,
)
from sympy.polys.factortools import dup_factor_list
from sympy.polys.polyutils import _sort_factors
from .ddm import DDM
from .sdm import SDM
from .dfm import DFM
from .rref import _dm_rref, _dm_rref_den
if GROUND_TYPES != 'flint':
__doctest_skip__ = ['DomainMatrix.to_dfm', 'DomainMatrix.to_dfm_or_ddm']
else:
__doctest_skip__ = ['DomainMatrix.from_list']
def DM(rows, domain):
"""Convenient alias for DomainMatrix.from_list
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DM
>>> DM([[1, 2], [3, 4]], ZZ)
DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ)
See Also
========
DomainMatrix.from_list
"""
return DomainMatrix.from_list(rows, domain)
class DomainMatrix:
r"""
Associate Matrix with :py:class:`~.Domain`
Explanation
===========
DomainMatrix uses :py:class:`~.Domain` for its internal representation
which makes it faster than the SymPy Matrix class (currently) for many
common operations, but this advantage makes it not entirely compatible
with Matrix. DomainMatrix are analogous to numpy arrays with "dtype".
In the DomainMatrix, each element has a domain such as :ref:`ZZ`
or :ref:`QQ(a)`.
Examples
========
Creating a DomainMatrix from the existing Matrix class:
>>> from sympy import Matrix
>>> from sympy.polys.matrices import DomainMatrix
>>> Matrix1 = Matrix([
... [1, 2],
... [3, 4]])
>>> A = DomainMatrix.from_Matrix(Matrix1)
>>> A
DomainMatrix({0: {0: 1, 1: 2}, 1: {0: 3, 1: 4}}, (2, 2), ZZ)
Directly forming a DomainMatrix:
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> A
DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ)
See Also
========
DDM
SDM
Domain
Poly
"""
rep: SDM | DDM | DFM
shape: tuple[int, int]
domain: Domain
def __new__(cls, rows, shape, domain, *, fmt=None):
"""
Creates a :py:class:`~.DomainMatrix`.
Parameters
==========
rows : Represents elements of DomainMatrix as list of lists
shape : Represents dimension of DomainMatrix
domain : Represents :py:class:`~.Domain` of DomainMatrix
Raises
======
TypeError
If any of rows, shape and domain are not provided
"""
if isinstance(rows, (DDM, SDM, DFM)):
raise TypeError("Use from_rep to initialise from SDM/DDM")
elif isinstance(rows, list):
rep = DDM(rows, shape, domain)
elif isinstance(rows, dict):
rep = SDM(rows, shape, domain)
else:
msg = "Input should be list-of-lists or dict-of-dicts"
raise TypeError(msg)
if fmt is not None:
if fmt == 'sparse':
rep = rep.to_sdm()
elif fmt == 'dense':
rep = rep.to_ddm()
else:
raise ValueError("fmt should be 'sparse' or 'dense'")
# Use python-flint for dense matrices if possible
if rep.fmt == 'dense' and DFM._supports_domain(domain):
rep = rep.to_dfm()
return cls.from_rep(rep)
def __reduce__(self):
rep = self.rep
if rep.fmt == 'dense':
arg = self.to_list()
elif rep.fmt == 'sparse':
arg = dict(rep)
else:
raise RuntimeError # pragma: no cover
args = (arg, rep.shape, rep.domain)
return (self.__class__, args)
def __getitem__(self, key):
i, j = key
m, n = self.shape
if not (isinstance(i, slice) or isinstance(j, slice)):
return DomainScalar(self.rep.getitem(i, j), self.domain)
if not isinstance(i, slice):
if not -m <= i < m:
raise IndexError("Row index out of range")
i = i % m
i = slice(i, i+1)
if not isinstance(j, slice):
if not -n <= j < n:
raise IndexError("Column index out of range")
j = j % n
j = slice(j, j+1)
return self.from_rep(self.rep.extract_slice(i, j))
def getitem_sympy(self, i, j):
return self.domain.to_sympy(self.rep.getitem(i, j))
def extract(self, rowslist, colslist):
return self.from_rep(self.rep.extract(rowslist, colslist))
def __setitem__(self, key, value):
i, j = key
if not self.domain.of_type(value):
raise TypeError
if isinstance(i, int) and isinstance(j, int):
self.rep.setitem(i, j, value)
else:
raise NotImplementedError
@classmethod
def from_rep(cls, rep):
"""Create a new DomainMatrix efficiently from DDM/SDM.
Examples
========
Create a :py:class:`~.DomainMatrix` with an dense internal
representation as :py:class:`~.DDM`:
>>> from sympy.polys.domains import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy.polys.matrices.ddm import DDM
>>> drep = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> dM = DomainMatrix.from_rep(drep)
>>> dM
DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ)
Create a :py:class:`~.DomainMatrix` with a sparse internal
representation as :py:class:`~.SDM`:
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy.polys.matrices.sdm import SDM
>>> from sympy import ZZ
>>> drep = SDM({0:{1:ZZ(1)},1:{0:ZZ(2)}}, (2, 2), ZZ)
>>> dM = DomainMatrix.from_rep(drep)
>>> dM
DomainMatrix({0: {1: 1}, 1: {0: 2}}, (2, 2), ZZ)
Parameters
==========
rep: SDM or DDM
The internal sparse or dense representation of the matrix.
Returns
=======
DomainMatrix
A :py:class:`~.DomainMatrix` wrapping *rep*.
Notes
=====
This takes ownership of rep as its internal representation. If rep is
being mutated elsewhere then a copy should be provided to
``from_rep``. Only minimal verification or checking is done on *rep*
as this is supposed to be an efficient internal routine.
"""
if not (isinstance(rep, (DDM, SDM)) or (DFM is not None and isinstance(rep, DFM))):
raise TypeError("rep should be of type DDM or SDM")
self = super().__new__(cls)
self.rep = rep
self.shape = rep.shape
self.domain = rep.domain
return self
@classmethod
@doctest_depends_on(ground_types=['python', 'gmpy'])
def from_list(cls, rows, domain):
r"""
Convert a list of lists into a DomainMatrix
Parameters
==========
rows: list of lists
Each element of the inner lists should be either the single arg,
or tuple of args, that would be passed to the domain constructor
in order to form an element of the domain. See examples.
Returns
=======
DomainMatrix containing elements defined in rows
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy import FF, QQ, ZZ
>>> A = DomainMatrix.from_list([[1, 0, 1], [0, 0, 1]], ZZ)
>>> A
DomainMatrix([[1, 0, 1], [0, 0, 1]], (2, 3), ZZ)
>>> B = DomainMatrix.from_list([[1, 0, 1], [0, 0, 1]], FF(7))
>>> B
DomainMatrix([[1 mod 7, 0 mod 7, 1 mod 7], [0 mod 7, 0 mod 7, 1 mod 7]], (2, 3), GF(7))
>>> C = DomainMatrix.from_list([[(1, 2), (3, 1)], [(1, 4), (5, 1)]], QQ)
>>> C
DomainMatrix([[1/2, 3], [1/4, 5]], (2, 2), QQ)
See Also
========
from_list_sympy
"""
nrows = len(rows)
ncols = 0 if not nrows else len(rows[0])
conv = lambda e: domain(*e) if isinstance(e, tuple) else domain(e)
domain_rows = [[conv(e) for e in row] for row in rows]
return DomainMatrix(domain_rows, (nrows, ncols), domain)
@classmethod
def from_list_sympy(cls, nrows, ncols, rows, **kwargs):
r"""
Convert a list of lists of Expr into a DomainMatrix using construct_domain
Parameters
==========
nrows: number of rows
ncols: number of columns
rows: list of lists
Returns
=======
DomainMatrix containing elements of rows
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy.abc import x, y, z
>>> A = DomainMatrix.from_list_sympy(1, 3, [[x, y, z]])
>>> A
DomainMatrix([[x, y, z]], (1, 3), ZZ[x,y,z])
See Also
========
sympy.polys.constructor.construct_domain, from_dict_sympy
"""
assert len(rows) == nrows
assert all(len(row) == ncols for row in rows)
items_sympy = [_sympify(item) for row in rows for item in row]
domain, items_domain = cls.get_domain(items_sympy, **kwargs)
domain_rows = [[items_domain[ncols*r + c] for c in range(ncols)] for r in range(nrows)]
return DomainMatrix(domain_rows, (nrows, ncols), domain)
@classmethod
def from_dict_sympy(cls, nrows, ncols, elemsdict, **kwargs):
"""
Parameters
==========
nrows: number of rows
ncols: number of cols
elemsdict: dict of dicts containing non-zero elements of the DomainMatrix
Returns
=======
DomainMatrix containing elements of elemsdict
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy.abc import x,y,z
>>> elemsdict = {0: {0:x}, 1:{1: y}, 2: {2: z}}
>>> A = DomainMatrix.from_dict_sympy(3, 3, elemsdict)
>>> A
DomainMatrix({0: {0: x}, 1: {1: y}, 2: {2: z}}, (3, 3), ZZ[x,y,z])
See Also
========
from_list_sympy
"""
if not all(0 <= r < nrows for r in elemsdict):
raise DMBadInputError("Row out of range")
if not all(0 <= c < ncols for row in elemsdict.values() for c in row):
raise DMBadInputError("Column out of range")
items_sympy = [_sympify(item) for row in elemsdict.values() for item in row.values()]
domain, items_domain = cls.get_domain(items_sympy, **kwargs)
idx = 0
items_dict = {}
for i, row in elemsdict.items():
items_dict[i] = {}
for j in row:
items_dict[i][j] = items_domain[idx]
idx += 1
return DomainMatrix(items_dict, (nrows, ncols), domain)
@classmethod
def from_Matrix(cls, M, fmt='sparse',**kwargs):
r"""
Convert Matrix to DomainMatrix
Parameters
==========
M: Matrix
Returns
=======
Returns DomainMatrix with identical elements as M
Examples
========
>>> from sympy import Matrix
>>> from sympy.polys.matrices import DomainMatrix
>>> M = Matrix([
... [1.0, 3.4],
... [2.4, 1]])
>>> A = DomainMatrix.from_Matrix(M)
>>> A
DomainMatrix({0: {0: 1.0, 1: 3.4}, 1: {0: 2.4, 1: 1.0}}, (2, 2), RR)
We can keep internal representation as ddm using fmt='dense'
>>> from sympy import Matrix, QQ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix.from_Matrix(Matrix([[QQ(1, 2), QQ(3, 4)], [QQ(0, 1), QQ(0, 1)]]), fmt='dense')
>>> A.rep
[[1/2, 3/4], [0, 0]]
See Also
========
Matrix
"""
if fmt == 'dense':
return cls.from_list_sympy(*M.shape, M.tolist(), **kwargs)
return cls.from_dict_sympy(*M.shape, M.todod(), **kwargs)
@classmethod
def get_domain(cls, items_sympy, **kwargs):
K, items_K = construct_domain(items_sympy, **kwargs)
return K, items_K
def choose_domain(self, **opts):
"""Convert to a domain found by :func:`~.construct_domain`.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DM
>>> M = DM([[1, 2], [3, 4]], ZZ)
>>> M
DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ)
>>> M.choose_domain(field=True)
DomainMatrix([[1, 2], [3, 4]], (2, 2), QQ)
>>> from sympy.abc import x
>>> M = DM([[1, x], [x**2, x**3]], ZZ[x])
>>> M.choose_domain(field=True).domain
ZZ(x)
Keyword arguments are passed to :func:`~.construct_domain`.
See Also
========
construct_domain
convert_to
"""
elements, data = self.to_sympy().to_flat_nz()
dom, elements_dom = construct_domain(elements, **opts)
return self.from_flat_nz(elements_dom, data, dom)
def copy(self):
return self.from_rep(self.rep.copy())
def convert_to(self, K):
r"""
Change the domain of DomainMatrix to desired domain or field
Parameters
==========
K : Represents the desired domain or field.
Alternatively, ``None`` may be passed, in which case this method
just returns a copy of this DomainMatrix.
Returns
=======
DomainMatrix
DomainMatrix with the desired domain or field
Examples
========
>>> from sympy import ZZ, ZZ_I
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> A.convert_to(ZZ_I)
DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ_I)
"""
if K == self.domain:
return self.copy()
rep = self.rep
# The DFM, DDM and SDM types do not do any implicit conversions so we
# manage switching between DDM and DFM here.
if rep.is_DFM and not DFM._supports_domain(K):
rep_K = rep.to_ddm().convert_to(K)
elif rep.is_DDM and DFM._supports_domain(K):
rep_K = rep.convert_to(K).to_dfm()
else:
rep_K = rep.convert_to(K)
return self.from_rep(rep_K)
def to_sympy(self):
return self.convert_to(EXRAW)
def to_field(self):
r"""
Returns a DomainMatrix with the appropriate field
Returns
=======
DomainMatrix
DomainMatrix with the appropriate field
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> A.to_field()
DomainMatrix([[1, 2], [3, 4]], (2, 2), QQ)
"""
K = self.domain.get_field()
return self.convert_to(K)
def to_sparse(self):
"""
Return a sparse DomainMatrix representation of *self*.
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy import QQ
>>> A = DomainMatrix([[1, 0],[0, 2]], (2, 2), QQ)
>>> A.rep
[[1, 0], [0, 2]]
>>> B = A.to_sparse()
>>> B.rep
{0: {0: 1}, 1: {1: 2}}
"""
if self.rep.fmt == 'sparse':
return self
return self.from_rep(self.rep.to_sdm())
def to_dense(self):
"""
Return a dense DomainMatrix representation of *self*.
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy import QQ
>>> A = DomainMatrix({0: {0: 1}, 1: {1: 2}}, (2, 2), QQ)
>>> A.rep
{0: {0: 1}, 1: {1: 2}}
>>> B = A.to_dense()
>>> B.rep
[[1, 0], [0, 2]]
"""
rep = self.rep
if rep.fmt == 'dense':
return self
return self.from_rep(rep.to_dfm_or_ddm())
def to_ddm(self):
"""
Return a :class:`~.DDM` representation of *self*.
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy import QQ
>>> A = DomainMatrix({0: {0: 1}, 1: {1: 2}}, (2, 2), QQ)
>>> ddm = A.to_ddm()
>>> ddm
[[1, 0], [0, 2]]
>>> type(ddm)
<class 'sympy.polys.matrices.ddm.DDM'>
See Also
========
to_sdm
to_dense
sympy.polys.matrices.ddm.DDM.to_sdm
"""
return self.rep.to_ddm()
def to_sdm(self):
"""
Return a :class:`~.SDM` representation of *self*.
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy import QQ
>>> A = DomainMatrix([[1, 0],[0, 2]], (2, 2), QQ)
>>> sdm = A.to_sdm()
>>> sdm
{0: {0: 1}, 1: {1: 2}}
>>> type(sdm)
<class 'sympy.polys.matrices.sdm.SDM'>
See Also
========
to_ddm
to_sparse
sympy.polys.matrices.sdm.SDM.to_ddm
"""
return self.rep.to_sdm()
@doctest_depends_on(ground_types=['flint'])
def to_dfm(self):
"""
Return a :class:`~.DFM` representation of *self*.
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy import QQ
>>> A = DomainMatrix([[1, 0],[0, 2]], (2, 2), QQ)
>>> dfm = A.to_dfm()
>>> dfm
[[1, 0], [0, 2]]
>>> type(dfm)
<class 'sympy.polys.matrices._dfm.DFM'>
See Also
========
to_ddm
to_dense
DFM
"""
return self.rep.to_dfm()
@doctest_depends_on(ground_types=['flint'])
def to_dfm_or_ddm(self):
"""
Return a :class:`~.DFM` or :class:`~.DDM` representation of *self*.
Explanation
===========
The :class:`~.DFM` representation can only be used if the ground types
are ``flint`` and the ground domain is supported by ``python-flint``.
This method will return a :class:`~.DFM` representation if possible,
but will return a :class:`~.DDM` representation otherwise.
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy import QQ
>>> A = DomainMatrix([[1, 0],[0, 2]], (2, 2), QQ)
>>> dfm = A.to_dfm_or_ddm()
>>> dfm
[[1, 0], [0, 2]]
>>> type(dfm) # Depends on the ground domain and ground types
<class 'sympy.polys.matrices._dfm.DFM'>
See Also
========
to_ddm: Always return a :class:`~.DDM` representation.
to_dfm: Returns a :class:`~.DFM` representation or raise an error.
to_dense: Convert internally to a :class:`~.DFM` or :class:`~.DDM`
DFM: The :class:`~.DFM` dense FLINT matrix representation.
DDM: The Python :class:`~.DDM` dense domain matrix representation.
"""
return self.rep.to_dfm_or_ddm()
@classmethod
def _unify_domain(cls, *matrices):
"""Convert matrices to a common domain"""
domains = {matrix.domain for matrix in matrices}
if len(domains) == 1:
return matrices
domain = reduce(lambda x, y: x.unify(y), domains)
return tuple(matrix.convert_to(domain) for matrix in matrices)
@classmethod
def _unify_fmt(cls, *matrices, fmt=None):
"""Convert matrices to the same format.
If all matrices have the same format, then return unmodified.
Otherwise convert both to the preferred format given as *fmt* which
should be 'dense' or 'sparse'.
"""
formats = {matrix.rep.fmt for matrix in matrices}
if len(formats) == 1:
return matrices
if fmt == 'sparse':
return tuple(matrix.to_sparse() for matrix in matrices)
elif fmt == 'dense':
return tuple(matrix.to_dense() for matrix in matrices)
else:
raise ValueError("fmt should be 'sparse' or 'dense'")
def unify(self, *others, fmt=None):
"""
Unifies the domains and the format of self and other
matrices.
Parameters
==========
others : DomainMatrix
fmt: string 'dense', 'sparse' or `None` (default)
The preferred format to convert to if self and other are not
already in the same format. If `None` or not specified then no
conversion if performed.
Returns
=======
Tuple[DomainMatrix]
Matrices with unified domain and format
Examples
========
Unify the domain of DomainMatrix that have different domains:
>>> from sympy import ZZ, QQ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([[ZZ(1), ZZ(2)]], (1, 2), ZZ)
>>> B = DomainMatrix([[QQ(1, 2), QQ(2)]], (1, 2), QQ)
>>> Aq, Bq = A.unify(B)
>>> Aq
DomainMatrix([[1, 2]], (1, 2), QQ)
>>> Bq
DomainMatrix([[1/2, 2]], (1, 2), QQ)
Unify the format (dense or sparse):
>>> A = DomainMatrix([[ZZ(1), ZZ(2)]], (1, 2), ZZ)
>>> B = DomainMatrix({0:{0: ZZ(1)}}, (2, 2), ZZ)
>>> B.rep
{0: {0: 1}}
>>> A2, B2 = A.unify(B, fmt='dense')
>>> B2.rep
[[1, 0], [0, 0]]
See Also
========
convert_to, to_dense, to_sparse
"""
matrices = (self,) + others
matrices = DomainMatrix._unify_domain(*matrices)
if fmt is not None:
matrices = DomainMatrix._unify_fmt(*matrices, fmt=fmt)
return matrices
def to_Matrix(self):
r"""
Convert DomainMatrix to Matrix
Returns
=======
Matrix
MutableDenseMatrix for the DomainMatrix
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> A.to_Matrix()
Matrix([
[1, 2],
[3, 4]])
See Also
========
from_Matrix
"""
from sympy.matrices.dense import MutableDenseMatrix
# XXX: If the internal representation of RepMatrix changes then this
# might need to be changed also.
if self.domain in (ZZ, QQ, EXRAW):
if self.rep.fmt == "sparse":
rep = self.copy()
else:
rep = self.to_sparse()
else:
rep = self.convert_to(EXRAW).to_sparse()
return MutableDenseMatrix._fromrep(rep)
def to_list(self):
"""
Convert :class:`DomainMatrix` to list of lists.
See Also
========
from_list
to_list_flat
to_flat_nz
to_dok
"""
return self.rep.to_list()
def to_list_flat(self):
"""
Convert :class:`DomainMatrix` to flat list.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> A.to_list_flat()
[1, 2, 3, 4]
See Also
========
from_list_flat
to_list
to_flat_nz
to_dok
"""
return self.rep.to_list_flat()
@classmethod
def from_list_flat(cls, elements, shape, domain):
"""
Create :class:`DomainMatrix` from flat list.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> element_list = [ZZ(1), ZZ(2), ZZ(3), ZZ(4)]
>>> A = DomainMatrix.from_list_flat(element_list, (2, 2), ZZ)
>>> A
DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ)
>>> A == A.from_list_flat(A.to_list_flat(), A.shape, A.domain)
True
See Also
========
to_list_flat
"""
ddm = DDM.from_list_flat(elements, shape, domain)
return cls.from_rep(ddm.to_dfm_or_ddm())
def to_flat_nz(self):
"""
Convert :class:`DomainMatrix` to list of nonzero elements and data.
Explanation
===========
Returns a tuple ``(elements, data)`` where ``elements`` is a list of
elements of the matrix with zeros possibly excluded. The matrix can be
reconstructed by passing these to :meth:`from_flat_nz`. The idea is to
be able to modify a flat list of the elements and then create a new
matrix of the same shape with the modified elements in the same
positions.
The format of ``data`` differs depending on whether the underlying
representation is dense or sparse but either way it represents the
positions of the elements in the list in a way that
:meth:`from_flat_nz` can use to reconstruct the matrix. The
:meth:`from_flat_nz` method should be called on the same
:class:`DomainMatrix` that was used to call :meth:`to_flat_nz`.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> elements, data = A.to_flat_nz()
>>> elements
[1, 2, 3, 4]
>>> A == A.from_flat_nz(elements, data, A.domain)
True
Create a matrix with the elements doubled:
>>> elements_doubled = [2*x for x in elements]
>>> A2 = A.from_flat_nz(elements_doubled, data, A.domain)
>>> A2 == 2*A
True
See Also
========
from_flat_nz
"""
return self.rep.to_flat_nz()
def from_flat_nz(self, elements, data, domain):
"""
Reconstruct :class:`DomainMatrix` after calling :meth:`to_flat_nz`.
See :meth:`to_flat_nz` for explanation.
See Also
========
to_flat_nz
"""
rep = self.rep.from_flat_nz(elements, data, domain)
return self.from_rep(rep)
def to_dod(self):
"""
Convert :class:`DomainMatrix` to dictionary of dictionaries (dod) format.
Explanation
===========
Returns a dictionary of dictionaries representing the matrix.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DM
>>> A = DM([[ZZ(1), ZZ(2), ZZ(0)], [ZZ(3), ZZ(0), ZZ(4)]], ZZ)
>>> A.to_dod()
{0: {0: 1, 1: 2}, 1: {0: 3, 2: 4}}
>>> A.to_sparse() == A.from_dod(A.to_dod(), A.shape, A.domain)
True
>>> A == A.from_dod_like(A.to_dod())
True
See Also
========
from_dod
from_dod_like
to_dok
to_list
to_list_flat
to_flat_nz
sympy.matrices.matrixbase.MatrixBase.todod
"""
return self.rep.to_dod()
@classmethod
def from_dod(cls, dod, shape, domain):
"""
Create sparse :class:`DomainMatrix` from dict of dict (dod) format.
See :meth:`to_dod` for explanation.
See Also
========
to_dod
from_dod_like
"""
return cls.from_rep(SDM.from_dod(dod, shape, domain))
def from_dod_like(self, dod, domain=None):
"""
Create :class:`DomainMatrix` like ``self`` from dict of dict (dod) format.
See :meth:`to_dod` for explanation.
See Also
========
to_dod
from_dod
"""
if domain is None:
domain = self.domain
return self.from_rep(self.rep.from_dod(dod, self.shape, domain))
def to_dok(self):
"""
Convert :class:`DomainMatrix` to dictionary of keys (dok) format.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(0)],
... [ZZ(0), ZZ(4)]], (2, 2), ZZ)
>>> A.to_dok()
{(0, 0): 1, (1, 1): 4}
The matrix can be reconstructed by calling :meth:`from_dok` although
the reconstructed matrix will always be in sparse format:
>>> A.to_sparse() == A.from_dok(A.to_dok(), A.shape, A.domain)
True
See Also
========
from_dok
to_list
to_list_flat
to_flat_nz
"""
return self.rep.to_dok()
@classmethod
def from_dok(cls, dok, shape, domain):
"""
Create :class:`DomainMatrix` from dictionary of keys (dok) format.
See :meth:`to_dok` for explanation.
See Also
========
to_dok
"""
return cls.from_rep(SDM.from_dok(dok, shape, domain))
def iter_values(self):
"""
Iterate over nonzero elements of the matrix.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([[ZZ(1), ZZ(0)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> list(A.iter_values())
[1, 3, 4]
See Also
========
iter_items
to_list_flat
sympy.matrices.matrixbase.MatrixBase.iter_values
"""
return self.rep.iter_values()
def iter_items(self):
"""
Iterate over indices and values of nonzero elements of the matrix.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([[ZZ(1), ZZ(0)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> list(A.iter_items())
[((0, 0), 1), ((1, 0), 3), ((1, 1), 4)]
See Also
========
iter_values
to_dok
sympy.matrices.matrixbase.MatrixBase.iter_items
"""
return self.rep.iter_items()
def nnz(self):
"""
Number of nonzero elements in the matrix.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DM
>>> A = DM([[1, 0], [0, 4]], ZZ)
>>> A.nnz()
2
"""
return self.rep.nnz()
def __repr__(self):
return 'DomainMatrix(%s, %r, %r)' % (str(self.rep), self.shape, self.domain)
def transpose(self):
"""Matrix transpose of ``self``"""
return self.from_rep(self.rep.transpose())
def flat(self):
rows, cols = self.shape
return [self[i,j].element for i in range(rows) for j in range(cols)]
@property
def is_zero_matrix(self):
return self.rep.is_zero_matrix()
@property
def is_upper(self):
"""
Says whether this matrix is upper-triangular. True can be returned
even if the matrix is not square.
"""
return self.rep.is_upper()
@property
def is_lower(self):
"""
Says whether this matrix is lower-triangular. True can be returned
even if the matrix is not square.
"""
return self.rep.is_lower()
@property
def is_diagonal(self):
"""
True if the matrix is diagonal.
Can return true for non-square matrices. A matrix is diagonal if
``M[i,j] == 0`` whenever ``i != j``.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DM
>>> M = DM([[ZZ(1), ZZ(0)], [ZZ(0), ZZ(1)]], ZZ)
>>> M.is_diagonal
True
See Also
========
is_upper
is_lower
is_square
diagonal
"""
return self.rep.is_diagonal()
def diagonal(self):
"""
Get the diagonal entries of the matrix as a list.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DM
>>> M = DM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], ZZ)
>>> M.diagonal()
[1, 4]
See Also
========
is_diagonal
diag
"""
return self.rep.diagonal()
@property
def is_square(self):
"""
True if the matrix is square.
"""
return self.shape[0] == self.shape[1]
def rank(self):
rref, pivots = self.rref()
return len(pivots)
def hstack(A, *B):
r"""Horizontally stack the given matrices.
Parameters
==========
B: DomainMatrix
Matrices to stack horizontally.
Returns
=======
DomainMatrix
DomainMatrix by stacking horizontally.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> B = DomainMatrix([[ZZ(5), ZZ(6)], [ZZ(7), ZZ(8)]], (2, 2), ZZ)
>>> A.hstack(B)
DomainMatrix([[1, 2, 5, 6], [3, 4, 7, 8]], (2, 4), ZZ)
>>> C = DomainMatrix([[ZZ(9), ZZ(10)], [ZZ(11), ZZ(12)]], (2, 2), ZZ)
>>> A.hstack(B, C)
DomainMatrix([[1, 2, 5, 6, 9, 10], [3, 4, 7, 8, 11, 12]], (2, 6), ZZ)
See Also
========
unify
"""
A, *B = A.unify(*B, fmt=A.rep.fmt)
return DomainMatrix.from_rep(A.rep.hstack(*(Bk.rep for Bk in B)))
def vstack(A, *B):
r"""Vertically stack the given matrices.
Parameters
==========
B: DomainMatrix
Matrices to stack vertically.
Returns
=======
DomainMatrix
DomainMatrix by stacking vertically.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> B = DomainMatrix([[ZZ(5), ZZ(6)], [ZZ(7), ZZ(8)]], (2, 2), ZZ)
>>> A.vstack(B)
DomainMatrix([[1, 2], [3, 4], [5, 6], [7, 8]], (4, 2), ZZ)
>>> C = DomainMatrix([[ZZ(9), ZZ(10)], [ZZ(11), ZZ(12)]], (2, 2), ZZ)
>>> A.vstack(B, C)
DomainMatrix([[1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12]], (6, 2), ZZ)
See Also
========
unify
"""
A, *B = A.unify(*B, fmt='dense')
return DomainMatrix.from_rep(A.rep.vstack(*(Bk.rep for Bk in B)))
def applyfunc(self, func, domain=None):
if domain is None:
domain = self.domain
return self.from_rep(self.rep.applyfunc(func, domain))
def __add__(A, B):
if not isinstance(B, DomainMatrix):
return NotImplemented
A, B = A.unify(B, fmt='dense')
return A.add(B)
def __sub__(A, B):
if not isinstance(B, DomainMatrix):
return NotImplemented
A, B = A.unify(B, fmt='dense')
return A.sub(B)
def __neg__(A):
return A.neg()
def __mul__(A, B):
"""A * B"""
if isinstance(B, DomainMatrix):
A, B = A.unify(B, fmt='dense')
return A.matmul(B)
elif B in A.domain:
return A.scalarmul(B)
elif isinstance(B, DomainScalar):
A, B = A.unify(B)
return A.scalarmul(B.element)
else:
return NotImplemented
def __rmul__(A, B):
if B in A.domain:
return A.rscalarmul(B)
elif isinstance(B, DomainScalar):
A, B = A.unify(B)
return A.rscalarmul(B.element)
else:
return NotImplemented
def __pow__(A, n):
"""A ** n"""
if not isinstance(n, int):
return NotImplemented
return A.pow(n)
def _check(a, op, b, ashape, bshape):
if a.domain != b.domain:
msg = "Domain mismatch: %s %s %s" % (a.domain, op, b.domain)
raise DMDomainError(msg)
if ashape != bshape:
msg = "Shape mismatch: %s %s %s" % (a.shape, op, b.shape)
raise DMShapeError(msg)
if a.rep.fmt != b.rep.fmt:
msg = "Format mismatch: %s %s %s" % (a.rep.fmt, op, b.rep.fmt)
raise DMFormatError(msg)
if type(a.rep) != type(b.rep):
msg = "Type mismatch: %s %s %s" % (type(a.rep), op, type(b.rep))
raise DMFormatError(msg)
def add(A, B):
r"""
Adds two DomainMatrix matrices of the same Domain
Parameters
==========
A, B: DomainMatrix
matrices to add
Returns
=======
DomainMatrix
DomainMatrix after Addition
Raises
======
DMShapeError
If the dimensions of the two DomainMatrix are not equal
ValueError
If the domain of the two DomainMatrix are not same
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> B = DomainMatrix([
... [ZZ(4), ZZ(3)],
... [ZZ(2), ZZ(1)]], (2, 2), ZZ)
>>> A.add(B)
DomainMatrix([[5, 5], [5, 5]], (2, 2), ZZ)
See Also
========
sub, matmul
"""
A._check('+', B, A.shape, B.shape)
return A.from_rep(A.rep.add(B.rep))
def sub(A, B):
r"""
Subtracts two DomainMatrix matrices of the same Domain
Parameters
==========
A, B: DomainMatrix
matrices to subtract
Returns
=======
DomainMatrix
DomainMatrix after Subtraction
Raises
======
DMShapeError
If the dimensions of the two DomainMatrix are not equal
ValueError
If the domain of the two DomainMatrix are not same
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> B = DomainMatrix([
... [ZZ(4), ZZ(3)],
... [ZZ(2), ZZ(1)]], (2, 2), ZZ)
>>> A.sub(B)
DomainMatrix([[-3, -1], [1, 3]], (2, 2), ZZ)
See Also
========
add, matmul
"""
A._check('-', B, A.shape, B.shape)
return A.from_rep(A.rep.sub(B.rep))
def neg(A):
r"""
Returns the negative of DomainMatrix
Parameters
==========
A : Represents a DomainMatrix
Returns
=======
DomainMatrix
DomainMatrix after Negation
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> A.neg()
DomainMatrix([[-1, -2], [-3, -4]], (2, 2), ZZ)
"""
return A.from_rep(A.rep.neg())
def mul(A, b):
r"""
Performs term by term multiplication for the second DomainMatrix
w.r.t first DomainMatrix. Returns a DomainMatrix whose rows are
list of DomainMatrix matrices created after term by term multiplication.
Parameters
==========
A, B: DomainMatrix
matrices to multiply term-wise
Returns
=======
DomainMatrix
DomainMatrix after term by term multiplication
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> b = ZZ(2)
>>> A.mul(b)
DomainMatrix([[2, 4], [6, 8]], (2, 2), ZZ)
See Also
========
matmul
"""
return A.from_rep(A.rep.mul(b))
def rmul(A, b):
return A.from_rep(A.rep.rmul(b))
def matmul(A, B):
r"""
Performs matrix multiplication of two DomainMatrix matrices
Parameters
==========
A, B: DomainMatrix
to multiply
Returns
=======
DomainMatrix
DomainMatrix after multiplication
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> B = DomainMatrix([
... [ZZ(1), ZZ(1)],
... [ZZ(0), ZZ(1)]], (2, 2), ZZ)
>>> A.matmul(B)
DomainMatrix([[1, 3], [3, 7]], (2, 2), ZZ)
See Also
========
mul, pow, add, sub
"""
A._check('*', B, A.shape[1], B.shape[0])
return A.from_rep(A.rep.matmul(B.rep))
def _scalarmul(A, lamda, reverse):
if lamda == A.domain.zero:
return DomainMatrix.zeros(A.shape, A.domain)
elif lamda == A.domain.one:
return A.copy()
elif reverse:
return A.rmul(lamda)
else:
return A.mul(lamda)
def scalarmul(A, lamda):
return A._scalarmul(lamda, reverse=False)
def rscalarmul(A, lamda):
return A._scalarmul(lamda, reverse=True)
def mul_elementwise(A, B):
assert A.domain == B.domain
return A.from_rep(A.rep.mul_elementwise(B.rep))
def __truediv__(A, lamda):
""" Method for Scalar Division"""
if isinstance(lamda, int) or ZZ.of_type(lamda):
lamda = DomainScalar(ZZ(lamda), ZZ)
elif A.domain.is_Field and lamda in A.domain:
K = A.domain
lamda = DomainScalar(K.convert(lamda), K)
if not isinstance(lamda, DomainScalar):
return NotImplemented
A, lamda = A.to_field().unify(lamda)
if lamda.element == lamda.domain.zero:
raise ZeroDivisionError
if lamda.element == lamda.domain.one:
return A
return A.mul(1 / lamda.element)
def pow(A, n):
r"""
Computes A**n
Parameters
==========
A : DomainMatrix
n : exponent for A
Returns
=======
DomainMatrix
DomainMatrix on computing A**n
Raises
======
NotImplementedError
if n is negative.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(1)],
... [ZZ(0), ZZ(1)]], (2, 2), ZZ)
>>> A.pow(2)
DomainMatrix([[1, 2], [0, 1]], (2, 2), ZZ)
See Also
========
matmul
"""
nrows, ncols = A.shape
if nrows != ncols:
raise DMNonSquareMatrixError('Power of a nonsquare matrix')
if n < 0:
raise NotImplementedError('Negative powers')
elif n == 0:
return A.eye(nrows, A.domain)
elif n == 1:
return A
elif n % 2 == 1:
return A * A**(n - 1)
else:
sqrtAn = A ** (n // 2)
return sqrtAn * sqrtAn
def scc(self):
"""Compute the strongly connected components of a DomainMatrix
Explanation
===========
A square matrix can be considered as the adjacency matrix for a
directed graph where the row and column indices are the vertices. In
this graph if there is an edge from vertex ``i`` to vertex ``j`` if
``M[i, j]`` is nonzero. This routine computes the strongly connected
components of that graph which are subsets of the rows and columns that
are connected by some nonzero element of the matrix. The strongly
connected components are useful because many operations such as the
determinant can be computed by working with the submatrices
corresponding to each component.
Examples
========
Find the strongly connected components of a matrix:
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> M = DomainMatrix([[ZZ(1), ZZ(0), ZZ(2)],
... [ZZ(0), ZZ(3), ZZ(0)],
... [ZZ(4), ZZ(6), ZZ(5)]], (3, 3), ZZ)
>>> M.scc()
[[1], [0, 2]]
Compute the determinant from the components:
>>> MM = M.to_Matrix()
>>> MM
Matrix([
[1, 0, 2],
[0, 3, 0],
[4, 6, 5]])
>>> MM[[1], [1]]
Matrix([[3]])
>>> MM[[0, 2], [0, 2]]
Matrix([
[1, 2],
[4, 5]])
>>> MM.det()
-9
>>> MM[[1], [1]].det() * MM[[0, 2], [0, 2]].det()
-9
The components are given in reverse topological order and represent a
permutation of the rows and columns that will bring the matrix into
block lower-triangular form:
>>> MM[[1, 0, 2], [1, 0, 2]]
Matrix([
[3, 0, 0],
[0, 1, 2],
[6, 4, 5]])
Returns
=======
List of lists of integers
Each list represents a strongly connected component.
See also
========
sympy.matrices.matrixbase.MatrixBase.strongly_connected_components
sympy.utilities.iterables.strongly_connected_components
"""
if not self.is_square:
raise DMNonSquareMatrixError('Matrix must be square for scc')
return self.rep.scc()
def clear_denoms(self, convert=False):
"""
Clear denominators, but keep the domain unchanged.
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DM
>>> A = DM([[(1,2), (1,3)], [(1,4), (1,5)]], QQ)
>>> den, Anum = A.clear_denoms()
>>> den.to_sympy()
60
>>> Anum.to_Matrix()
Matrix([
[30, 20],
[15, 12]])
>>> den * A == Anum
True
The numerator matrix will be in the same domain as the original matrix
unless ``convert`` is set to ``True``:
>>> A.clear_denoms()[1].domain
QQ
>>> A.clear_denoms(convert=True)[1].domain
ZZ
The denominator is always in the associated ring:
>>> A.clear_denoms()[0].domain
ZZ
>>> A.domain.get_ring()
ZZ
See Also
========
sympy.polys.polytools.Poly.clear_denoms
clear_denoms_rowwise
"""
elems0, data = self.to_flat_nz()
K0 = self.domain
K1 = K0.get_ring() if K0.has_assoc_Ring else K0
den, elems1 = dup_clear_denoms(elems0, K0, K1, convert=convert)
if convert:
Kden, Knum = K1, K1
else:
Kden, Knum = K1, K0
den = DomainScalar(den, Kden)
num = self.from_flat_nz(elems1, data, Knum)
return den, num
def clear_denoms_rowwise(self, convert=False):
"""
Clear denominators from each row of the matrix.
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DM
>>> A = DM([[(1,2), (1,3), (1,4)], [(1,5), (1,6), (1,7)]], QQ)
>>> den, Anum = A.clear_denoms_rowwise()
>>> den.to_Matrix()
Matrix([
[12, 0],
[ 0, 210]])
>>> Anum.to_Matrix()
Matrix([
[ 6, 4, 3],
[42, 35, 30]])
The denominator matrix is a diagonal matrix with the denominators of
each row on the diagonal. The invariants are:
>>> den * A == Anum
True
>>> A == den.to_field().inv() * Anum
True
The numerator matrix will be in the same domain as the original matrix
unless ``convert`` is set to ``True``:
>>> A.clear_denoms_rowwise()[1].domain
QQ
>>> A.clear_denoms_rowwise(convert=True)[1].domain
ZZ
The domain of the denominator matrix is the associated ring:
>>> A.clear_denoms_rowwise()[0].domain
ZZ
See Also
========
sympy.polys.polytools.Poly.clear_denoms
clear_denoms
"""
dod = self.to_dod()
K0 = self.domain
K1 = K0.get_ring() if K0.has_assoc_Ring else K0
diagonals = [K0.one] * self.shape[0]
dod_num = {}
for i, rowi in dod.items():
indices, elems = zip(*rowi.items())
den, elems_num = dup_clear_denoms(elems, K0, K1, convert=convert)
rowi_num = dict(zip(indices, elems_num))
diagonals[i] = den
dod_num[i] = rowi_num
if convert:
Kden, Knum = K1, K1
else:
Kden, Knum = K1, K0
den = self.diag(diagonals, Kden)
num = self.from_dod_like(dod_num, Knum)
return den, num
def cancel_denom(self, denom):
"""
Cancel factors between a matrix and a denominator.
Returns a matrix and denominator on lowest terms.
Requires ``gcd`` in the ground domain.
Methods like :meth:`solve_den`, :meth:`inv_den` and :meth:`rref_den`
return a matrix and denominator but not necessarily on lowest terms.
Reduction to lowest terms without fractions can be performed with
:meth:`cancel_denom`.
Examples
========
>>> from sympy.polys.matrices import DM
>>> from sympy import ZZ
>>> M = DM([[2, 2, 0],
... [0, 2, 2],
... [0, 0, 2]], ZZ)
>>> Minv, den = M.inv_den()
>>> Minv.to_Matrix()
Matrix([
[1, -1, 1],
[0, 1, -1],
[0, 0, 1]])
>>> den
2
>>> Minv_reduced, den_reduced = Minv.cancel_denom(den)
>>> Minv_reduced.to_Matrix()
Matrix([
[1, -1, 1],
[0, 1, -1],
[0, 0, 1]])
>>> den_reduced
2
>>> Minv_reduced.to_field() / den_reduced == Minv.to_field() / den
True
The denominator is made canonical with respect to units (e.g. a
negative denominator is made positive):
>>> M = DM([[2, 2, 0]], ZZ)
>>> den = ZZ(-4)
>>> M.cancel_denom(den)
(DomainMatrix([[-1, -1, 0]], (1, 3), ZZ), 2)
Any factor common to _all_ elements will be cancelled but there can
still be factors in common between _some_ elements of the matrix and
the denominator. To cancel factors between each element and the
denominator, use :meth:`cancel_denom_elementwise` or otherwise convert
to a field and use division:
>>> M = DM([[4, 6]], ZZ)
>>> den = ZZ(12)
>>> M.cancel_denom(den)
(DomainMatrix([[2, 3]], (1, 2), ZZ), 6)
>>> numers, denoms = M.cancel_denom_elementwise(den)
>>> numers
DomainMatrix([[1, 1]], (1, 2), ZZ)
>>> denoms
DomainMatrix([[3, 2]], (1, 2), ZZ)
>>> M.to_field() / den
DomainMatrix([[1/3, 1/2]], (1, 2), QQ)
See Also
========
solve_den
inv_den
rref_den
cancel_denom_elementwise
"""
M = self
K = self.domain
if K.is_zero(denom):
raise ZeroDivisionError('denominator is zero')
elif K.is_one(denom):
return (M.copy(), denom)
elements, data = M.to_flat_nz()
# First canonicalize the denominator (e.g. multiply by -1).
if K.is_negative(denom):
u = -K.one
else:
u = K.canonical_unit(denom)
# Often after e.g. solve_den the denominator will be much more
# complicated than the elements of the numerator. Hopefully it will be
# quicker to find the gcd of the numerator and if there is no content
# then we do not need to look at the denominator at all.
content = dup_content(elements, K)
common = K.gcd(content, denom)
if not K.is_one(content):
common = K.gcd(content, denom)
if not K.is_one(common):
elements = dup_quo_ground(elements, common, K)
denom = K.quo(denom, common)
if not K.is_one(u):
elements = dup_mul_ground(elements, u, K)
denom = u * denom
elif K.is_one(common):
return (M.copy(), denom)
M_cancelled = M.from_flat_nz(elements, data, K)
return M_cancelled, denom
def cancel_denom_elementwise(self, denom):
"""
Cancel factors between the elements of a matrix and a denominator.
Returns a matrix of numerators and matrix of denominators.
Requires ``gcd`` in the ground domain.
Examples
========
>>> from sympy.polys.matrices import DM
>>> from sympy import ZZ
>>> M = DM([[2, 3], [4, 12]], ZZ)
>>> denom = ZZ(6)
>>> numers, denoms = M.cancel_denom_elementwise(denom)
>>> numers.to_Matrix()
Matrix([
[1, 1],
[2, 2]])
>>> denoms.to_Matrix()
Matrix([
[3, 2],
[3, 1]])
>>> M_frac = (M.to_field() / denom).to_Matrix()
>>> M_frac
Matrix([
[1/3, 1/2],
[2/3, 2]])
>>> denoms_inverted = denoms.to_Matrix().applyfunc(lambda e: 1/e)
>>> numers.to_Matrix().multiply_elementwise(denoms_inverted) == M_frac
True
Use :meth:`cancel_denom` to cancel factors between the matrix and the
denominator while preserving the form of a matrix with a scalar
denominator.
See Also
========
cancel_denom
"""
K = self.domain
M = self
if K.is_zero(denom):
raise ZeroDivisionError('denominator is zero')
elif K.is_one(denom):
M_numers = M.copy()
M_denoms = M.ones(M.shape, M.domain)
return (M_numers, M_denoms)
elements, data = M.to_flat_nz()
cofactors = [K.cofactors(numer, denom) for numer in elements]
gcds, numers, denoms = zip(*cofactors)
M_numers = M.from_flat_nz(list(numers), data, K)
M_denoms = M.from_flat_nz(list(denoms), data, K)
return (M_numers, M_denoms)
def content(self):
"""
Return the gcd of the elements of the matrix.
Requires ``gcd`` in the ground domain.
Examples
========
>>> from sympy.polys.matrices import DM
>>> from sympy import ZZ
>>> M = DM([[2, 4], [4, 12]], ZZ)
>>> M.content()
2
See Also
========
primitive
cancel_denom
"""
K = self.domain
elements, _ = self.to_flat_nz()
return dup_content(elements, K)
def primitive(self):
"""
Factor out gcd of the elements of a matrix.
Requires ``gcd`` in the ground domain.
Examples
========
>>> from sympy.polys.matrices import DM
>>> from sympy import ZZ
>>> M = DM([[2, 4], [4, 12]], ZZ)
>>> content, M_primitive = M.primitive()
>>> content
2
>>> M_primitive
DomainMatrix([[1, 2], [2, 6]], (2, 2), ZZ)
>>> content * M_primitive == M
True
>>> M_primitive.content() == ZZ(1)
True
See Also
========
content
cancel_denom
"""
K = self.domain
elements, data = self.to_flat_nz()
content, prims = dup_primitive(elements, K)
M_primitive = self.from_flat_nz(prims, data, K)
return content, M_primitive
def rref(self, *, method='auto'):
r"""
Returns reduced-row echelon form (RREF) and list of pivots.
If the domain is not a field then it will be converted to a field. See
:meth:`rref_den` for the fraction-free version of this routine that
returns RREF with denominator instead.
The domain must either be a field or have an associated fraction field
(see :meth:`to_field`).
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [QQ(2), QQ(-1), QQ(0)],
... [QQ(-1), QQ(2), QQ(-1)],
... [QQ(0), QQ(0), QQ(2)]], (3, 3), QQ)
>>> rref_matrix, rref_pivots = A.rref()
>>> rref_matrix
DomainMatrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]], (3, 3), QQ)
>>> rref_pivots
(0, 1, 2)
Parameters
==========
method : str, optional (default: 'auto')
The method to use to compute the RREF. The default is ``'auto'``,
which will attempt to choose the fastest method. The other options
are:
- ``A.rref(method='GJ')`` uses Gauss-Jordan elimination with
division. If the domain is not a field then it will be converted
to a field with :meth:`to_field` first and RREF will be computed
by inverting the pivot elements in each row. This is most
efficient for very sparse matrices or for matrices whose elements
have complex denominators.
- ``A.rref(method='FF')`` uses fraction-free Gauss-Jordan
elimination. Elimination is performed using exact division
(``exquo``) to control the growth of the coefficients. In this
case the current domain is always used for elimination but if
the domain is not a field then it will be converted to a field
at the end and divided by the denominator. This is most efficient
for dense matrices or for matrices with simple denominators.
- ``A.rref(method='CD')`` clears the denominators before using
fraction-free Gauss-Jordan elimination in the associated ring.
This is most efficient for dense matrices with very simple
denominators.
- ``A.rref(method='GJ_dense')``, ``A.rref(method='FF_dense')``, and
``A.rref(method='CD_dense')`` are the same as the above methods
except that the dense implementations of the algorithms are used.
By default ``A.rref(method='auto')`` will usually choose the
sparse implementations for RREF.
Regardless of which algorithm is used the returned matrix will
always have the same format (sparse or dense) as the input and its
domain will always be the field of fractions of the input domain.
Returns
=======
(DomainMatrix, list)
reduced-row echelon form and list of pivots for the DomainMatrix
See Also
========
rref_den
RREF with denominator
sympy.polys.matrices.sdm.sdm_irref
Sparse implementation of ``method='GJ'``.
sympy.polys.matrices.sdm.sdm_rref_den
Sparse implementation of ``method='FF'`` and ``method='CD'``.
sympy.polys.matrices.dense.ddm_irref
Dense implementation of ``method='GJ'``.
sympy.polys.matrices.dense.ddm_irref_den
Dense implementation of ``method='FF'`` and ``method='CD'``.
clear_denoms
Clear denominators from a matrix, used by ``method='CD'`` and
by ``method='GJ'`` when the original domain is not a field.
"""
return _dm_rref(self, method=method)
def rref_den(self, *, method='auto', keep_domain=True):
r"""
Returns reduced-row echelon form with denominator and list of pivots.
Requires exact division in the ground domain (``exquo``).
Examples
========
>>> from sympy import ZZ, QQ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(2), ZZ(-1), ZZ(0)],
... [ZZ(-1), ZZ(2), ZZ(-1)],
... [ZZ(0), ZZ(0), ZZ(2)]], (3, 3), ZZ)
>>> A_rref, denom, pivots = A.rref_den()
>>> A_rref
DomainMatrix([[6, 0, 0], [0, 6, 0], [0, 0, 6]], (3, 3), ZZ)
>>> denom
6
>>> pivots
(0, 1, 2)
>>> A_rref.to_field() / denom
DomainMatrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]], (3, 3), QQ)
>>> A_rref.to_field() / denom == A.convert_to(QQ).rref()[0]
True
Parameters
==========
method : str, optional (default: 'auto')
The method to use to compute the RREF. The default is ``'auto'``,
which will attempt to choose the fastest method. The other options
are:
- ``A.rref(method='FF')`` uses fraction-free Gauss-Jordan
elimination. Elimination is performed using exact division
(``exquo``) to control the growth of the coefficients. In this
case the current domain is always used for elimination and the
result is always returned as a matrix over the current domain.
This is most efficient for dense matrices or for matrices with
simple denominators.
- ``A.rref(method='CD')`` clears denominators before using
fraction-free Gauss-Jordan elimination in the associated ring.
The result will be converted back to the original domain unless
``keep_domain=False`` is passed in which case the result will be
over the ring used for elimination. This is most efficient for
dense matrices with very simple denominators.
- ``A.rref(method='GJ')`` uses Gauss-Jordan elimination with
division. If the domain is not a field then it will be converted
to a field with :meth:`to_field` first and RREF will be computed
by inverting the pivot elements in each row. The result is
converted back to the original domain by clearing denominators
unless ``keep_domain=False`` is passed in which case the result
will be over the field used for elimination. This is most
efficient for very sparse matrices or for matrices whose elements
have complex denominators.
- ``A.rref(method='GJ_dense')``, ``A.rref(method='FF_dense')``, and
``A.rref(method='CD_dense')`` are the same as the above methods
except that the dense implementations of the algorithms are used.
By default ``A.rref(method='auto')`` will usually choose the
sparse implementations for RREF.
Regardless of which algorithm is used the returned matrix will
always have the same format (sparse or dense) as the input and if
``keep_domain=True`` its domain will always be the same as the
input.
keep_domain : bool, optional
If True (the default), the domain of the returned matrix and
denominator are the same as the domain of the input matrix. If
False, the domain of the returned matrix might be changed to an
associated ring or field if the algorithm used a different domain.
This is useful for efficiency if the caller does not need the
result to be in the original domain e.g. it avoids clearing
denominators in the case of ``A.rref(method='GJ')``.
Returns
=======
(DomainMatrix, scalar, list)
Reduced-row echelon form, denominator and list of pivot indices.
See Also
========
rref
RREF without denominator for field domains.
sympy.polys.matrices.sdm.sdm_irref
Sparse implementation of ``method='GJ'``.
sympy.polys.matrices.sdm.sdm_rref_den
Sparse implementation of ``method='FF'`` and ``method='CD'``.
sympy.polys.matrices.dense.ddm_irref
Dense implementation of ``method='GJ'``.
sympy.polys.matrices.dense.ddm_irref_den
Dense implementation of ``method='FF'`` and ``method='CD'``.
clear_denoms
Clear denominators from a matrix, used by ``method='CD'``.
"""
return _dm_rref_den(self, method=method, keep_domain=keep_domain)
def columnspace(self):
r"""
Returns the columnspace for the DomainMatrix
Returns
=======
DomainMatrix
The columns of this matrix form a basis for the columnspace.
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [QQ(1), QQ(-1)],
... [QQ(2), QQ(-2)]], (2, 2), QQ)
>>> A.columnspace()
DomainMatrix([[1], [2]], (2, 1), QQ)
"""
if not self.domain.is_Field:
raise DMNotAField('Not a field')
rref, pivots = self.rref()
rows, cols = self.shape
return self.extract(range(rows), pivots)
def rowspace(self):
r"""
Returns the rowspace for the DomainMatrix
Returns
=======
DomainMatrix
The rows of this matrix form a basis for the rowspace.
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [QQ(1), QQ(-1)],
... [QQ(2), QQ(-2)]], (2, 2), QQ)
>>> A.rowspace()
DomainMatrix([[1, -1]], (1, 2), QQ)
"""
if not self.domain.is_Field:
raise DMNotAField('Not a field')
rref, pivots = self.rref()
rows, cols = self.shape
return self.extract(range(len(pivots)), range(cols))
def nullspace(self, divide_last=False):
r"""
Returns the nullspace for the DomainMatrix
Returns
=======
DomainMatrix
The rows of this matrix form a basis for the nullspace.
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DM
>>> A = DM([
... [QQ(2), QQ(-2)],
... [QQ(4), QQ(-4)]], QQ)
>>> A.nullspace()
DomainMatrix([[1, 1]], (1, 2), QQ)
The returned matrix is a basis for the nullspace:
>>> A_null = A.nullspace().transpose()
>>> A * A_null
DomainMatrix([[0], [0]], (2, 1), QQ)
>>> rows, cols = A.shape
>>> nullity = rows - A.rank()
>>> A_null.shape == (cols, nullity)
True
Nullspace can also be computed for non-field rings. If the ring is not
a field then division is not used. Setting ``divide_last`` to True will
raise an error in this case:
>>> from sympy import ZZ
>>> B = DM([[6, -3],
... [4, -2]], ZZ)
>>> B.nullspace()
DomainMatrix([[3, 6]], (1, 2), ZZ)
>>> B.nullspace(divide_last=True)
Traceback (most recent call last):
...
DMNotAField: Cannot normalize vectors over a non-field
Over a ring with ``gcd`` defined the nullspace can potentially be
reduced with :meth:`primitive`:
>>> B.nullspace().primitive()
(3, DomainMatrix([[1, 2]], (1, 2), ZZ))
A matrix over a ring can often be normalized by converting it to a
field but it is often a bad idea to do so:
>>> from sympy.abc import a, b, c
>>> from sympy import Matrix
>>> M = Matrix([[ a*b, b + c, c],
... [ a - b, b*c, c**2],
... [a*b + a - b, b*c + b + c, c**2 + c]])
>>> M.to_DM().domain
ZZ[a,b,c]
>>> M.to_DM().nullspace().to_Matrix().transpose()
Matrix([
[ c**3],
[ -a*b*c**2 + a*c - b*c],
[a*b**2*c - a*b - a*c + b**2 + b*c]])
The unnormalized form here is nicer than the normalized form that
spreads a large denominator throughout the matrix:
>>> M.to_DM().to_field().nullspace(divide_last=True).to_Matrix().transpose()
Matrix([
[ c**3/(a*b**2*c - a*b - a*c + b**2 + b*c)],
[(-a*b*c**2 + a*c - b*c)/(a*b**2*c - a*b - a*c + b**2 + b*c)],
[ 1]])
Parameters
==========
divide_last : bool, optional
If False (the default), the vectors are not normalized and the RREF
is computed using :meth:`rref_den` and the denominator is
discarded. If True, then each row is divided by its final element;
the domain must be a field in this case.
See Also
========
nullspace_from_rref
rref
rref_den
rowspace
"""
A = self
K = A.domain
if divide_last and not K.is_Field:
raise DMNotAField("Cannot normalize vectors over a non-field")
if divide_last:
A_rref, pivots = A.rref()
else:
A_rref, den, pivots = A.rref_den()
# Ensure that the sign is canonical before discarding the
# denominator. Then M.nullspace().primitive() is canonical.
u = K.canonical_unit(den)
if u != K.one:
A_rref *= u
A_null = A_rref.nullspace_from_rref(pivots)
return A_null
def nullspace_from_rref(self, pivots=None):
"""
Compute nullspace from rref and pivots.
The domain of the matrix can be any domain.
The matrix must be in reduced row echelon form already. Otherwise the
result will be incorrect. Use :meth:`rref` or :meth:`rref_den` first
to get the reduced row echelon form or use :meth:`nullspace` instead.
See Also
========
nullspace
rref
rref_den
sympy.polys.matrices.sdm.SDM.nullspace_from_rref
sympy.polys.matrices.ddm.DDM.nullspace_from_rref
"""
null_rep, nonpivots = self.rep.nullspace_from_rref(pivots)
return self.from_rep(null_rep)
def inv(self):
r"""
Finds the inverse of the DomainMatrix if exists
Returns
=======
DomainMatrix
DomainMatrix after inverse
Raises
======
ValueError
If the domain of DomainMatrix not a Field
DMNonSquareMatrixError
If the DomainMatrix is not a not Square DomainMatrix
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [QQ(2), QQ(-1), QQ(0)],
... [QQ(-1), QQ(2), QQ(-1)],
... [QQ(0), QQ(0), QQ(2)]], (3, 3), QQ)
>>> A.inv()
DomainMatrix([[2/3, 1/3, 1/6], [1/3, 2/3, 1/3], [0, 0, 1/2]], (3, 3), QQ)
See Also
========
neg
"""
if not self.domain.is_Field:
raise DMNotAField('Not a field')
m, n = self.shape
if m != n:
raise DMNonSquareMatrixError
inv = self.rep.inv()
return self.from_rep(inv)
def det(self):
r"""
Returns the determinant of a square :class:`DomainMatrix`.
Returns
=======
determinant: DomainElement
Determinant of the matrix.
Raises
======
ValueError
If the domain of DomainMatrix is not a Field
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> A.det()
-2
"""
m, n = self.shape
if m != n:
raise DMNonSquareMatrixError
return self.rep.det()
def adj_det(self):
"""
Adjugate and determinant of a square :class:`DomainMatrix`.
Returns
=======
(adjugate, determinant) : (DomainMatrix, DomainScalar)
The adjugate matrix and determinant of this matrix.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DM
>>> A = DM([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], ZZ)
>>> adjA, detA = A.adj_det()
>>> adjA
DomainMatrix([[4, -2], [-3, 1]], (2, 2), ZZ)
>>> detA
-2
See Also
========
adjugate
Returns only the adjugate matrix.
det
Returns only the determinant.
inv_den
Returns a matrix/denominator pair representing the inverse matrix
but perhaps differing from the adjugate and determinant by a common
factor.
"""
m, n = self.shape
I_m = self.eye((m, m), self.domain)
adjA, detA = self.solve_den_charpoly(I_m, check=False)
if self.rep.fmt == "dense":
adjA = adjA.to_dense()
return adjA, detA
def adjugate(self):
"""
Adjugate of a square :class:`DomainMatrix`.
The adjugate matrix is the transpose of the cofactor matrix and is
related to the inverse by::
adj(A) = det(A) * A.inv()
Unlike the inverse matrix the adjugate matrix can be computed and
expressed without division or fractions in the ground domain.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DM
>>> A = DM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], ZZ)
>>> A.adjugate()
DomainMatrix([[4, -2], [-3, 1]], (2, 2), ZZ)
Returns
=======
DomainMatrix
The adjugate matrix of this matrix with the same domain.
See Also
========
adj_det
"""
adjA, detA = self.adj_det()
return adjA
def inv_den(self, method=None):
"""
Return the inverse as a :class:`DomainMatrix` with denominator.
Returns
=======
(inv, den) : (:class:`DomainMatrix`, :class:`~.DomainElement`)
The inverse matrix and its denominator.
This is more or less equivalent to :meth:`adj_det` except that ``inv``
and ``den`` are not guaranteed to be the adjugate and inverse. The
ratio ``inv/den`` is equivalent to ``adj/det`` but some factors
might be cancelled between ``inv`` and ``den``. In simple cases this
might just be a minus sign so that ``(inv, den) == (-adj, -det)`` but
factors more complicated than ``-1`` can also be cancelled.
Cancellation is not guaranteed to be complete so ``inv`` and ``den``
may not be on lowest terms. The denominator ``den`` will be zero if and
only if the determinant is zero.
If the actual adjugate and determinant are needed, use :meth:`adj_det`
instead. If the intention is to compute the inverse matrix or solve a
system of equations then :meth:`inv_den` is more efficient.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(2), ZZ(-1), ZZ(0)],
... [ZZ(-1), ZZ(2), ZZ(-1)],
... [ZZ(0), ZZ(0), ZZ(2)]], (3, 3), ZZ)
>>> Ainv, den = A.inv_den()
>>> den
6
>>> Ainv
DomainMatrix([[4, 2, 1], [2, 4, 2], [0, 0, 3]], (3, 3), ZZ)
>>> A * Ainv == den * A.eye(A.shape, A.domain).to_dense()
True
Parameters
==========
method : str, optional
The method to use to compute the inverse. Can be one of ``None``,
``'rref'`` or ``'charpoly'``. If ``None`` then the method is
chosen automatically (see :meth:`solve_den` for details).
See Also
========
inv
det
adj_det
solve_den
"""
I = self.eye(self.shape, self.domain)
return self.solve_den(I, method=method)
def solve_den(self, b, method=None):
"""
Solve matrix equation $Ax = b$ without fractions in the ground domain.
Examples
========
Solve a matrix equation over the integers:
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DM
>>> A = DM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], ZZ)
>>> b = DM([[ZZ(5)], [ZZ(6)]], ZZ)
>>> xnum, xden = A.solve_den(b)
>>> xden
-2
>>> xnum
DomainMatrix([[8], [-9]], (2, 1), ZZ)
>>> A * xnum == xden * b
True
Solve a matrix equation over a polynomial ring:
>>> from sympy import ZZ
>>> from sympy.abc import x, y, z, a, b
>>> R = ZZ[x, y, z, a, b]
>>> M = DM([[x*y, x*z], [y*z, x*z]], R)
>>> b = DM([[a], [b]], R)
>>> M.to_Matrix()
Matrix([
[x*y, x*z],
[y*z, x*z]])
>>> b.to_Matrix()
Matrix([
[a],
[b]])
>>> xnum, xden = M.solve_den(b)
>>> xden
x**2*y*z - x*y*z**2
>>> xnum.to_Matrix()
Matrix([
[ a*x*z - b*x*z],
[-a*y*z + b*x*y]])
>>> M * xnum == xden * b
True
The solution can be expressed over a fraction field which will cancel
gcds between the denominator and the elements of the numerator:
>>> xsol = xnum.to_field() / xden
>>> xsol.to_Matrix()
Matrix([
[ (a - b)/(x*y - y*z)],
[(-a*z + b*x)/(x**2*z - x*z**2)]])
>>> (M * xsol).to_Matrix() == b.to_Matrix()
True
When solving a large system of equations this cancellation step might
be a lot slower than :func:`solve_den` itself. The solution can also be
expressed as a ``Matrix`` without attempting any polynomial
cancellation between the numerator and denominator giving a less
simplified result more quickly:
>>> xsol_uncancelled = xnum.to_Matrix() / xnum.domain.to_sympy(xden)
>>> xsol_uncancelled
Matrix([
[ (a*x*z - b*x*z)/(x**2*y*z - x*y*z**2)],
[(-a*y*z + b*x*y)/(x**2*y*z - x*y*z**2)]])
>>> from sympy import cancel
>>> cancel(xsol_uncancelled) == xsol.to_Matrix()
True
Parameters
==========
self : :class:`DomainMatrix`
The ``m x n`` matrix $A$ in the equation $Ax = b$. Underdetermined
systems are not supported so ``m >= n``: $A$ should be square or
have more rows than columns.
b : :class:`DomainMatrix`
The ``n x m`` matrix $b$ for the rhs.
cp : list of :class:`~.DomainElement`, optional
The characteristic polynomial of the matrix $A$. If not given, it
will be computed using :meth:`charpoly`.
method: str, optional
The method to use for solving the system. Can be one of ``None``,
``'charpoly'`` or ``'rref'``. If ``None`` (the default) then the
method will be chosen automatically.
The ``charpoly`` method uses :meth:`solve_den_charpoly` and can
only be used if the matrix is square. This method is division free
and can be used with any domain.
The ``rref`` method is fraction free but requires exact division
in the ground domain (``exquo``). This is also suitable for most
domains. This method can be used with overdetermined systems (more
equations than unknowns) but not underdetermined systems as a
unique solution is sought.
Returns
=======
(xnum, xden) : (DomainMatrix, DomainElement)
The solution of the equation $Ax = b$ as a pair consisting of an
``n x m`` matrix numerator ``xnum`` and a scalar denominator
``xden``.
The solution $x$ is given by ``x = xnum / xden``. The division free
invariant is ``A * xnum == xden * b``. If $A$ is square then the
denominator ``xden`` will be a divisor of the determinant $det(A)$.
Raises
======
DMNonInvertibleMatrixError
If the system $Ax = b$ does not have a unique solution.
See Also
========
solve_den_charpoly
solve_den_rref
inv_den
"""
m, n = self.shape
bm, bn = b.shape
if m != bm:
raise DMShapeError("Matrix equation shape mismatch.")
if method is None:
method = 'rref'
elif method == 'charpoly' and m != n:
raise DMNonSquareMatrixError("method='charpoly' requires a square matrix.")
if method == 'charpoly':
xnum, xden = self.solve_den_charpoly(b)
elif method == 'rref':
xnum, xden = self.solve_den_rref(b)
else:
raise DMBadInputError("method should be 'rref' or 'charpoly'")
return xnum, xden
def solve_den_rref(self, b):
"""
Solve matrix equation $Ax = b$ using fraction-free RREF
Solves the matrix equation $Ax = b$ for $x$ and returns the solution
as a numerator/denominator pair.
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DM
>>> A = DM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], ZZ)
>>> b = DM([[ZZ(5)], [ZZ(6)]], ZZ)
>>> xnum, xden = A.solve_den_rref(b)
>>> xden
-2
>>> xnum
DomainMatrix([[8], [-9]], (2, 1), ZZ)
>>> A * xnum == xden * b
True
See Also
========
solve_den
solve_den_charpoly
"""
A = self
m, n = A.shape
bm, bn = b.shape
if m != bm:
raise DMShapeError("Matrix equation shape mismatch.")
if m < n:
raise DMShapeError("Underdetermined matrix equation.")
Aaug = A.hstack(b)
Aaug_rref, denom, pivots = Aaug.rref_den()
# XXX: We check here if there are pivots after the last column. If
# there were than it possibly means that rref_den performed some
# unnecessary elimination. It would be better if rref methods had a
# parameter indicating how many columns should be used for elimination.
if len(pivots) != n or pivots and pivots[-1] >= n:
raise DMNonInvertibleMatrixError("Non-unique solution.")
xnum = Aaug_rref[:n, n:]
xden = denom
return xnum, xden
def solve_den_charpoly(self, b, cp=None, check=True):
"""
Solve matrix equation $Ax = b$ using the characteristic polynomial.
This method solves the square matrix equation $Ax = b$ for $x$ using
the characteristic polynomial without any division or fractions in the
ground domain.
Examples
========
Solve a matrix equation over the integers:
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DM
>>> A = DM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], ZZ)
>>> b = DM([[ZZ(5)], [ZZ(6)]], ZZ)
>>> xnum, detA = A.solve_den_charpoly(b)
>>> detA
-2
>>> xnum
DomainMatrix([[8], [-9]], (2, 1), ZZ)
>>> A * xnum == detA * b
True
Parameters
==========
self : DomainMatrix
The ``n x n`` matrix `A` in the equation `Ax = b`. Must be square
and invertible.
b : DomainMatrix
The ``n x m`` matrix `b` for the rhs.
cp : list, optional
The characteristic polynomial of the matrix `A` if known. If not
given, it will be computed using :meth:`charpoly`.
check : bool, optional
If ``True`` (the default) check that the determinant is not zero
and raise an error if it is. If ``False`` then if the determinant
is zero the return value will be equal to ``(A.adjugate()*b, 0)``.
Returns
=======
(xnum, detA) : (DomainMatrix, DomainElement)
The solution of the equation `Ax = b` as a matrix numerator and
scalar denominator pair. The denominator is equal to the
determinant of `A` and the numerator is ``adj(A)*b``.
The solution $x$ is given by ``x = xnum / detA``. The division free
invariant is ``A * xnum == detA * b``.
If ``b`` is the identity matrix, then ``xnum`` is the adjugate matrix
and we have ``A * adj(A) == detA * I``.
See Also
========
solve_den
Main frontend for solving matrix equations with denominator.
solve_den_rref
Solve matrix equations using fraction-free RREF.
inv_den
Invert a matrix using the characteristic polynomial.
"""
A, b = self.unify(b)
m, n = self.shape
mb, nb = b.shape
if m != n:
raise DMNonSquareMatrixError("Matrix must be square")
if mb != m:
raise DMShapeError("Matrix and vector must have the same number of rows")
f, detA = self.adj_poly_det(cp=cp)
if check and not detA:
raise DMNonInvertibleMatrixError("Matrix is not invertible")
# Compute adj(A)*b = det(A)*inv(A)*b using Horner's method without
# constructing inv(A) explicitly.
adjA_b = self.eval_poly_mul(f, b)
return (adjA_b, detA)
def adj_poly_det(self, cp=None):
"""
Return the polynomial $p$ such that $p(A) = adj(A)$ and also the
determinant of $A$.
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DM
>>> A = DM([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], QQ)
>>> p, detA = A.adj_poly_det()
>>> p
[-1, 5]
>>> p_A = A.eval_poly(p)
>>> p_A
DomainMatrix([[4, -2], [-3, 1]], (2, 2), QQ)
>>> p[0]*A**1 + p[1]*A**0 == p_A
True
>>> p_A == A.adjugate()
True
>>> A * A.adjugate() == detA * A.eye(A.shape, A.domain).to_dense()
True
See Also
========
adjugate
eval_poly
adj_det
"""
# Cayley-Hamilton says that a matrix satisfies its own minimal
# polynomial
#
# p[0]*A^n + p[1]*A^(n-1) + ... + p[n]*I = 0
#
# with p[0]=1 and p[n]=(-1)^n*det(A) or
#
# det(A)*I = -(-1)^n*(p[0]*A^(n-1) + p[1]*A^(n-2) + ... + p[n-1]*A).
#
# Define a new polynomial f with f[i] = -(-1)^n*p[i] for i=0..n-1. Then
#
# det(A)*I = f[0]*A^n + f[1]*A^(n-1) + ... + f[n-1]*A.
#
# Multiplying on the right by inv(A) gives
#
# det(A)*inv(A) = f[0]*A^(n-1) + f[1]*A^(n-2) + ... + f[n-1].
#
# So adj(A) = det(A)*inv(A) = f(A)
A = self
m, n = self.shape
if m != n:
raise DMNonSquareMatrixError("Matrix must be square")
if cp is None:
cp = A.charpoly()
if len(cp) % 2:
# n is even
detA = cp[-1]
f = [-cpi for cpi in cp[:-1]]
else:
# n is odd
detA = -cp[-1]
f = cp[:-1]
return f, detA
def eval_poly(self, p):
"""
Evaluate polynomial function of a matrix $p(A)$.
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DM
>>> A = DM([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], QQ)
>>> p = [QQ(1), QQ(2), QQ(3)]
>>> p_A = A.eval_poly(p)
>>> p_A
DomainMatrix([[12, 14], [21, 33]], (2, 2), QQ)
>>> p_A == p[0]*A**2 + p[1]*A + p[2]*A**0
True
See Also
========
eval_poly_mul
"""
A = self
m, n = A.shape
if m != n:
raise DMNonSquareMatrixError("Matrix must be square")
if not p:
return self.zeros(self.shape, self.domain)
elif len(p) == 1:
return p[0] * self.eye(self.shape, self.domain)
# Evaluate p(A) using Horner's method:
# XXX: Use Paterson-Stockmeyer method?
I = A.eye(A.shape, A.domain)
p_A = p[0] * I
for pi in p[1:]:
p_A = A*p_A + pi*I
return p_A
def eval_poly_mul(self, p, B):
r"""
Evaluate polynomial matrix product $p(A) \times B$.
Evaluate the polynomial matrix product $p(A) \times B$ using Horner's
method without creating the matrix $p(A)$ explicitly. If $B$ is a
column matrix then this method will only use matrix-vector multiplies
and no matrix-matrix multiplies are needed.
If $B$ is square or wide or if $A$ can be represented in a simpler
domain than $B$ then it might be faster to evaluate $p(A)$ explicitly
(see :func:`eval_poly`) and then multiply with $B$.
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DM
>>> A = DM([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], QQ)
>>> b = DM([[QQ(5)], [QQ(6)]], QQ)
>>> p = [QQ(1), QQ(2), QQ(3)]
>>> p_A_b = A.eval_poly_mul(p, b)
>>> p_A_b
DomainMatrix([[144], [303]], (2, 1), QQ)
>>> p_A_b == p[0]*A**2*b + p[1]*A*b + p[2]*b
True
>>> A.eval_poly_mul(p, b) == A.eval_poly(p)*b
True
See Also
========
eval_poly
solve_den_charpoly
"""
A = self
m, n = A.shape
mb, nb = B.shape
if m != n:
raise DMNonSquareMatrixError("Matrix must be square")
if mb != n:
raise DMShapeError("Matrices are not aligned")
if A.domain != B.domain:
raise DMDomainError("Matrices must have the same domain")
# Given a polynomial p(x) = p[0]*x^n + p[1]*x^(n-1) + ... + p[n-1]
# and matrices A and B we want to find
#
# p(A)*B = p[0]*A^n*B + p[1]*A^(n-1)*B + ... + p[n-1]*B
#
# Factoring out A term by term we get
#
# p(A)*B = A*(...A*(A*(A*(p[0]*B) + p[1]*B) + p[2]*B) + ...) + p[n-1]*B
#
# where each pair of brackets represents one iteration of the loop
# below starting from the innermost p[0]*B. If B is a column matrix
# then products like A*(...) are matrix-vector multiplies and products
# like p[i]*B are scalar-vector multiplies so there are no
# matrix-matrix multiplies.
if not p:
return B.zeros(B.shape, B.domain, fmt=B.rep.fmt)
p_A_B = p[0]*B
for p_i in p[1:]:
p_A_B = A*p_A_B + p_i*B
return p_A_B
def lu(self):
r"""
Returns Lower and Upper decomposition of the DomainMatrix
Returns
=======
(L, U, exchange)
L, U are Lower and Upper decomposition of the DomainMatrix,
exchange is the list of indices of rows exchanged in the
decomposition.
Raises
======
ValueError
If the domain of DomainMatrix not a Field
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [QQ(1), QQ(-1)],
... [QQ(2), QQ(-2)]], (2, 2), QQ)
>>> L, U, exchange = A.lu()
>>> L
DomainMatrix([[1, 0], [2, 1]], (2, 2), QQ)
>>> U
DomainMatrix([[1, -1], [0, 0]], (2, 2), QQ)
>>> exchange
[]
See Also
========
lu_solve
"""
if not self.domain.is_Field:
raise DMNotAField('Not a field')
L, U, swaps = self.rep.lu()
return self.from_rep(L), self.from_rep(U), swaps
def qr(self):
r"""
QR decomposition of the DomainMatrix.
Explanation
===========
The QR decomposition expresses a matrix as the product of an orthogonal
matrix (Q) and an upper triangular matrix (R). In this implementation,
Q is not orthonormal: its columns are orthogonal but not normalized to
unit vectors. This avoids unnecessary divisions and is particularly
suited for exact arithmetic domains.
Note
====
This implementation is valid only for matrices over real domains. For
matrices over complex domains, a proper QR decomposition would require
handling conjugation to ensure orthogonality.
Returns
=======
(Q, R)
Q is the orthogonal matrix, and R is the upper triangular matrix
resulting from the QR decomposition of the DomainMatrix.
Raises
======
DMDomainError
If the domain of the DomainMatrix is not a field (e.g., QQ).
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([[1, 2], [3, 4], [5, 6]], (3, 2), QQ)
>>> Q, R = A.qr()
>>> Q
DomainMatrix([[1, 26/35], [3, 8/35], [5, -2/7]], (3, 2), QQ)
>>> R
DomainMatrix([[1, 44/35], [0, 1]], (2, 2), QQ)
>>> Q * R == A
True
>>> (Q.transpose() * Q).is_diagonal
True
>>> R.is_upper
True
See Also
========
lu
"""
ddm_q, ddm_r = self.rep.qr()
Q = self.from_rep(ddm_q)
R = self.from_rep(ddm_r)
return Q, R
def lu_solve(self, rhs):
r"""
Solver for DomainMatrix x in the A*x = B
Parameters
==========
rhs : DomainMatrix B
Returns
=======
DomainMatrix
x in A*x = B
Raises
======
DMShapeError
If the DomainMatrix A and rhs have different number of rows
ValueError
If the domain of DomainMatrix A not a Field
Examples
========
>>> from sympy import QQ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [QQ(1), QQ(2)],
... [QQ(3), QQ(4)]], (2, 2), QQ)
>>> B = DomainMatrix([
... [QQ(1), QQ(1)],
... [QQ(0), QQ(1)]], (2, 2), QQ)
>>> A.lu_solve(B)
DomainMatrix([[-2, -1], [3/2, 1]], (2, 2), QQ)
See Also
========
lu
"""
if self.shape[0] != rhs.shape[0]:
raise DMShapeError("Shape")
if not self.domain.is_Field:
raise DMNotAField('Not a field')
sol = self.rep.lu_solve(rhs.rep)
return self.from_rep(sol)
def fflu(self):
"""
Fraction-free LU decomposition of DomainMatrix.
Explanation
===========
This method computes the PLDU decomposition
using Gauss-Bareiss elimination in a fraction-free manner,
it ensures that all intermediate results remain in
the domain of the input matrix. Unlike standard
LU decomposition, which introduces division, this approach
avoids fractions, making it particularly suitable
for exact arithmetic over integers or polynomials.
This method satisfies the invariant:
P * A = L * inv(D) * U
Returns
=======
(P, L, D, U)
- P (Permutation matrix)
- L (Lower triangular matrix)
- D (Diagonal matrix)
- U (Upper triangular matrix)
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ)
>>> P, L, D, U = A.fflu()
>>> P
DomainMatrix([[1, 0], [0, 1]], (2, 2), ZZ)
>>> L
DomainMatrix([[1, 0], [3, -2]], (2, 2), ZZ)
>>> D
DomainMatrix([[1, 0], [0, -2]], (2, 2), ZZ)
>>> U
DomainMatrix([[1, 2], [0, -2]], (2, 2), ZZ)
>>> L.is_lower and U.is_upper and D.is_diagonal
True
>>> L * D.to_field().inv() * U == P * A.to_field()
True
>>> I, d = D.inv_den()
>>> L * I * U == d * P * A
True
See Also
========
sympy.polys.matrices.ddm.DDM.fflu
References
==========
.. [1] Nakos, G. C., Turner, P. R., & Williams, R. M. (1997). Fraction-free
algorithms for linear and polynomial equations. ACM SIGSAM Bulletin,
31(3), 11-19. https://doi.org/10.1145/271130.271133
.. [2] Middeke, J.; Jeffrey, D.J.; Koutschan, C. (2020), "Common Factors
in Fraction-Free Matrix Decompositions", Mathematics in Computer Science,
15 (4): 589–608, arXiv:2005.12380, doi:10.1007/s11786-020-00495-9
.. [3] https://en.wikipedia.org/wiki/Bareiss_algorithm
"""
from_rep = self.from_rep
P, L, D, U = self.rep.fflu()
return from_rep(P), from_rep(L), from_rep(D), from_rep(U)
def _solve(A, b):
# XXX: Not sure about this method or its signature. It is just created
# because it is needed by the holonomic module.
if A.shape[0] != b.shape[0]:
raise DMShapeError("Shape")
if A.domain != b.domain or not A.domain.is_Field:
raise DMNotAField('Not a field')
Aaug = A.hstack(b)
Arref, pivots = Aaug.rref()
particular = Arref.from_rep(Arref.rep.particular())
nullspace_rep, nonpivots = Arref[:,:-1].rep.nullspace()
nullspace = Arref.from_rep(nullspace_rep)
return particular, nullspace
def charpoly(self):
r"""
Characteristic polynomial of a square matrix.
Computes the characteristic polynomial in a fully expanded form using
division free arithmetic. If a factorization of the characteristic
polynomial is needed then it is more efficient to call
:meth:`charpoly_factor_list` than calling :meth:`charpoly` and then
factorizing the result.
Returns
=======
list: list of DomainElement
coefficients of the characteristic polynomial
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> A.charpoly()
[1, -5, -2]
See Also
========
charpoly_factor_list
Compute the factorisation of the characteristic polynomial.
charpoly_factor_blocks
A partial factorisation of the characteristic polynomial that can
be computed more efficiently than either the full factorisation or
the fully expanded polynomial.
"""
M = self
K = M.domain
factors = M.charpoly_factor_blocks()
cp = [K.one]
for f, mult in factors:
for _ in range(mult):
cp = dup_mul(cp, f, K)
return cp
def charpoly_factor_list(self):
"""
Full factorization of the characteristic polynomial.
Examples
========
>>> from sympy.polys.matrices import DM
>>> from sympy import ZZ
>>> M = DM([[6, -1, 0, 0],
... [9, 12, 0, 0],
... [0, 0, 1, 2],
... [0, 0, 5, 6]], ZZ)
Compute the factorization of the characteristic polynomial:
>>> M.charpoly_factor_list()
[([1, -9], 2), ([1, -7, -4], 1)]
Use :meth:`charpoly` to get the unfactorized characteristic polynomial:
>>> M.charpoly()
[1, -25, 203, -495, -324]
The same calculations with ``Matrix``:
>>> M.to_Matrix().charpoly().as_expr()
lambda**4 - 25*lambda**3 + 203*lambda**2 - 495*lambda - 324
>>> M.to_Matrix().charpoly().as_expr().factor()
(lambda - 9)**2*(lambda**2 - 7*lambda - 4)
Returns
=======
list: list of pairs (factor, multiplicity)
A full factorization of the characteristic polynomial.
See Also
========
charpoly
Expanded form of the characteristic polynomial.
charpoly_factor_blocks
A partial factorisation of the characteristic polynomial that can
be computed more efficiently.
"""
M = self
K = M.domain
# It is more efficient to start from the partial factorization provided
# for free by M.charpoly_factor_blocks than the expanded M.charpoly.
factors = M.charpoly_factor_blocks()
factors_irreducible = []
for factor_i, mult_i in factors:
_, factors_list = dup_factor_list(factor_i, K)
for factor_j, mult_j in factors_list:
factors_irreducible.append((factor_j, mult_i * mult_j))
return _collect_factors(factors_irreducible)
def charpoly_factor_blocks(self):
"""
Partial factorisation of the characteristic polynomial.
This factorisation arises from a block structure of the matrix (if any)
and so the factors are not guaranteed to be irreducible. The
:meth:`charpoly_factor_blocks` method is the most efficient way to get
a representation of the characteristic polynomial but the result is
neither fully expanded nor fully factored.
Examples
========
>>> from sympy.polys.matrices import DM
>>> from sympy import ZZ
>>> M = DM([[6, -1, 0, 0],
... [9, 12, 0, 0],
... [0, 0, 1, 2],
... [0, 0, 5, 6]], ZZ)
This computes a partial factorization using only the block structure of
the matrix to reveal factors:
>>> M.charpoly_factor_blocks()
[([1, -18, 81], 1), ([1, -7, -4], 1)]
These factors correspond to the two diagonal blocks in the matrix:
>>> DM([[6, -1], [9, 12]], ZZ).charpoly()
[1, -18, 81]
>>> DM([[1, 2], [5, 6]], ZZ).charpoly()
[1, -7, -4]
Use :meth:`charpoly_factor_list` to get a complete factorization into
irreducibles:
>>> M.charpoly_factor_list()
[([1, -9], 2), ([1, -7, -4], 1)]
Use :meth:`charpoly` to get the expanded characteristic polynomial:
>>> M.charpoly()
[1, -25, 203, -495, -324]
Returns
=======
list: list of pairs (factor, multiplicity)
A partial factorization of the characteristic polynomial.
See Also
========
charpoly
Compute the fully expanded characteristic polynomial.
charpoly_factor_list
Compute a full factorization of the characteristic polynomial.
"""
M = self
if not M.is_square:
raise DMNonSquareMatrixError("not square")
# scc returns indices that permute the matrix into block triangular
# form and can extract the diagonal blocks. M.charpoly() is equal to
# the product of the diagonal block charpolys.
components = M.scc()
block_factors = []
for indices in components:
block = M.extract(indices, indices)
block_factors.append((block.charpoly_base(), 1))
return _collect_factors(block_factors)
def charpoly_base(self):
"""
Base case for :meth:`charpoly_factor_blocks` after block decomposition.
This method is used internally by :meth:`charpoly_factor_blocks` as the
base case for computing the characteristic polynomial of a block. It is
more efficient to call :meth:`charpoly_factor_blocks`, :meth:`charpoly`
or :meth:`charpoly_factor_list` rather than call this method directly.
This will use either the dense or the sparse implementation depending
on the sparsity of the matrix and will clear denominators if possible
before calling :meth:`charpoly_berk` to compute the characteristic
polynomial using the Berkowitz algorithm.
See Also
========
charpoly
charpoly_factor_list
charpoly_factor_blocks
charpoly_berk
"""
M = self
K = M.domain
# It seems that the sparse implementation is always faster for random
# matrices with fewer than 50% non-zero entries. This does not seem to
# depend on domain, size, bit count etc.
density = self.nnz() / self.shape[0]**2
if density < 0.5:
M = M.to_sparse()
else:
M = M.to_dense()
# Clearing denominators is always more efficient if it can be done.
# Doing it here after block decomposition is good because each block
# might have a smaller denominator. However it might be better for
# charpoly and charpoly_factor_list to restore the denominators only at
# the very end so that they can call e.g. dup_factor_list before
# restoring the denominators. The methods would need to be changed to
# return (poly, denom) pairs to make that work though.
clear_denoms = K.is_Field and K.has_assoc_Ring
if clear_denoms:
clear_denoms = True
d, M = M.clear_denoms(convert=True)
d = d.element
K_f = K
K_r = M.domain
# Berkowitz algorithm over K_r.
cp = M.charpoly_berk()
if clear_denoms:
# Restore the denominator in the charpoly over K_f.
#
# If M = N/d then p_M(x) = p_N(x*d)/d^n.
cp = dup_convert(cp, K_r, K_f)
p = [K_f.one, K_f.zero]
q = [K_f.one/d]
cp = dup_transform(cp, p, q, K_f)
return cp
def charpoly_berk(self):
"""Compute the characteristic polynomial using the Berkowitz algorithm.
This method directly calls the underlying implementation of the
Berkowitz algorithm (:meth:`sympy.polys.matrices.dense.ddm_berk` or
:meth:`sympy.polys.matrices.sdm.sdm_berk`).
This is used by :meth:`charpoly` and other methods as the base case for
for computing the characteristic polynomial. However those methods will
apply other optimizations such as block decomposition, clearing
denominators and converting between dense and sparse representations
before calling this method. It is more efficient to call those methods
instead of this one but this method is provided for direct access to
the Berkowitz algorithm.
Examples
========
>>> from sympy.polys.matrices import DM
>>> from sympy import QQ
>>> M = DM([[6, -1, 0, 0],
... [9, 12, 0, 0],
... [0, 0, 1, 2],
... [0, 0, 5, 6]], QQ)
>>> M.charpoly_berk()
[1, -25, 203, -495, -324]
See Also
========
charpoly
charpoly_base
charpoly_factor_list
charpoly_factor_blocks
sympy.polys.matrices.dense.ddm_berk
sympy.polys.matrices.sdm.sdm_berk
"""
return self.rep.charpoly()
@classmethod
def eye(cls, shape, domain):
r"""
Return identity matrix of size n or shape (m, n).
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy import QQ
>>> DomainMatrix.eye(3, QQ)
DomainMatrix({0: {0: 1}, 1: {1: 1}, 2: {2: 1}}, (3, 3), QQ)
"""
if isinstance(shape, int):
shape = (shape, shape)
return cls.from_rep(SDM.eye(shape, domain))
@classmethod
def diag(cls, diagonal, domain, shape=None):
r"""
Return diagonal matrix with entries from ``diagonal``.
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy import ZZ
>>> DomainMatrix.diag([ZZ(5), ZZ(6)], ZZ)
DomainMatrix({0: {0: 5}, 1: {1: 6}}, (2, 2), ZZ)
"""
if shape is None:
N = len(diagonal)
shape = (N, N)
return cls.from_rep(SDM.diag(diagonal, domain, shape))
@classmethod
def zeros(cls, shape, domain, *, fmt='sparse'):
"""Returns a zero DomainMatrix of size shape, belonging to the specified domain
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy import QQ
>>> DomainMatrix.zeros((2, 3), QQ)
DomainMatrix({}, (2, 3), QQ)
"""
return cls.from_rep(SDM.zeros(shape, domain))
@classmethod
def ones(cls, shape, domain):
"""Returns a DomainMatrix of 1s, of size shape, belonging to the specified domain
Examples
========
>>> from sympy.polys.matrices import DomainMatrix
>>> from sympy import QQ
>>> DomainMatrix.ones((2,3), QQ)
DomainMatrix([[1, 1, 1], [1, 1, 1]], (2, 3), QQ)
"""
return cls.from_rep(DDM.ones(shape, domain).to_dfm_or_ddm())
def __eq__(A, B):
r"""
Checks for two DomainMatrix matrices to be equal or not
Parameters
==========
A, B: DomainMatrix
to check equality
Returns
=======
Boolean
True for equal, else False
Raises
======
NotImplementedError
If B is not a DomainMatrix
Examples
========
>>> from sympy import ZZ
>>> from sympy.polys.matrices import DomainMatrix
>>> A = DomainMatrix([
... [ZZ(1), ZZ(2)],
... [ZZ(3), ZZ(4)]], (2, 2), ZZ)
>>> B = DomainMatrix([
... [ZZ(1), ZZ(1)],
... [ZZ(0), ZZ(1)]], (2, 2), ZZ)
>>> A.__eq__(A)
True
>>> A.__eq__(B)
False
"""
if not isinstance(A, type(B)):
return NotImplemented
return A.domain == B.domain and A.rep == B.rep
def unify_eq(A, B):
if A.shape != B.shape:
return False
if A.domain != B.domain:
A, B = A.unify(B)
return A == B
def lll(A, delta=QQ(3, 4)):
"""
Performs the Lenstra–Lenstra–Lovász (LLL) basis reduction algorithm.
See [1]_ and [2]_.
Parameters
==========
delta : QQ, optional
The Lovász parameter. Must be in the interval (0.25, 1), with larger
values producing a more reduced basis. The default is 0.75 for
historical reasons.
Returns
=======
The reduced basis as a DomainMatrix over ZZ.
Throws
======
DMValueError: if delta is not in the range (0.25, 1)
DMShapeError: if the matrix is not of shape (m, n) with m <= n
DMDomainError: if the matrix domain is not ZZ
DMRankError: if the matrix contains linearly dependent rows
Examples
========
>>> from sympy.polys.domains import ZZ, QQ
>>> from sympy.polys.matrices import DM
>>> x = DM([[1, 0, 0, 0, -20160],
... [0, 1, 0, 0, 33768],
... [0, 0, 1, 0, 39578],
... [0, 0, 0, 1, 47757]], ZZ)
>>> y = DM([[10, -3, -2, 8, -4],
... [3, -9, 8, 1, -11],
... [-3, 13, -9, -3, -9],
... [-12, -7, -11, 9, -1]], ZZ)
>>> assert x.lll(delta=QQ(5, 6)) == y
Notes
=====
The implementation is derived from the Maple code given in Figures 4.3
and 4.4 of [3]_ (pp.68-69). It uses the efficient method of only calculating
state updates as they are required.
See also
========
lll_transform
References
==========
.. [1] https://en.wikipedia.org/wiki/Lenstra%E2%80%93Lenstra%E2%80%93Lov%C3%A1sz_lattice_basis_reduction_algorithm
.. [2] https://web.archive.org/web/20221029115428/https://web.cs.elte.hu/~lovasz/scans/lll.pdf
.. [3] Murray R. Bremner, "Lattice Basis Reduction: An Introduction to the LLL Algorithm and Its Applications"
"""
return DomainMatrix.from_rep(A.rep.lll(delta=delta))
def lll_transform(A, delta=QQ(3, 4)):
"""
Performs the Lenstra–Lenstra–Lovász (LLL) basis reduction algorithm
and returns the reduced basis and transformation matrix.
Explanation
===========
Parameters, algorithm and basis are the same as for :meth:`lll` except that
the return value is a tuple `(B, T)` with `B` the reduced basis and
`T` a transformation matrix. The original basis `A` is transformed to
`B` with `T*A == B`. If only `B` is needed then :meth:`lll` should be
used as it is a little faster.
Examples
========
>>> from sympy.polys.domains import ZZ, QQ
>>> from sympy.polys.matrices import DM
>>> X = DM([[1, 0, 0, 0, -20160],
... [0, 1, 0, 0, 33768],
... [0, 0, 1, 0, 39578],
... [0, 0, 0, 1, 47757]], ZZ)
>>> B, T = X.lll_transform(delta=QQ(5, 6))
>>> T * X == B
True
See also
========
lll
"""
reduced, transform = A.rep.lll_transform(delta=delta)
return DomainMatrix.from_rep(reduced), DomainMatrix.from_rep(transform)
def _collect_factors(factors_list):
"""
Collect repeating factors and sort.
>>> from sympy.polys.matrices.domainmatrix import _collect_factors
>>> _collect_factors([([1, 2], 2), ([1, 4], 3), ([1, 2], 5)])
[([1, 4], 3), ([1, 2], 7)]
"""
factors = Counter()
for factor, exponent in factors_list:
factors[tuple(factor)] += exponent
factors_list = [(list(f), e) for f, e in factors.items()]
return _sort_factors(factors_list)
Reference in a new issue Copy permalink