From fa986676deaf95c6daf7cd2ab5e9c9db37cda5fc Mon Sep 17 00:00:00 2001 From: Oscar Benjamin Date: Thu, 7 Jan 2021 19:26:04 +0000 Subject: [PATCH 1/6] maint(polys): add sympy.polys.matrices package Refactor domainmatrix.py into a package sympy.polys.matrices --- sympy/matrices/eigen.py | 4 +- sympy/polys/domainmatrix.py | 791 +----------- sympy/polys/matrices/__init__.py | 5 + sympy/polys/matrices/ddm.py | 204 ++++ sympy/polys/matrices/dense.py | 301 +++++ sympy/polys/matrices/domainmatrix.py | 202 +++ sympy/polys/matrices/eigen.py | 85 ++ sympy/polys/matrices/exceptions.py | 29 + sympy/polys/matrices/tests/__init__.py | 0 sympy/polys/matrices/tests/test_ddm.py | 341 ++++++ sympy/polys/matrices/tests/test_dense.py | 341 ++++++ .../polys/matrices/tests/test_domainmatrix.py | 422 +++++++ sympy/polys/tests/test_domainmatrix.py | 1085 ----------------- 13 files changed, 1937 insertions(+), 1873 deletions(-) create mode 100644 sympy/polys/matrices/__init__.py create mode 100644 sympy/polys/matrices/ddm.py create mode 100644 sympy/polys/matrices/dense.py create mode 100644 sympy/polys/matrices/domainmatrix.py create mode 100644 sympy/polys/matrices/eigen.py create mode 100644 sympy/polys/matrices/exceptions.py create mode 100644 sympy/polys/matrices/tests/__init__.py create mode 100644 sympy/polys/matrices/tests/test_ddm.py create mode 100644 sympy/polys/matrices/tests/test_dense.py create mode 100644 sympy/polys/matrices/tests/test_domainmatrix.py delete mode 100644 sympy/polys/tests/test_domainmatrix.py diff --git a/sympy/matrices/eigen.py b/sympy/matrices/eigen.py index 04efdda461a4..829d1f7542bc 100644 --- a/sympy/matrices/eigen.py +++ b/sympy/matrices/eigen.py @@ -11,8 +11,8 @@ from sympy.core.sympify import _sympify from sympy.functions.elementary.miscellaneous import sqrt from sympy.polys import roots, CRootOf, EX -from sympy.polys.domainmatrix import ( - DomainMatrix, dom_eigenvects, dom_eigenvects_to_sympy) +from sympy.polys.matrices import DomainMatrix +from sympy.polys.matrices.eigen import dom_eigenvects, dom_eigenvects_to_sympy from sympy.simplify import nsimplify, simplify as _simplify from sympy.utilities.exceptions import SymPyDeprecationWarning diff --git a/sympy/polys/domainmatrix.py b/sympy/polys/domainmatrix.py index dd97d0ba8b3b..2414264bdcf6 100644 --- a/sympy/polys/domainmatrix.py +++ b/sympy/polys/domainmatrix.py @@ -1,787 +1,6 @@ -from operator import mul +# +# Stub to expose DomainMatrix which has now moved to +# sympy.polys.matrices package. It should really be imported from there. +# -from sympy.core.symbol import Dummy -from sympy.core.sympify import _sympify - -from sympy.matrices.common import (NonInvertibleMatrixError, - NonSquareMatrixError, ShapeError) -from sympy.polys import Poly -from sympy.polys.agca.extensions import FiniteExtension -from sympy.polys.constructor import construct_domain -from sympy.polys.factortools import dup_factor_list -from sympy.polys.polyroots import roots -from sympy.polys.rootoftools import CRootOf - - -class DDMError(Exception): - """Base class for errors raised by DDM""" - pass - - -class DDMBadInputError(DDMError): - """list of lists is inconsistent with shape""" - pass - - -class DDMDomainError(DDMError): - """domains do not match""" - pass - - -class DDMShapeError(DDMError): - """shapes are inconsistent""" - pass - - -class DDM(list): - """Dense matrix based on polys domain elements - - This is a list subclass and is a wrapper for a list of lists that supports - basic matrix arithmetic +, -, *, **. - """ - def __init__(self, rowslist, shape, domain): - super().__init__(rowslist) - self.shape = self.rows, self.cols = m, n = shape - self.domain = domain - - if not (len(self) == m and all(len(row) == n for row in self)): - raise DDMBadInputError("Inconsistent row-list/shape") - - def __str__(self): - cls = type(self).__name__ - rows = list.__str__(self) - return '%s(%s, %s, %s)' % (cls, rows, self.shape, self.domain) - - def __eq__(self, other): - if not isinstance(other, DDM): - return False - return (super().__eq__(other) and self.domain == other.domain) - - def __ne__(self, other): - return not self.__eq__(other) - - @classmethod - def zeros(cls, shape, domain): - z = domain.zero - m, n = shape - rowslist = ([z] * n for _ in range(m)) - return DDM(rowslist, shape, domain) - - @classmethod - def eye(cls, size, domain): - one = domain.one - ddm = cls.zeros((size, size), domain) - for i in range(size): - ddm[i][i] = one - return ddm - - def copy(self): - copyrows = (row[:] for row in self) - return DDM(copyrows, self.shape, self.domain) - - def __add__(a, b): - if not isinstance(b, DDM): - return NotImplemented - return a.add(b) - - def __sub__(a, b): - if not isinstance(b, DDM): - return NotImplemented - return a.sub(b) - - def __neg__(a): - return a.neg() - - def __mul__(a, b): - if b in a.domain: - return a.mul(b) - else: - return NotImplemented - - def __matmul__(a, b): - if isinstance(b, DDM): - return a.matmul(b) - else: - return NotImplemented - - @classmethod - def _check(cls, a, op, b, ashape, bshape): - if a.domain != b.domain: - msg = "Domain mismatch: %s %s %s" % (a.domain, op, b.domain) - raise DDMDomainError(msg) - if ashape != bshape: - msg = "Shape mismatch: %s %s %s" % (a.shape, op, b.shape) - raise DDMShapeError(msg) - - def add(a, b): - """a + b""" - a._check(a, '+', b, a.shape, b.shape) - c = a.copy() - ddm_iadd(c, b) - return c - - def sub(a, b): - """a - b""" - a._check(a, '-', b, a.shape, b.shape) - c = a.copy() - ddm_isub(c, b) - return c - - def neg(a): - """-a""" - b = a.copy() - ddm_ineg(b) - return b - - def mul(a, b): - c = a.copy() - ddm_imul(c, b) - return c - - def matmul(a, b): - """a @ b (matrix product)""" - m, o = a.shape - o2, n = b.shape - a._check(a, '*', b, o, o2) - c = a.zeros((m, n), a.domain) - ddm_imatmul(c, a, b) - return c - - def rref(a): - """Reduced-row echelon form of a and list of pivots""" - b = a.copy() - pivots = ddm_irref(b) - return b, pivots - - def nullspace(a): - rref, pivots = a.rref() - rows, cols = a.shape - domain = a.domain - - basis = [] - for i in range(cols): - if i in pivots: - continue - vec = [domain.one if i == j else domain.zero for j in range(cols)] - for ii, jj in enumerate(pivots): - vec[jj] -= rref[ii][i] - basis.append(vec) - - return DDM(basis, (len(basis), cols), domain) - - def det(a): - """Determinant of a""" - m, n = a.shape - if m != n: - raise DDMShapeError("Determinant of non-square matrix") - b = a.copy() - K = b.domain - deta = ddm_idet(b, K) - return deta - - def inv(a): - """Inverse of a""" - m, n = a.shape - if m != n: - raise DDMShapeError("Determinant of non-square matrix") - ainv = a.copy() - K = a.domain - ddm_iinv(ainv, a, K) - return ainv - - def lu(a): - """L, U decomposition of a""" - m, n = a.shape - K = a.domain - - U = a.copy() - L = a.eye(m, K) - swaps = ddm_ilu_split(L, U, K) - - return L, U, swaps - - def lu_solve(a, b): - """x where a*x = b""" - m, n = a.shape - m2, o = b.shape - a._check(a, 'lu_solve', b, m, m2) - - L, U, swaps = a.lu() - x = a.zeros((n, o), a.domain) - ddm_ilu_solve(x, L, U, swaps, b) - return x - - def charpoly(a): - """Coefficients of characteristic polynomial of a""" - K = a.domain - m, n = a.shape - if m != n: - raise DDMShapeError("Charpoly of non-square matrix") - vec = ddm_berk(a, K) - coeffs = [vec[i][0] for i in range(n+1)] - return coeffs - - -def ddm_iadd(a, b): - """a += b""" - for ai, bi in zip(a, b): - for j, bij in enumerate(bi): - ai[j] += bij - - -def ddm_isub(a, b): - """a -= b""" - for ai, bi in zip(a, b): - for j, bij in enumerate(bi): - ai[j] -= bij - - -def ddm_ineg(a): - """a <-- -a""" - for ai in a: - for j, aij in enumerate(ai): - ai[j] = -aij - - -def ddm_imul(a, b): - for ai in a: - for j, aij in enumerate(ai): - ai[j] = b * aij - - -def ddm_imatmul(a, b, c): - """a += b @ c""" - cT = list(zip(*c)) - - for bi, ai in zip(b, a): - for j, cTj in enumerate(cT): - ai[j] = sum(map(mul, bi, cTj), ai[j]) - - -def ddm_irref(a): - """a <-- rref(a)""" - # a is (m x n) - m = len(a) - if not m: - return [] - n = len(a[0]) - - i = 0 - pivots = [] - - for j in range(n): - # pivot - aij = a[i][j] - - # zero-pivot - if not aij: - for ip in range(i+1, m): - aij = a[ip][j] - # row-swap - if aij: - a[i], a[ip] = a[ip], a[i] - break - else: - # next column - continue - - # normalise row - ai = a[i] - aijinv = aij**-1 - for l in range(j, n): - ai[l] *= aijinv # ai[j] = one - - # eliminate above and below to the right - for k, ak in enumerate(a): - if k == i or not ak[j]: - continue - akj = ak[j] - ak[j] -= akj # ak[j] = zero - for l in range(j+1, n): - ak[l] -= akj * ai[l] - - # next row - pivots.append(j) - i += 1 - - # no more rows? - if i >= m: - break - - return pivots - - -def ddm_idet(a, K): - """a <-- echelon(a); return det""" - # Fraction-free Gaussian elimination - # https://www.math.usm.edu/perry/Research/Thesis_DRL.pdf - - # a is (m x n) - m = len(a) - if not m: - return K.one - n = len(a[0]) - - is_field = K.is_Field - # uf keeps track of the effect of row swaps and multiplies - uf = K.one - for j in range(n-1): - # if zero on the diagonal need to swap - if not a[j][j]: - for l in range(j+1, n): - if a[l][j]: - a[j], a[l] = a[l], a[j] - uf = -uf - break - else: - # unable to swap: det = 0 - return K.zero - for i in range(j+1, n): - if a[i][j]: - if not is_field: - d = K.gcd(a[j][j], a[i][j]) - b = a[j][j] // d - c = a[i][j] // d - else: - b = a[j][j] - c = a[i][j] - # account for multiplying row i by b - uf = b * uf - for k in range(j+1, n): - a[i][k] = b*a[i][k] - c*a[j][k] - - # triangular det is product of diagonal - prod = K.one - for i in range(n): - prod = prod * a[i][i] - # incorporate swaps and multiplies - if not is_field: - D = prod // uf - else: - D = prod / uf - return D - - -def ddm_iinv(ainv, a, K): - if not K.is_Field: - raise ValueError('Not a field') - - # a is (m x n) - m = len(a) - if not m: - return - n = len(a[0]) - if m != n: - raise NonSquareMatrixError - - eye = [[K.one if i==j else K.zero for j in range(n)] for i in range(n)] - Aaug = [row + eyerow for row, eyerow in zip(a, eye)] - pivots = ddm_irref(Aaug) - if pivots != list(range(n)): - raise NonInvertibleMatrixError('Matrix det == 0; not invertible.') - ainv[:] = [row[n:] for row in Aaug] - - -def ddm_ilu_split(L, U, K): - """L, U <-- LU(U)""" - m = len(U) - if not m: - return [] - n = len(U[0]) - - swaps = ddm_ilu(U) - - zeros = [K.zero] * min(m, n) - for i in range(1, m): - j = min(i, n) - L[i][:j] = U[i][:j] - U[i][:j] = zeros[:j] - - return swaps - - -def ddm_ilu(a): - """a <-- LU(a)""" - m = len(a) - if not m: - return [] - n = len(a[0]) - - swaps = [] - - for i in range(min(m, n)): - if not a[i][i]: - for ip in range(i+1, m): - if a[ip][i]: - swaps.append((i, ip)) - a[i], a[ip] = a[ip], a[i] - break - else: - # M = Matrix([[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 1, 2]]) - continue - for j in range(i+1, m): - l_ji = a[j][i] / a[i][i] - a[j][i] = l_ji - for k in range(i+1, n): - a[j][k] -= l_ji * a[i][k] - - return swaps - - -def ddm_ilu_solve(x, L, U, swaps, b): - """x <-- solve(L*U*x = swaps(b))""" - m = len(U) - if not m: - return - n = len(U[0]) - - m2 = len(b) - if not m2: - raise DDMShapeError("Shape mismtch") - o = len(b[0]) - - if m != m2: - raise DDMShapeError("Shape mismtch") - if m < n: - raise NotImplementedError("Underdetermined") - - if swaps: - b = [row[:] for row in b] - for i1, i2 in swaps: - b[i1], b[i2] = b[i2], b[i1] - - # solve Ly = b - y = [[None] * o for _ in range(m)] - for k in range(o): - for i in range(m): - rhs = b[i][k] - for j in range(i): - rhs -= L[i][j] * y[j][k] - y[i][k] = rhs - - if m > n: - for i in range(n, m): - for j in range(o): - if y[i][j]: - raise NonInvertibleMatrixError - - # Solve Ux = y - for k in range(o): - for i in reversed(range(n)): - if not U[i][i]: - raise NonInvertibleMatrixError - rhs = y[i][k] - for j in range(i+1, n): - rhs -= U[i][j] * x[j][k] - x[i][k] = rhs / U[i][i] - - -def ddm_berk(M, K): - m = len(M) - if not m: - return [[K.one]] - n = len(M[0]) - - if m != n: - raise DDMShapeError("Not square") - - if n == 1: - return [[K.one], [-M[0][0]]] - - a = M[0][0] - R = [M[0][1:]] - C = [[row[0]] for row in M[1:]] - A = [row[1:] for row in M[1:]] - - q = ddm_berk(A, K) - - T = [[K.zero] * n for _ in range(n+1)] - for i in range(n): - T[i][i] = K.one - T[i+1][i] = -a - for i in range(2, n+1): - if i == 2: - AnC = C - else: - C = AnC - AnC = [[K.zero] for row in C] - ddm_imatmul(AnC, A, C) - RAnC = [[K.zero]] - ddm_imatmul(RAnC, R, AnC) - for j in range(0, n+1-i): - T[i+j][j] = -RAnC[0][0] - - qout = [[K.zero] for _ in range(n+1)] - ddm_imatmul(qout, T, q) - return qout - - -class DomainMatrix: - - def __init__(self, rows, shape, domain): - self.rep = DDM(rows, shape, domain) - self.shape = shape - self.domain = domain - - @classmethod - def from_ddm(cls, ddm): - return cls(ddm, ddm.shape, ddm.domain) - - @classmethod - def from_list_sympy(cls, nrows, ncols, rows, **kwargs): - 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_Matrix(cls, M, **kwargs): - return cls.from_list_sympy(*M.shape, M.tolist(), **kwargs) - - @classmethod - def get_domain(cls, items_sympy, **kwargs): - K, items_K = construct_domain(items_sympy, **kwargs) - return K, items_K - - def convert_to(self, K): - Kold = self.domain - if K == Kold: - return self.from_ddm(self.rep.copy()) - new_rows = [[K.convert_from(e, Kold) for e in row] for row in self.rep] - return DomainMatrix(new_rows, self.shape, K) - - def to_field(self): - K = self.domain.get_field() - return self.convert_to(K) - - def unify(self, other): - K1 = self.domain - K2 = other.domain - if K1 == K2: - return self, other - K = K1.unify(K2) - if K1 != K: - self = self.convert_to(K) - if K2 != K: - other = other.convert_to(K) - return self, other - - def to_Matrix(self): - from sympy.matrices.dense import MutableDenseMatrix - rows_sympy = [[self.domain.to_sympy(e) for e in row] for row in self.rep] - return MutableDenseMatrix(rows_sympy) - - def __repr__(self): - rows_str = ['[%s]' % (', '.join(map(str, row))) for row in self.rep] - rowstr = '[%s]' % ', '.join(rows_str) - return 'DomainMatrix(%s, %r, %r)' % (rowstr, self.shape, self.domain) - - def __add__(A, B): - if not isinstance(B, DomainMatrix): - return NotImplemented - return A.add(B) - - def __sub__(A, B): - if not isinstance(B, DomainMatrix): - return NotImplemented - return A.sub(B) - - def __neg__(A): - return A.neg() - - def __mul__(A, B): - """A * B""" - if isinstance(B, DomainMatrix): - return A.matmul(B) - elif B in A.domain: - return A.from_ddm(A.rep * B) - else: - return NotImplemented - - def __rmul__(A, B): - if B in A.domain: - return A.from_ddm(A.rep * B) - else: - return NotImplemented - - def __pow__(A, n): - """A ** n""" - if not isinstance(n, int): - return NotImplemented - return A.pow(n) - - def add(A, B): - if A.shape != B.shape: - raise ShapeError("shape") - if A.domain != B.domain: - raise ValueError("domain") - return A.from_ddm(A.rep.add(B.rep)) - - def sub(A, B): - if A.shape != B.shape: - raise ShapeError("shape") - if A.domain != B.domain: - raise ValueError("domain") - return A.from_ddm(A.rep.sub(B.rep)) - - def neg(A): - return A.from_ddm(A.rep.neg()) - - def mul(A, b): - return A.from_ddm(A.rep.mul(b)) - - def matmul(A, B): - return A.from_ddm(A.rep.matmul(B.rep)) - - def pow(A, n): - if n < 0: - raise NotImplementedError('Negative powers') - elif n == 0: - m, n = A.shape - rows = [[A.domain.zero] * m for _ in range(m)] - for i in range(m): - rows[i][i] = A.domain.one - return type(A)(rows, A.shape, 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 rref(self): - if not self.domain.is_Field: - raise ValueError('Not a field') - rref_ddm, pivots = self.rep.rref() - return self.from_ddm(rref_ddm), tuple(pivots) - - def nullspace(self): - return self.from_ddm(self.rep.nullspace()) - - def inv(self): - if not self.domain.is_Field: - raise ValueError('Not a field') - m, n = self.shape - if m != n: - raise NonSquareMatrixError - inv = self.rep.inv() - return self.from_ddm(inv) - - def det(self): - m, n = self.shape - if m != n: - raise NonSquareMatrixError - return self.rep.det() - - def lu(self): - if not self.domain.is_Field: - raise ValueError('Not a field') - L, U, swaps = self.rep.lu() - return self.from_ddm(L), self.from_ddm(U), swaps - - def lu_solve(self, rhs): - if self.shape[0] != rhs.shape[0]: - raise ShapeError("Shape") - if not self.domain.is_Field: - raise ValueError('Not a field') - sol = self.rep.lu_solve(rhs.rep) - return self.from_ddm(sol) - - def charpoly(self): - m, n = self.shape - if m != n: - raise NonSquareMatrixError("not square") - return self.rep.charpoly() - - @classmethod - def eye(cls, n, domain): - return cls.from_ddm(DDM.eye(n, domain)) - - def __eq__(A, B): - """A == B""" - if not isinstance(B, DomainMatrix): - return NotImplemented - return A.rep == B.rep - - -def dom_eigenvects(A, l=Dummy('lambda')): - charpoly = A.charpoly() - rows, cols = A.shape - domain = A.domain - _, factors = dup_factor_list(charpoly, domain) - - rational_eigenvects = [] - algebraic_eigenvects = [] - for base, exp in factors: - if len(base) == 2: - field = domain - eigenval = -base[1] / base[0] - - EE_items = [ - [eigenval if i == j else field.zero for j in range(cols)] - for i in range(rows)] - EE = DomainMatrix(EE_items, (rows, cols), field) - - basis = (A - EE).nullspace() - rational_eigenvects.append((field, eigenval, exp, basis)) - else: - minpoly = Poly.from_list(base, l, domain=domain) - field = FiniteExtension(minpoly) - eigenval = field(l) - - AA_items = [ - [Poly.from_list([item], l, domain=domain).rep for item in row] - for row in A.rep] - AA_items = [[field(item) for item in row] for row in AA_items] - AA = DomainMatrix(AA_items, (rows, cols), field) - EE_items = [ - [eigenval if i == j else field.zero for j in range(cols)] - for i in range(rows)] - EE = DomainMatrix(EE_items, (rows, cols), field) - - basis = (AA - EE).nullspace() - algebraic_eigenvects.append((field, minpoly, exp, basis)) - - return rational_eigenvects, algebraic_eigenvects - - -def dom_eigenvects_to_sympy( - rational_eigenvects, algebraic_eigenvects, - Matrix, **kwargs -): - result = [] - - for field, eigenvalue, multiplicity, eigenvects in rational_eigenvects: - eigenvects = eigenvects.rep - eigenvalue = field.to_sympy(eigenvalue) - new_eigenvects = [ - Matrix([field.to_sympy(x) for x in vect]) - for vect in eigenvects] - result.append((eigenvalue, multiplicity, new_eigenvects)) - - for field, minpoly, multiplicity, eigenvects in algebraic_eigenvects: - eigenvects = eigenvects.rep - l = minpoly.gens[0] - - eigenvects = [[field.to_sympy(x) for x in vect] for vect in eigenvects] - - degree = minpoly.degree() - minpoly = minpoly.as_expr() - eigenvals = roots(minpoly, l, **kwargs) - if len(eigenvals) != degree: - eigenvals = [CRootOf(minpoly, l, idx) for idx in range(degree)] - - for eigenvalue in eigenvals: - new_eigenvects = [ - Matrix([x.subs(l, eigenvalue) for x in vect]) - for vect in eigenvects] - result.append((eigenvalue, multiplicity, new_eigenvects)) - - return result +from sympy.polys.matrices.domainmatrix import DomainMatrix diff --git a/sympy/polys/matrices/__init__.py b/sympy/polys/matrices/__init__.py new file mode 100644 index 000000000000..1d58d4174344 --- /dev/null +++ b/sympy/polys/matrices/__init__.py @@ -0,0 +1,5 @@ +from .domainmatrix import DomainMatrix + +__all__ = [ + 'DomainMatrix', +] diff --git a/sympy/polys/matrices/ddm.py b/sympy/polys/matrices/ddm.py new file mode 100644 index 000000000000..b0652fea4237 --- /dev/null +++ b/sympy/polys/matrices/ddm.py @@ -0,0 +1,204 @@ +from .exceptions import DDMBadInputError, DDMShapeError, DDMDomainError + +from .dense import ( + ddm_iadd, + ddm_isub, + ddm_ineg, + ddm_imul, + ddm_imatmul, + ddm_irref, + ddm_idet, + ddm_iinv, + ddm_ilu_split, + ddm_ilu_solve, + ddm_berk, + ) + + +class DDM(list): + """Dense matrix based on polys domain elements + + This is a list subclass and is a wrapper for a list of lists that supports + basic matrix arithmetic +, -, *, **. + """ + def __init__(self, rowslist, shape, domain): + super().__init__(rowslist) + self.shape = self.rows, self.cols = m, n = shape + self.domain = domain + + if not (len(self) == m and all(len(row) == n for row in self)): + raise DDMBadInputError("Inconsistent row-list/shape") + + def __str__(self): + cls = type(self).__name__ + rows = list.__str__(self) + return '%s(%s, %s, %s)' % (cls, rows, self.shape, self.domain) + + def __eq__(self, other): + if not isinstance(other, DDM): + return False + return (super().__eq__(other) and self.domain == other.domain) + + def __ne__(self, other): + return not self.__eq__(other) + + @classmethod + def zeros(cls, shape, domain): + z = domain.zero + m, n = shape + rowslist = ([z] * n for _ in range(m)) + return DDM(rowslist, shape, domain) + + @classmethod + def eye(cls, size, domain): + one = domain.one + ddm = cls.zeros((size, size), domain) + for i in range(size): + ddm[i][i] = one + return ddm + + def copy(self): + copyrows = (row[:] for row in self) + return DDM(copyrows, self.shape, self.domain) + + def __add__(a, b): + if not isinstance(b, DDM): + return NotImplemented + return a.add(b) + + def __sub__(a, b): + if not isinstance(b, DDM): + return NotImplemented + return a.sub(b) + + def __neg__(a): + return a.neg() + + def __mul__(a, b): + if b in a.domain: + return a.mul(b) + else: + return NotImplemented + + def __matmul__(a, b): + if isinstance(b, DDM): + return a.matmul(b) + else: + return NotImplemented + + @classmethod + def _check(cls, a, op, b, ashape, bshape): + if a.domain != b.domain: + msg = "Domain mismatch: %s %s %s" % (a.domain, op, b.domain) + raise DDMDomainError(msg) + if ashape != bshape: + msg = "Shape mismatch: %s %s %s" % (a.shape, op, b.shape) + raise DDMShapeError(msg) + + def add(a, b): + """a + b""" + a._check(a, '+', b, a.shape, b.shape) + c = a.copy() + ddm_iadd(c, b) + return c + + def sub(a, b): + """a - b""" + a._check(a, '-', b, a.shape, b.shape) + c = a.copy() + ddm_isub(c, b) + return c + + def neg(a): + """-a""" + b = a.copy() + ddm_ineg(b) + return b + + def mul(a, b): + c = a.copy() + ddm_imul(c, b) + return c + + def matmul(a, b): + """a @ b (matrix product)""" + m, o = a.shape + o2, n = b.shape + a._check(a, '*', b, o, o2) + c = a.zeros((m, n), a.domain) + ddm_imatmul(c, a, b) + return c + + def rref(a): + """Reduced-row echelon form of a and list of pivots""" + b = a.copy() + pivots = ddm_irref(b) + return b, pivots + + def nullspace(a): + rref, pivots = a.rref() + rows, cols = a.shape + domain = a.domain + + basis = [] + for i in range(cols): + if i in pivots: + continue + vec = [domain.one if i == j else domain.zero for j in range(cols)] + for ii, jj in enumerate(pivots): + vec[jj] -= rref[ii][i] + basis.append(vec) + + return DDM(basis, (len(basis), cols), domain) + + def det(a): + """Determinant of a""" + m, n = a.shape + if m != n: + raise DDMShapeError("Determinant of non-square matrix") + b = a.copy() + K = b.domain + deta = ddm_idet(b, K) + return deta + + def inv(a): + """Inverse of a""" + m, n = a.shape + if m != n: + raise DDMShapeError("Determinant of non-square matrix") + ainv = a.copy() + K = a.domain + ddm_iinv(ainv, a, K) + return ainv + + def lu(a): + """L, U decomposition of a""" + m, n = a.shape + K = a.domain + + U = a.copy() + L = a.eye(m, K) + swaps = ddm_ilu_split(L, U, K) + + return L, U, swaps + + def lu_solve(a, b): + """x where a*x = b""" + m, n = a.shape + m2, o = b.shape + a._check(a, 'lu_solve', b, m, m2) + + L, U, swaps = a.lu() + x = a.zeros((n, o), a.domain) + ddm_ilu_solve(x, L, U, swaps, b) + return x + + def charpoly(a): + """Coefficients of characteristic polynomial of a""" + K = a.domain + m, n = a.shape + if m != n: + raise DDMShapeError("Charpoly of non-square matrix") + vec = ddm_berk(a, K) + coeffs = [vec[i][0] for i in range(n+1)] + return coeffs diff --git a/sympy/polys/matrices/dense.py b/sympy/polys/matrices/dense.py new file mode 100644 index 000000000000..01d5dfc42fb5 --- /dev/null +++ b/sympy/polys/matrices/dense.py @@ -0,0 +1,301 @@ +from operator import mul + +from .exceptions import ( + DDMShapeError, + NonInvertibleMatrixError, + NonSquareMatrixError, + ) + + +def ddm_iadd(a, b): + """a += b""" + for ai, bi in zip(a, b): + for j, bij in enumerate(bi): + ai[j] += bij + + +def ddm_isub(a, b): + """a -= b""" + for ai, bi in zip(a, b): + for j, bij in enumerate(bi): + ai[j] -= bij + + +def ddm_ineg(a): + """a <-- -a""" + for ai in a: + for j, aij in enumerate(ai): + ai[j] = -aij + + +def ddm_imul(a, b): + for ai in a: + for j, aij in enumerate(ai): + ai[j] = b * aij + + +def ddm_imatmul(a, b, c): + """a += b @ c""" + cT = list(zip(*c)) + + for bi, ai in zip(b, a): + for j, cTj in enumerate(cT): + ai[j] = sum(map(mul, bi, cTj), ai[j]) + + +def ddm_irref(a): + """a <-- rref(a)""" + # a is (m x n) + m = len(a) + if not m: + return [] + n = len(a[0]) + + i = 0 + pivots = [] + + for j in range(n): + # pivot + aij = a[i][j] + + # zero-pivot + if not aij: + for ip in range(i+1, m): + aij = a[ip][j] + # row-swap + if aij: + a[i], a[ip] = a[ip], a[i] + break + else: + # next column + continue + + # normalise row + ai = a[i] + aijinv = aij**-1 + for l in range(j, n): + ai[l] *= aijinv # ai[j] = one + + # eliminate above and below to the right + for k, ak in enumerate(a): + if k == i or not ak[j]: + continue + akj = ak[j] + ak[j] -= akj # ak[j] = zero + for l in range(j+1, n): + ak[l] -= akj * ai[l] + + # next row + pivots.append(j) + i += 1 + + # no more rows? + if i >= m: + break + + return pivots + + +def ddm_idet(a, K): + """a <-- echelon(a); return det""" + # Fraction-free Gaussian elimination + # https://www.math.usm.edu/perry/Research/Thesis_DRL.pdf + + # a is (m x n) + m = len(a) + if not m: + return K.one + n = len(a[0]) + + is_field = K.is_Field + # uf keeps track of the effect of row swaps and multiplies + uf = K.one + for j in range(n-1): + # if zero on the diagonal need to swap + if not a[j][j]: + for l in range(j+1, n): + if a[l][j]: + a[j], a[l] = a[l], a[j] + uf = -uf + break + else: + # unable to swap: det = 0 + return K.zero + for i in range(j+1, n): + if a[i][j]: + if not is_field: + d = K.gcd(a[j][j], a[i][j]) + b = a[j][j] // d + c = a[i][j] // d + else: + b = a[j][j] + c = a[i][j] + # account for multiplying row i by b + uf = b * uf + for k in range(j+1, n): + a[i][k] = b*a[i][k] - c*a[j][k] + + # triangular det is product of diagonal + prod = K.one + for i in range(n): + prod = prod * a[i][i] + # incorporate swaps and multiplies + if not is_field: + D = prod // uf + else: + D = prod / uf + return D + + +def ddm_iinv(ainv, a, K): + if not K.is_Field: + raise ValueError('Not a field') + + # a is (m x n) + m = len(a) + if not m: + return + n = len(a[0]) + if m != n: + raise NonSquareMatrixError + + eye = [[K.one if i==j else K.zero for j in range(n)] for i in range(n)] + Aaug = [row + eyerow for row, eyerow in zip(a, eye)] + pivots = ddm_irref(Aaug) + if pivots != list(range(n)): + raise NonInvertibleMatrixError('Matrix det == 0; not invertible.') + ainv[:] = [row[n:] for row in Aaug] + + +def ddm_ilu_split(L, U, K): + """L, U <-- LU(U)""" + m = len(U) + if not m: + return [] + n = len(U[0]) + + swaps = ddm_ilu(U) + + zeros = [K.zero] * min(m, n) + for i in range(1, m): + j = min(i, n) + L[i][:j] = U[i][:j] + U[i][:j] = zeros[:j] + + return swaps + + +def ddm_ilu(a): + """a <-- LU(a)""" + m = len(a) + if not m: + return [] + n = len(a[0]) + + swaps = [] + + for i in range(min(m, n)): + if not a[i][i]: + for ip in range(i+1, m): + if a[ip][i]: + swaps.append((i, ip)) + a[i], a[ip] = a[ip], a[i] + break + else: + # M = Matrix([[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 1, 2]]) + continue + for j in range(i+1, m): + l_ji = a[j][i] / a[i][i] + a[j][i] = l_ji + for k in range(i+1, n): + a[j][k] -= l_ji * a[i][k] + + return swaps + + +def ddm_ilu_solve(x, L, U, swaps, b): + """x <-- solve(L*U*x = swaps(b))""" + m = len(U) + if not m: + return + n = len(U[0]) + + m2 = len(b) + if not m2: + raise DDMShapeError("Shape mismtch") + o = len(b[0]) + + if m != m2: + raise DDMShapeError("Shape mismtch") + if m < n: + raise NotImplementedError("Underdetermined") + + if swaps: + b = [row[:] for row in b] + for i1, i2 in swaps: + b[i1], b[i2] = b[i2], b[i1] + + # solve Ly = b + y = [[None] * o for _ in range(m)] + for k in range(o): + for i in range(m): + rhs = b[i][k] + for j in range(i): + rhs -= L[i][j] * y[j][k] + y[i][k] = rhs + + if m > n: + for i in range(n, m): + for j in range(o): + if y[i][j]: + raise NonInvertibleMatrixError + + # Solve Ux = y + for k in range(o): + for i in reversed(range(n)): + if not U[i][i]: + raise NonInvertibleMatrixError + rhs = y[i][k] + for j in range(i+1, n): + rhs -= U[i][j] * x[j][k] + x[i][k] = rhs / U[i][i] + + +def ddm_berk(M, K): + m = len(M) + if not m: + return [[K.one]] + n = len(M[0]) + + if m != n: + raise DDMShapeError("Not square") + + if n == 1: + return [[K.one], [-M[0][0]]] + + a = M[0][0] + R = [M[0][1:]] + C = [[row[0]] for row in M[1:]] + A = [row[1:] for row in M[1:]] + + q = ddm_berk(A, K) + + T = [[K.zero] * n for _ in range(n+1)] + for i in range(n): + T[i][i] = K.one + T[i+1][i] = -a + for i in range(2, n+1): + if i == 2: + AnC = C + else: + C = AnC + AnC = [[K.zero] for row in C] + ddm_imatmul(AnC, A, C) + RAnC = [[K.zero]] + ddm_imatmul(RAnC, R, AnC) + for j in range(0, n+1-i): + T[i+j][j] = -RAnC[0][0] + + qout = [[K.zero] for _ in range(n+1)] + ddm_imatmul(qout, T, q) + return qout diff --git a/sympy/polys/matrices/domainmatrix.py b/sympy/polys/matrices/domainmatrix.py new file mode 100644 index 000000000000..808fdd861919 --- /dev/null +++ b/sympy/polys/matrices/domainmatrix.py @@ -0,0 +1,202 @@ +from sympy.core.sympify import _sympify + +from ..constructor import construct_domain + +from .exceptions import NonSquareMatrixError, ShapeError + +from .ddm import DDM + + +class DomainMatrix: + + def __init__(self, rows, shape, domain): + self.rep = DDM(rows, shape, domain) + self.shape = shape + self.domain = domain + + @classmethod + def from_ddm(cls, ddm): + return cls(ddm, ddm.shape, ddm.domain) + + @classmethod + def from_list_sympy(cls, nrows, ncols, rows, **kwargs): + 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_Matrix(cls, M, **kwargs): + return cls.from_list_sympy(*M.shape, M.tolist(), **kwargs) + + @classmethod + def get_domain(cls, items_sympy, **kwargs): + K, items_K = construct_domain(items_sympy, **kwargs) + return K, items_K + + def convert_to(self, K): + Kold = self.domain + if K == Kold: + return self.from_ddm(self.rep.copy()) + new_rows = [[K.convert_from(e, Kold) for e in row] for row in self.rep] + return DomainMatrix(new_rows, self.shape, K) + + def to_field(self): + K = self.domain.get_field() + return self.convert_to(K) + + def unify(self, other): + K1 = self.domain + K2 = other.domain + if K1 == K2: + return self, other + K = K1.unify(K2) + if K1 != K: + self = self.convert_to(K) + if K2 != K: + other = other.convert_to(K) + return self, other + + def to_Matrix(self): + from sympy.matrices.dense import MutableDenseMatrix + rows_sympy = [[self.domain.to_sympy(e) for e in row] for row in self.rep] + return MutableDenseMatrix(rows_sympy) + + def __repr__(self): + rows_str = ['[%s]' % (', '.join(map(str, row))) for row in self.rep] + rowstr = '[%s]' % ', '.join(rows_str) + return 'DomainMatrix(%s, %r, %r)' % (rowstr, self.shape, self.domain) + + def __add__(A, B): + if not isinstance(B, DomainMatrix): + return NotImplemented + return A.add(B) + + def __sub__(A, B): + if not isinstance(B, DomainMatrix): + return NotImplemented + return A.sub(B) + + def __neg__(A): + return A.neg() + + def __mul__(A, B): + """A * B""" + if isinstance(B, DomainMatrix): + return A.matmul(B) + elif B in A.domain: + return A.from_ddm(A.rep * B) + else: + return NotImplemented + + def __rmul__(A, B): + if B in A.domain: + return A.from_ddm(A.rep * B) + else: + return NotImplemented + + def __pow__(A, n): + """A ** n""" + if not isinstance(n, int): + return NotImplemented + return A.pow(n) + + def add(A, B): + if A.shape != B.shape: + raise ShapeError("shape") + if A.domain != B.domain: + raise ValueError("domain") + return A.from_ddm(A.rep.add(B.rep)) + + def sub(A, B): + if A.shape != B.shape: + raise ShapeError("shape") + if A.domain != B.domain: + raise ValueError("domain") + return A.from_ddm(A.rep.sub(B.rep)) + + def neg(A): + return A.from_ddm(A.rep.neg()) + + def mul(A, b): + return A.from_ddm(A.rep.mul(b)) + + def matmul(A, B): + return A.from_ddm(A.rep.matmul(B.rep)) + + def pow(A, n): + if n < 0: + raise NotImplementedError('Negative powers') + elif n == 0: + m, n = A.shape + rows = [[A.domain.zero] * m for _ in range(m)] + for i in range(m): + rows[i][i] = A.domain.one + return type(A)(rows, A.shape, 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 rref(self): + if not self.domain.is_Field: + raise ValueError('Not a field') + rref_ddm, pivots = self.rep.rref() + return self.from_ddm(rref_ddm), tuple(pivots) + + def nullspace(self): + return self.from_ddm(self.rep.nullspace()) + + def inv(self): + if not self.domain.is_Field: + raise ValueError('Not a field') + m, n = self.shape + if m != n: + raise NonSquareMatrixError + inv = self.rep.inv() + return self.from_ddm(inv) + + def det(self): + m, n = self.shape + if m != n: + raise NonSquareMatrixError + return self.rep.det() + + def lu(self): + if not self.domain.is_Field: + raise ValueError('Not a field') + L, U, swaps = self.rep.lu() + return self.from_ddm(L), self.from_ddm(U), swaps + + def lu_solve(self, rhs): + if self.shape[0] != rhs.shape[0]: + raise ShapeError("Shape") + if not self.domain.is_Field: + raise ValueError('Not a field') + sol = self.rep.lu_solve(rhs.rep) + return self.from_ddm(sol) + + def charpoly(self): + m, n = self.shape + if m != n: + raise NonSquareMatrixError("not square") + return self.rep.charpoly() + + @classmethod + def eye(cls, n, domain): + return cls.from_ddm(DDM.eye(n, domain)) + + def __eq__(A, B): + """A == B""" + if not isinstance(B, DomainMatrix): + return NotImplemented + return A.rep == B.rep diff --git a/sympy/polys/matrices/eigen.py b/sympy/polys/matrices/eigen.py new file mode 100644 index 000000000000..6cbbd4fbf26e --- /dev/null +++ b/sympy/polys/matrices/eigen.py @@ -0,0 +1,85 @@ +from sympy.core.symbol import Dummy + +from ..agca.extensions import FiniteExtension +from ..factortools import dup_factor_list +from ..polyroots import roots +from ..polytools import Poly +from ..rootoftools import CRootOf + +from .domainmatrix import DomainMatrix + + +def dom_eigenvects(A, l=Dummy('lambda')): + charpoly = A.charpoly() + rows, cols = A.shape + domain = A.domain + _, factors = dup_factor_list(charpoly, domain) + + rational_eigenvects = [] + algebraic_eigenvects = [] + for base, exp in factors: + if len(base) == 2: + field = domain + eigenval = -base[1] / base[0] + + EE_items = [ + [eigenval if i == j else field.zero for j in range(cols)] + for i in range(rows)] + EE = DomainMatrix(EE_items, (rows, cols), field) + + basis = (A - EE).nullspace() + rational_eigenvects.append((field, eigenval, exp, basis)) + else: + minpoly = Poly.from_list(base, l, domain=domain) + field = FiniteExtension(minpoly) + eigenval = field(l) + + AA_items = [ + [Poly.from_list([item], l, domain=domain).rep for item in row] + for row in A.rep] + AA_items = [[field(item) for item in row] for row in AA_items] + AA = DomainMatrix(AA_items, (rows, cols), field) + EE_items = [ + [eigenval if i == j else field.zero for j in range(cols)] + for i in range(rows)] + EE = DomainMatrix(EE_items, (rows, cols), field) + + basis = (AA - EE).nullspace() + algebraic_eigenvects.append((field, minpoly, exp, basis)) + + return rational_eigenvects, algebraic_eigenvects + + +def dom_eigenvects_to_sympy( + rational_eigenvects, algebraic_eigenvects, + Matrix, **kwargs +): + result = [] + + for field, eigenvalue, multiplicity, eigenvects in rational_eigenvects: + eigenvects = eigenvects.rep + eigenvalue = field.to_sympy(eigenvalue) + new_eigenvects = [ + Matrix([field.to_sympy(x) for x in vect]) + for vect in eigenvects] + result.append((eigenvalue, multiplicity, new_eigenvects)) + + for field, minpoly, multiplicity, eigenvects in algebraic_eigenvects: + eigenvects = eigenvects.rep + l = minpoly.gens[0] + + eigenvects = [[field.to_sympy(x) for x in vect] for vect in eigenvects] + + degree = minpoly.degree() + minpoly = minpoly.as_expr() + eigenvals = roots(minpoly, l, **kwargs) + if len(eigenvals) != degree: + eigenvals = [CRootOf(minpoly, l, idx) for idx in range(degree)] + + for eigenvalue in eigenvals: + new_eigenvects = [ + Matrix([x.subs(l, eigenvalue) for x in vect]) + for vect in eigenvects] + result.append((eigenvalue, multiplicity, new_eigenvects)) + + return result diff --git a/sympy/polys/matrices/exceptions.py b/sympy/polys/matrices/exceptions.py new file mode 100644 index 000000000000..dc13dccb5d67 --- /dev/null +++ b/sympy/polys/matrices/exceptions.py @@ -0,0 +1,29 @@ +from sympy.matrices.common import (NonInvertibleMatrixError, + NonSquareMatrixError, ShapeError) + + +class DDMError(Exception): + """Base class for errors raised by DDM""" + pass + + +class DDMBadInputError(DDMError): + """list of lists is inconsistent with shape""" + pass + + +class DDMDomainError(DDMError): + """domains do not match""" + pass + + +class DDMShapeError(DDMError): + """shapes are inconsistent""" + pass + + +__all__ = [ + 'DDMError', 'DDMShapeError', 'DDMDomainError', + + 'NonSquareMatrixError', 'NonInvertibleMatrixError', 'ShapeError', +] diff --git a/sympy/polys/matrices/tests/__init__.py b/sympy/polys/matrices/tests/__init__.py new file mode 100644 index 000000000000..e69de29bb2d1 diff --git a/sympy/polys/matrices/tests/test_ddm.py b/sympy/polys/matrices/tests/test_ddm.py new file mode 100644 index 000000000000..f2700aa357d1 --- /dev/null +++ b/sympy/polys/matrices/tests/test_ddm.py @@ -0,0 +1,341 @@ +from sympy.testing.pytest import raises +from sympy.core.compatibility import HAS_GMPY + +from sympy.polys import ZZ, QQ + +from sympy.polys.matrices.ddm import DDM +from sympy.polys.matrices.exceptions import ( + DDMShapeError, NonInvertibleMatrixError, DDMDomainError, + DDMBadInputError) + + +def test_DDM_init(): + items = [[ZZ(0), ZZ(1), ZZ(2)], [ZZ(3), ZZ(4), ZZ(5)]] + shape = (2, 3) + ddm = DDM(items, shape, ZZ) + assert ddm.shape == shape + assert ddm.rows == 2 + assert ddm.cols == 3 + assert ddm.domain == ZZ + + raises(DDMBadInputError, lambda: DDM([[ZZ(2), ZZ(3)]], (2, 2), ZZ)) + raises(DDMBadInputError, lambda: DDM([[ZZ(1)], [ZZ(2), ZZ(3)]], (2, 2), ZZ)) + + +def test_DDM_getsetitem(): + ddm = DDM([[ZZ(2), ZZ(3)], [ZZ(4), ZZ(5)]], (2, 2), ZZ) + + assert ddm[0][0] == ZZ(2) + assert ddm[0][1] == ZZ(3) + assert ddm[1][0] == ZZ(4) + assert ddm[1][1] == ZZ(5) + + raises(IndexError, lambda: ddm[2][0]) + raises(IndexError, lambda: ddm[0][2]) + + ddm[0][0] = ZZ(-1) + assert ddm[0][0] == ZZ(-1) + + +def test_DDM_str(): + ddm = DDM([[ZZ(0), ZZ(1)], [ZZ(2), ZZ(3)]], (2, 2), ZZ) + if HAS_GMPY: + assert str(ddm) == 'DDM([[mpz(0), mpz(1)], [mpz(2), mpz(3)]], (2, 2), ZZ)' + else: + assert str(ddm) == 'DDM([[0, 1], [2, 3]], (2, 2), ZZ)' + + +def test_DDM_eq(): + items = [[ZZ(0), ZZ(1)], [ZZ(2), ZZ(3)]] + ddm1 = DDM(items, (2, 2), ZZ) + ddm2 = DDM(items, (2, 2), ZZ) + + assert (ddm1 == ddm1) is True + assert (ddm1 == items) is False + assert (items == ddm1) is False + assert (ddm1 == ddm2) is True + assert (ddm2 == ddm1) is True + + assert (ddm1 != ddm1) is False + assert (ddm1 != items) is True + assert (items != ddm1) is True + assert (ddm1 != ddm2) is False + assert (ddm2 != ddm1) is False + + ddm3 = DDM([[ZZ(0), ZZ(1)], [ZZ(3), ZZ(3)]], (2, 2), ZZ) + ddm3 = DDM(items, (2, 2), QQ) + + assert (ddm1 == ddm3) is False + assert (ddm3 == ddm1) is False + assert (ddm1 != ddm3) is True + assert (ddm3 != ddm1) is True + + +def test_DDM_zeros(): + ddmz = DDM.zeros((3, 4), QQ) + assert list(ddmz) == [[QQ(0)] * 4] * 3 + assert ddmz.shape == (3, 4) + assert ddmz.domain == QQ + + +def test_DDM_eye(): + ddmz = DDM.eye(3, QQ) + f = lambda i, j: QQ(1) if i == j else QQ(0) + assert list(ddmz) == [[f(i, j) for i in range(3)] for j in range(3)] + assert ddmz.shape == (3, 3) + assert ddmz.domain == QQ + + +def test_DDM_copy(): + ddm1 = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) + ddm2 = ddm1.copy() + assert (ddm1 == ddm2) is True + ddm1[0][0] = QQ(-1) + assert (ddm1 == ddm2) is False + ddm2[0][0] = QQ(-1) + assert (ddm1 == ddm2) is True + + +def test_DDM_add(): + A = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) + B = DDM([[ZZ(3)], [ZZ(4)]], (2, 1), ZZ) + C = DDM([[ZZ(4)], [ZZ(6)]], (2, 1), ZZ) + AQ = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) + assert A + B == A.add(B) == C + + raises(DDMShapeError, lambda: A + DDM([[ZZ(5)]], (1, 1), ZZ)) + raises(TypeError, lambda: A + ZZ(1)) + raises(TypeError, lambda: ZZ(1) + A) + raises(DDMDomainError, lambda: A + AQ) + raises(DDMDomainError, lambda: AQ + A) + + +def test_DDM_sub(): + A = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) + B = DDM([[ZZ(3)], [ZZ(4)]], (2, 1), ZZ) + C = DDM([[ZZ(-2)], [ZZ(-2)]], (2, 1), ZZ) + AQ = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) + D = DDM([[ZZ(5)]], (1, 1), ZZ) + assert A - B == A.sub(B) == C + + raises(TypeError, lambda: A - ZZ(1)) + raises(TypeError, lambda: ZZ(1) - A) + raises(DDMShapeError, lambda: A - D) + raises(DDMShapeError, lambda: D - A) + raises(DDMShapeError, lambda: A.sub(D)) + raises(DDMShapeError, lambda: D.sub(A)) + raises(DDMDomainError, lambda: A - AQ) + raises(DDMDomainError, lambda: AQ - A) + raises(DDMDomainError, lambda: A.sub(AQ)) + raises(DDMDomainError, lambda: AQ.sub(A)) + + +def test_DDM_neg(): + A = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) + An = DDM([[ZZ(-1)], [ZZ(-2)]], (2, 1), ZZ) + assert -A == A.neg() == An + assert -An == An.neg() == A + + +def test_DDM_mul(): + A = DDM([[ZZ(1)]], (1, 1), ZZ) + raises(TypeError, lambda: [[1]] * A) + raises(TypeError, lambda: A * [[1]]) + + +def test_DDM_matmul(): + A = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) + B = DDM([[ZZ(3), ZZ(4)]], (1, 2), ZZ) + AB = DDM([[ZZ(3), ZZ(4)], [ZZ(6), ZZ(8)]], (2, 2), ZZ) + BA = DDM([[ZZ(11)]], (1, 1), ZZ) + + assert A @ B == A.matmul(B) == AB + assert B @ A == B.matmul(A) == BA + + raises(TypeError, lambda: A @ 1) + raises(TypeError, lambda: A @ [[3, 4]]) + + Bq = DDM([[QQ(3), QQ(4)]], (1, 2), QQ) + + raises(DDMDomainError, lambda: A @ Bq) + raises(DDMDomainError, lambda: Bq @ A) + + C = DDM([[ZZ(1)]], (1, 1), ZZ) + + assert A @ C == A.matmul(C) == A + + raises(DDMShapeError, lambda: C @ A) + raises(DDMShapeError, lambda: C.matmul(A)) + + Z04 = DDM([], (0, 4), ZZ) + Z40 = DDM([[]]*4, (4, 0), ZZ) + Z50 = DDM([[]]*5, (5, 0), ZZ) + Z05 = DDM([], (0, 5), ZZ) + Z45 = DDM([[0] * 5] * 4, (4, 5), ZZ) + Z54 = DDM([[0] * 4] * 5, (5, 4), ZZ) + Z00 = DDM([], (0, 0), ZZ) + + assert Z04 @ Z45 == Z04.matmul(Z45) == Z05 + assert Z45 @ Z50 == Z45.matmul(Z50) == Z40 + assert Z00 @ Z04 == Z00.matmul(Z04) == Z04 + assert Z50 @ Z00 == Z50.matmul(Z00) == Z50 + assert Z00 @ Z00 == Z00.matmul(Z00) == Z00 + assert Z50 @ Z04 == Z50.matmul(Z04) == Z54 + + raises(DDMShapeError, lambda: Z05 @ Z40) + raises(DDMShapeError, lambda: Z05.matmul(Z40)) + + +def test_DDM_rref(): + + A = DDM([], (0, 4), QQ) + assert A.rref() == (A, []) + + A = DDM([[QQ(0), QQ(1)], [QQ(1), QQ(1)]], (2, 2), QQ) + Ar = DDM([[QQ(1), QQ(0)], [QQ(0), QQ(1)]], (2, 2), QQ) + pivots = [0, 1] + assert A.rref() == (Ar, pivots) + + A = DDM([[QQ(1), QQ(2), QQ(1)], [QQ(3), QQ(4), QQ(1)]], (2, 3), QQ) + Ar = DDM([[QQ(1), QQ(0), QQ(-1)], [QQ(0), QQ(1), QQ(1)]], (2, 3), QQ) + pivots = [0, 1] + assert A.rref() == (Ar, pivots) + + A = DDM([[QQ(3), QQ(4), QQ(1)], [QQ(1), QQ(2), QQ(1)]], (2, 3), QQ) + Ar = DDM([[QQ(1), QQ(0), QQ(-1)], [QQ(0), QQ(1), QQ(1)]], (2, 3), QQ) + pivots = [0, 1] + assert A.rref() == (Ar, pivots) + + A = DDM([[QQ(1), QQ(0)], [QQ(1), QQ(3)], [QQ(0), QQ(1)]], (3, 2), QQ) + Ar = DDM([[QQ(1), QQ(0)], [QQ(0), QQ(1)], [QQ(0), QQ(0)]], (3, 2), QQ) + pivots = [0, 1] + assert A.rref() == (Ar, pivots) + + A = DDM([[QQ(1), QQ(0), QQ(1)], [QQ(3), QQ(0), QQ(1)]], (2, 3), QQ) + Ar = DDM([[QQ(1), QQ(0), QQ(0)], [QQ(0), QQ(0), QQ(1)]], (2, 3), QQ) + pivots = [0, 2] + assert A.rref() == (Ar, pivots) + + +def test_DDM_det(): + # 0x0 case + A = DDM([], (0, 0), ZZ) + assert A.det() == ZZ(1) + + # 1x1 case + A = DDM([[ZZ(2)]], (1, 1), ZZ) + assert A.det() == ZZ(2) + + # 2x2 case + A = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + assert A.det() == ZZ(-2) + + # 3x3 with swap + A = DDM([[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(2), ZZ(5)]], (3, 3), ZZ) + assert A.det() == ZZ(0) + + # 2x2 QQ case + A = DDM([[QQ(1, 2), QQ(1, 2)], [QQ(1, 3), QQ(1, 4)]], (2, 2), QQ) + assert A.det() == QQ(-1, 24) + + # Nonsquare error + A = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) + raises(DDMShapeError, lambda: A.det()) + + # Nonsquare error with empty matrix + A = DDM([], (0, 1), ZZ) + raises(DDMShapeError, lambda: A.det()) + + +def test_DDM_inv(): + A = DDM([[QQ(1, 1), QQ(2, 1)], [QQ(3, 1), QQ(4, 1)]], (2, 2), QQ) + Ainv = DDM([[QQ(-2, 1), QQ(1, 1)], [QQ(3, 2), QQ(-1, 2)]], (2, 2), QQ) + assert A.inv() == Ainv + + A = DDM([[QQ(1), QQ(2)]], (1, 2), QQ) + raises(DDMShapeError, lambda: A.inv()) + + A = DDM([[ZZ(2)]], (1, 1), ZZ) + raises(ValueError, lambda: A.inv()) + + A = DDM([], (0, 0), QQ) + assert A.inv() == A + + A = DDM([[QQ(1), QQ(2)], [QQ(2), QQ(4)]], (2, 2), QQ) + raises(NonInvertibleMatrixError, lambda: A.inv()) + + +def test_DDM_lu(): + A = DDM([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + L, U, swaps = A.lu() + assert L == DDM([[QQ(1), QQ(0)], [QQ(3), QQ(1)]], (2, 2), QQ) + assert U == DDM([[QQ(1), QQ(2)], [QQ(0), QQ(-2)]], (2, 2), QQ) + assert swaps == [] + + A = [[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 1, 2]] + Lexp = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 1, 1]] + Uexp = [[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]] + to_dom = lambda rows, dom: [[dom(e) for e in row] for row in rows] + A = DDM(to_dom(A, QQ), (4, 4), QQ) + Lexp = DDM(to_dom(Lexp, QQ), (4, 4), QQ) + Uexp = DDM(to_dom(Uexp, QQ), (4, 4), QQ) + L, U, swaps = A.lu() + assert L == Lexp + assert U == Uexp + assert swaps == [] + + +def test_DDM_lu_solve(): + # Basic example + A = DDM([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) + x = DDM([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) + assert A.lu_solve(b) == x + + # Example with swaps + A = DDM([[QQ(0), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + assert A.lu_solve(b) == x + + # Overdetermined, consistent + A = DDM([[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]], (3, 2), QQ) + b = DDM([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ) + assert A.lu_solve(b) == x + + # Overdetermined, inconsistent + b = DDM([[QQ(1)], [QQ(2)], [QQ(4)]], (3, 1), QQ) + raises(NonInvertibleMatrixError, lambda: A.lu_solve(b)) + + # Square, noninvertible + A = DDM([[QQ(1), QQ(2)], [QQ(1), QQ(2)]], (2, 2), QQ) + b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) + raises(NonInvertibleMatrixError, lambda: A.lu_solve(b)) + + # Underdetermined + A = DDM([[QQ(1), QQ(2)]], (1, 2), QQ) + b = DDM([[QQ(3)]], (1, 1), QQ) + raises(NotImplementedError, lambda: A.lu_solve(b)) + + # Domain mismatch + bz = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) + raises(DDMDomainError, lambda: A.lu_solve(bz)) + + # Shape mismatch + b3 = DDM([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ) + raises(DDMShapeError, lambda: A.lu_solve(b3)) + + +def test_DDM_charpoly(): + A = DDM([], (0, 0), ZZ) + assert A.charpoly() == [ZZ(1)] + + A = DDM([ + [ZZ(1), ZZ(2), ZZ(3)], + [ZZ(4), ZZ(5), ZZ(6)], + [ZZ(7), ZZ(8), ZZ(9)]], (3, 3), ZZ) + Avec = [ZZ(1), ZZ(-15), ZZ(-18), ZZ(0)] + assert A.charpoly() == Avec + + A = DDM([[ZZ(1), ZZ(2)]], (1, 2), ZZ) + raises(DDMShapeError, lambda: A.charpoly()) + + diff --git a/sympy/polys/matrices/tests/test_dense.py b/sympy/polys/matrices/tests/test_dense.py new file mode 100644 index 000000000000..16f21b622e5d --- /dev/null +++ b/sympy/polys/matrices/tests/test_dense.py @@ -0,0 +1,341 @@ +from sympy.testing.pytest import raises + +from sympy.polys import ZZ, QQ + +from sympy.polys.matrices.ddm import DDM +from sympy.polys.matrices.dense import ( + ddm_iadd, ddm_isub, ddm_ineg, ddm_imatmul, ddm_imul, ddm_irref, + ddm_idet, ddm_iinv, ddm_ilu, ddm_ilu_split, ddm_ilu_solve, ddm_berk) +from sympy.polys.matrices.exceptions import ( + DDMShapeError, NonInvertibleMatrixError, NonSquareMatrixError) + + +def test_ddm_iadd(): + a = [[1, 2], [3, 4]] + b = [[5, 6], [7, 8]] + ddm_iadd(a, b) + assert a == [[6, 8], [10, 12]] + + +def test_ddm_isub(): + a = [[1, 2], [3, 4]] + b = [[5, 6], [7, 8]] + ddm_isub(a, b) + assert a == [[-4, -4], [-4, -4]] + + +def test_ddm_ineg(): + a = [[1, 2], [3, 4]] + ddm_ineg(a) + assert a == [[-1, -2], [-3, -4]] + + +def test_ddm_matmul(): + a = [[1, 2], [3, 4]] + ddm_imul(a, 2) + assert a == [[2, 4], [6, 8]] + + a = [[1, 2], [3, 4]] + ddm_imul(a, 0) + assert a == [[0, 0], [0, 0]] + + +def test_ddm_imatmul(): + a = [[1, 2, 3], [4, 5, 6]] + b = [[1, 2], [3, 4], [5, 6]] + + c1 = [[0, 0], [0, 0]] + ddm_imatmul(c1, a, b) + assert c1 == [[22, 28], [49, 64]] + + c2 = [[0, 0, 0], [0, 0, 0], [0, 0, 0]] + ddm_imatmul(c2, b, a) + assert c2 == [[9, 12, 15], [19, 26, 33], [29, 40, 51]] + + b3 = [[1], [2], [3]] + c3 = [[0], [0]] + ddm_imatmul(c3, a, b3) + assert c3 == [[14], [32]] + + +def test_ddm_irref(): + # Empty matrix + A = [] + Ar = [] + pivots = [] + assert ddm_irref(A) == pivots + assert A == Ar + + # Standard square case + A = [[QQ(0), QQ(1)], [QQ(1), QQ(1)]] + Ar = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]] + pivots = [0, 1] + assert ddm_irref(A) == pivots + assert A == Ar + + # m < n case + A = [[QQ(1), QQ(2), QQ(1)], [QQ(3), QQ(4), QQ(1)]] + Ar = [[QQ(1), QQ(0), QQ(-1)], [QQ(0), QQ(1), QQ(1)]] + pivots = [0, 1] + assert ddm_irref(A) == pivots + assert A == Ar + + # same m < n but reversed + A = [[QQ(3), QQ(4), QQ(1)], [QQ(1), QQ(2), QQ(1)]] + Ar = [[QQ(1), QQ(0), QQ(-1)], [QQ(0), QQ(1), QQ(1)]] + pivots = [0, 1] + assert ddm_irref(A) == pivots + assert A == Ar + + # m > n case + A = [[QQ(1), QQ(0)], [QQ(1), QQ(3)], [QQ(0), QQ(1)]] + Ar = [[QQ(1), QQ(0)], [QQ(0), QQ(1)], [QQ(0), QQ(0)]] + pivots = [0, 1] + assert ddm_irref(A) == pivots + assert A == Ar + + # Example with missing pivot + A = [[QQ(1), QQ(0), QQ(1)], [QQ(3), QQ(0), QQ(1)]] + Ar = [[QQ(1), QQ(0), QQ(0)], [QQ(0), QQ(0), QQ(1)]] + pivots = [0, 2] + assert ddm_irref(A) == pivots + assert A == Ar + + # Example with missing pivot and no replacement + A = [[QQ(0), QQ(1)], [QQ(0), QQ(2)], [QQ(1), QQ(0)]] + Ar = [[QQ(1), QQ(0)], [QQ(0), QQ(1)], [QQ(0), QQ(0)]] + pivots = [0, 1] + assert ddm_irref(A) == pivots + assert A == Ar + + +def test_ddm_idet(): + A = [] + assert ddm_idet(A, ZZ) == ZZ(1) + + A = [[ZZ(2)]] + assert ddm_idet(A, ZZ) == ZZ(2) + + A = [[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]] + assert ddm_idet(A, ZZ) == ZZ(-2) + + A = [[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(3), ZZ(5)]] + assert ddm_idet(A, ZZ) == ZZ(-1) + + A = [[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(2), ZZ(5)]] + assert ddm_idet(A, ZZ) == ZZ(0) + + A = [[QQ(1, 2), QQ(1, 2)], [QQ(1, 3), QQ(1, 4)]] + assert ddm_idet(A, QQ) == QQ(-1, 24) + + +def test_ddm_inv(): + A = [] + Ainv = [] + ddm_iinv(Ainv, A, QQ) + assert Ainv == A + + A = [] + Ainv = [] + raises(ValueError, lambda: ddm_iinv(Ainv, A, ZZ)) + + A = [[QQ(1), QQ(2)]] + Ainv = [[QQ(0), QQ(0)]] + raises(NonSquareMatrixError, lambda: ddm_iinv(Ainv, A, QQ)) + + A = [[QQ(1, 1), QQ(2, 1)], [QQ(3, 1), QQ(4, 1)]] + Ainv = [[QQ(0), QQ(0)], [QQ(0), QQ(0)]] + Ainv_expected = [[QQ(-2, 1), QQ(1, 1)], [QQ(3, 2), QQ(-1, 2)]] + ddm_iinv(Ainv, A, QQ) + assert Ainv == Ainv_expected + + A = [[QQ(1, 1), QQ(2, 1)], [QQ(2, 1), QQ(4, 1)]] + Ainv = [[QQ(0), QQ(0)], [QQ(0), QQ(0)]] + raises(NonInvertibleMatrixError, lambda: ddm_iinv(Ainv, A, QQ)) + + +def test_ddm_ilu(): + A = [] + Alu = [] + swaps = ddm_ilu(A) + assert A == Alu + assert swaps == [] + + A = [[]] + Alu = [[]] + swaps = ddm_ilu(A) + assert A == Alu + assert swaps == [] + + A = [[QQ(1), QQ(2)], [QQ(3), QQ(4)]] + Alu = [[QQ(1), QQ(2)], [QQ(3), QQ(-2)]] + swaps = ddm_ilu(A) + assert A == Alu + assert swaps == [] + + A = [[QQ(0), QQ(2)], [QQ(3), QQ(4)]] + Alu = [[QQ(3), QQ(4)], [QQ(0), QQ(2)]] + swaps = ddm_ilu(A) + assert A == Alu + assert swaps == [(0, 1)] + + A = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)], [QQ(7), QQ(8), QQ(9)]] + Alu = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(-3), QQ(-6)], [QQ(7), QQ(2), QQ(0)]] + swaps = ddm_ilu(A) + assert A == Alu + assert swaps == [] + + A = [[QQ(0), QQ(1), QQ(2)], [QQ(0), QQ(1), QQ(3)], [QQ(1), QQ(1), QQ(2)]] + Alu = [[QQ(1), QQ(1), QQ(2)], [QQ(0), QQ(1), QQ(3)], [QQ(0), QQ(1), QQ(-1)]] + swaps = ddm_ilu(A) + assert A == Alu + assert swaps == [(0, 2)] + + A = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)]] + Alu = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(-3), QQ(-6)]] + swaps = ddm_ilu(A) + assert A == Alu + assert swaps == [] + + A = [[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]] + Alu = [[QQ(1), QQ(2)], [QQ(3), QQ(-2)], [QQ(5), QQ(2)]] + swaps = ddm_ilu(A) + assert A == Alu + assert swaps == [] + + +def test_ddm_ilu_split(): + U = [] + L = [] + Uexp = [] + Lexp = [] + swaps = ddm_ilu_split(L, U, QQ) + assert U == Uexp + assert L == Lexp + assert swaps == [] + + U = [[]] + L = [[QQ(1)]] + Uexp = [[]] + Lexp = [[QQ(1)]] + swaps = ddm_ilu_split(L, U, QQ) + assert U == Uexp + assert L == Lexp + assert swaps == [] + + U = [[QQ(1), QQ(2)], [QQ(3), QQ(4)]] + L = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]] + Uexp = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)]] + Lexp = [[QQ(1), QQ(0)], [QQ(3), QQ(1)]] + swaps = ddm_ilu_split(L, U, QQ) + assert U == Uexp + assert L == Lexp + assert swaps == [] + + U = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)]] + L = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]] + Uexp = [[QQ(1), QQ(2), QQ(3)], [QQ(0), QQ(-3), QQ(-6)]] + Lexp = [[QQ(1), QQ(0)], [QQ(4), QQ(1)]] + swaps = ddm_ilu_split(L, U, QQ) + assert U == Uexp + assert L == Lexp + assert swaps == [] + + U = [[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]] + L = [[QQ(1), QQ(0), QQ(0)], [QQ(0), QQ(1), QQ(0)], [QQ(0), QQ(0), QQ(1)]] + Uexp = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)], [QQ(0), QQ(0)]] + Lexp = [[QQ(1), QQ(0), QQ(0)], [QQ(3), QQ(1), QQ(0)], [QQ(5), QQ(2), QQ(1)]] + swaps = ddm_ilu_split(L, U, QQ) + assert U == Uexp + assert L == Lexp + assert swaps == [] + + +def test_ddm_ilu_solve(): + # Basic example + # A = [[QQ(1), QQ(2)], [QQ(3), QQ(4)]] + U = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)]] + L = [[QQ(1), QQ(0)], [QQ(3), QQ(1)]] + swaps = [] + b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) + x = DDM([[QQ(0)], [QQ(0)]], (2, 1), QQ) + xexp = DDM([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) + ddm_ilu_solve(x, L, U, swaps, b) + assert x == xexp + + # Example with swaps + # A = [[QQ(0), QQ(2)], [QQ(3), QQ(4)]] + U = [[QQ(3), QQ(4)], [QQ(0), QQ(2)]] + L = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]] + swaps = [(0, 1)] + b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) + x = DDM([[QQ(0)], [QQ(0)]], (2, 1), QQ) + xexp = DDM([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) + ddm_ilu_solve(x, L, U, swaps, b) + assert x == xexp + + # Overdetermined, consistent + # A = DDM([[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]], (3, 2), QQ) + U = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)], [QQ(0), QQ(0)]] + L = [[QQ(1), QQ(0), QQ(0)], [QQ(3), QQ(1), QQ(0)], [QQ(5), QQ(2), QQ(1)]] + swaps = [] + b = DDM([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ) + x = DDM([[QQ(0)], [QQ(0)]], (2, 1), QQ) + xexp = DDM([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) + ddm_ilu_solve(x, L, U, swaps, b) + assert x == xexp + + # Overdetermined, inconsistent + b = DDM([[QQ(1)], [QQ(2)], [QQ(4)]], (3, 1), QQ) + raises(NonInvertibleMatrixError, lambda: ddm_ilu_solve(x, L, U, swaps, b)) + + # Square, noninvertible + # A = DDM([[QQ(1), QQ(2)], [QQ(1), QQ(2)]], (2, 2), QQ) + U = [[QQ(1), QQ(2)], [QQ(0), QQ(0)]] + L = [[QQ(1), QQ(0)], [QQ(1), QQ(1)]] + swaps = [] + b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) + raises(NonInvertibleMatrixError, lambda: ddm_ilu_solve(x, L, U, swaps, b)) + + # Underdetermined + # A = DDM([[QQ(1), QQ(2)]], (1, 2), QQ) + U = [[QQ(1), QQ(2)]] + L = [[QQ(1)]] + swaps = [] + b = DDM([[QQ(3)]], (1, 1), QQ) + raises(NotImplementedError, lambda: ddm_ilu_solve(x, L, U, swaps, b)) + + # Shape mismatch + b3 = DDM([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ) + raises(DDMShapeError, lambda: ddm_ilu_solve(x, L, U, swaps, b3)) + + # Empty shape mismatch + U = [[QQ(1)]] + L = [[QQ(1)]] + swaps = [] + x = [[QQ(1)]] + b = [] + raises(DDMShapeError, lambda: ddm_ilu_solve(x, L, U, swaps, b)) + + # Empty system + U = [] + L = [] + swaps = [] + b = [] + x = [] + ddm_ilu_solve(x, L, U, swaps, b) + assert x == [] + + +def test_ddm_charpoly(): + A = [] + assert ddm_berk(A, ZZ) == [[ZZ(1)]] + + A = [[ZZ(1), ZZ(2), ZZ(3)], [ZZ(4), ZZ(5), ZZ(6)], [ZZ(7), ZZ(8), ZZ(9)]] + Avec = [[ZZ(1)], [ZZ(-15)], [ZZ(-18)], [ZZ(0)]] + assert ddm_berk(A, ZZ) == Avec + + A = DDM([[ZZ(1), ZZ(2)]], (1, 2), ZZ) + raises(DDMShapeError, lambda: ddm_berk(A, ZZ)) + + diff --git a/sympy/polys/matrices/tests/test_domainmatrix.py b/sympy/polys/matrices/tests/test_domainmatrix.py new file mode 100644 index 000000000000..375e5ea606c8 --- /dev/null +++ b/sympy/polys/matrices/tests/test_domainmatrix.py @@ -0,0 +1,422 @@ +from sympy.testing.pytest import raises + +from sympy.core.numbers import Rational +from sympy.functions import sqrt + +from sympy.matrices.common import (NonInvertibleMatrixError, + NonSquareMatrixError, ShapeError) +from sympy.matrices.dense import Matrix +from sympy.polys import ZZ, QQ + +from sympy.polys.matrices.domainmatrix import DomainMatrix +from sympy.polys.matrices.exceptions import DDMBadInputError, DDMDomainError +from sympy.polys.matrices.ddm import DDM + + +def test_DomainMatrix_init(): + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + assert A.rep == DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + assert A.shape == (2, 2) + assert A.domain == ZZ + + raises(DDMBadInputError, lambda: DomainMatrix([[ZZ(1), ZZ(2)]], (2, 2), ZZ)) + + +def test_DomainMatrix_from_ddm(): + ddm = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + A = DomainMatrix.from_ddm(ddm) + assert A.rep == ddm + assert A.shape == (2, 2) + assert A.domain == ZZ + + +def test_DomainMatrix_from_list_sympy(): + ddm = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + A = DomainMatrix.from_list_sympy(2, 2, [[1, 2], [3, 4]]) + assert A.rep == ddm + assert A.shape == (2, 2) + assert A.domain == ZZ + + K = QQ.algebraic_field(sqrt(2)) + ddm = DDM( + [[K.convert(1 + sqrt(2)), K.convert(2 + sqrt(2))], + [K.convert(3 + sqrt(2)), K.convert(4 + sqrt(2))]], + (2, 2), + K + ) + A = DomainMatrix.from_list_sympy( + 2, 2, [[1 + sqrt(2), 2 + sqrt(2)], [3 + sqrt(2), 4 + sqrt(2)]], + extension=True) + assert A.rep == ddm + assert A.shape == (2, 2) + assert A.domain == K + + +def test_DomainMatrix_from_Matrix(): + ddm = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + A = DomainMatrix.from_Matrix(Matrix([[1, 2], [3, 4]])) + assert A.rep == ddm + assert A.shape == (2, 2) + assert A.domain == ZZ + + K = QQ.algebraic_field(sqrt(2)) + ddm = DDM( + [[K.convert(1 + sqrt(2)), K.convert(2 + sqrt(2))], + [K.convert(3 + sqrt(2)), K.convert(4 + sqrt(2))]], + (2, 2), + K + ) + A = DomainMatrix.from_Matrix( + Matrix([[1 + sqrt(2), 2 + sqrt(2)], [3 + sqrt(2), 4 + sqrt(2)]]), + extension=True) + assert A.rep == ddm + assert A.shape == (2, 2) + assert A.domain == K + + +def test_DomainMatrix_eq(): + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + assert A == A + B = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(1)]], (2, 2), ZZ) + assert A != B + C = [[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]] + assert A != C + + +def test_DomainMatrix_get_domain(): + K, items = DomainMatrix.get_domain([1, 2, 3, 4]) + assert items == [ZZ(1), ZZ(2), ZZ(3), ZZ(4)] + assert K == ZZ + + K, items = DomainMatrix.get_domain([1, 2, 3, Rational(1, 2)]) + assert items == [QQ(1), QQ(2), QQ(3), QQ(1, 2)] + assert K == QQ + + +def test_DomainMatrix_convert_to(): + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + Aq = A.convert_to(QQ) + assert Aq == DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + + +def test_DomainMatrix_to_field(): + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + Aq = A.to_field() + assert Aq == DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + + +def test_DomainMatrix_unify(): + Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + Aq = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + assert Az.unify(Az) == (Az, Az) + assert Az.unify(Aq) == (Aq, Aq) + assert Aq.unify(Az) == (Aq, Aq) + assert Aq.unify(Aq) == (Aq, Aq) + + +def test_DomainMatrix_to_Matrix(): + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + assert A.to_Matrix() == Matrix([[1, 2], [3, 4]]) + + +def test_DomainMatrix_repr(): + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + assert repr(A) == 'DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ)' + + +def test_DomainMatrix_add(): + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + B = DomainMatrix([[ZZ(2), ZZ(4)], [ZZ(6), ZZ(8)]], (2, 2), ZZ) + assert A + A == A.add(A) == B + + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + L = [[2, 3], [3, 4]] + raises(TypeError, lambda: A + L) + raises(TypeError, lambda: L + A) + + A1 = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + A2 = DomainMatrix([[ZZ(1), ZZ(2)]], (1, 2), ZZ) + raises(ShapeError, lambda: A1 + A2) + raises(ShapeError, lambda: A2 + A1) + raises(ShapeError, lambda: A1.add(A2)) + raises(ShapeError, lambda: A2.add(A1)) + + Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + Aq = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + raises(ValueError, lambda: Az + Aq) + raises(ValueError, lambda: Aq + Az) + raises(ValueError, lambda: Az.add(Aq)) + raises(ValueError, lambda: Aq.add(Az)) + + +def test_DomainMatrix_sub(): + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + B = DomainMatrix([[ZZ(0), ZZ(0)], [ZZ(0), ZZ(0)]], (2, 2), ZZ) + assert A - A == A.sub(A) == B + + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + L = [[2, 3], [3, 4]] + raises(TypeError, lambda: A - L) + raises(TypeError, lambda: L - A) + + A1 = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + A2 = DomainMatrix([[ZZ(1), ZZ(2)]], (1, 2), ZZ) + raises(ShapeError, lambda: A1 - A2) + raises(ShapeError, lambda: A2 - A1) + raises(ShapeError, lambda: A1.sub(A2)) + raises(ShapeError, lambda: A2.sub(A1)) + + Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + Aq = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + raises(ValueError, lambda: Az - Aq) + raises(ValueError, lambda: Aq - Az) + raises(ValueError, lambda: Az.sub(Aq)) + raises(ValueError, lambda: Aq.sub(Az)) + + +def test_DomainMatrix_neg(): + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + Aneg = DomainMatrix([[ZZ(-1), ZZ(-2)], [ZZ(-3), ZZ(-4)]], (2, 2), ZZ) + assert -A == A.neg() == Aneg + + +def test_DomainMatrix_mul(): + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + A2 = DomainMatrix([[ZZ(7), ZZ(10)], [ZZ(15), ZZ(22)]], (2, 2), ZZ) + assert A*A == A.matmul(A) == A2 + + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + L = [[1, 2], [3, 4]] + raises(TypeError, lambda: A * L) + raises(TypeError, lambda: L * A) + + Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + Aq = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + raises(DDMDomainError, lambda: Az * Aq) + raises(DDMDomainError, lambda: Aq * Az) + raises(DDMDomainError, lambda: Az.matmul(Aq)) + raises(DDMDomainError, lambda: Aq.matmul(Az)) + + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + AA = DomainMatrix([[ZZ(2), ZZ(4)], [ZZ(6), ZZ(8)]], (2, 2), ZZ) + x = ZZ(2) + assert A * x == x * A == A.mul(x) == AA + + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + AA = DomainMatrix([[ZZ(0), ZZ(0)], [ZZ(0), ZZ(0)]], (2, 2), ZZ) + x = ZZ(0) + assert A * x == x * A == A.mul(x) == AA + + +def test_DomainMatrix_pow(): + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + A2 = DomainMatrix([[ZZ(7), ZZ(10)], [ZZ(15), ZZ(22)]], (2, 2), ZZ) + A3 = DomainMatrix([[ZZ(37), ZZ(54)], [ZZ(81), ZZ(118)]], (2, 2), ZZ) + eye = DomainMatrix([[ZZ(1), ZZ(0)], [ZZ(0), ZZ(1)]], (2, 2), ZZ) + assert A**0 == A.pow(0) == eye + assert A**1 == A.pow(1) == A + assert A**2 == A.pow(2) == A2 + assert A**3 == A.pow(3) == A3 + + raises(TypeError, lambda: A ** Rational(1, 2)) + raises(NotImplementedError, lambda: A ** -1) + raises(NotImplementedError, lambda: A.pow(-1)) + + +def test_DomainMatrix_rref(): + A = DomainMatrix([], (0, 1), QQ) + assert A.rref() == (A, ()) + + A = DomainMatrix([[QQ(1)]], (1, 1), QQ) + assert A.rref() == (A, (0,)) + + A = DomainMatrix([[QQ(0)]], (1, 1), QQ) + assert A.rref() == (A, ()) + + A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + Ar, pivots = A.rref() + assert Ar == DomainMatrix([[QQ(1), QQ(0)], [QQ(0), QQ(1)]], (2, 2), QQ) + assert pivots == (0, 1) + + A = DomainMatrix([[QQ(0), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + Ar, pivots = A.rref() + assert Ar == DomainMatrix([[QQ(1), QQ(0)], [QQ(0), QQ(1)]], (2, 2), QQ) + assert pivots == (0, 1) + + A = DomainMatrix([[QQ(0), QQ(2)], [QQ(0), QQ(4)]], (2, 2), QQ) + Ar, pivots = A.rref() + assert Ar == DomainMatrix([[QQ(0), QQ(1)], [QQ(0), QQ(0)]], (2, 2), QQ) + assert pivots == (1,) + + Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + raises(ValueError, lambda: Az.rref()) + + +def test_DomainMatrix_inv(): + A = DomainMatrix([], (0, 0), QQ) + assert A.inv() == A + + A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + Ainv = DomainMatrix([[QQ(-2), QQ(1)], [QQ(3, 2), QQ(-1, 2)]], (2, 2), QQ) + assert A.inv() == Ainv + + Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + raises(ValueError, lambda: Az.inv()) + + Ans = DomainMatrix([[QQ(1), QQ(2)]], (1, 2), QQ) + raises(NonSquareMatrixError, lambda: Ans.inv()) + + Aninv = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(6)]], (2, 2), QQ) + raises(NonInvertibleMatrixError, lambda: Aninv.inv()) + + +def test_DomainMatrix_det(): + A = DomainMatrix([], (0, 0), ZZ) + assert A.det() == 1 + + A = DomainMatrix([[1]], (1, 1), ZZ) + assert A.det() == 1 + + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + assert A.det() == ZZ(-2) + + A = DomainMatrix([[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(3), ZZ(5)]], (3, 3), ZZ) + assert A.det() == ZZ(-1) + + A = DomainMatrix([[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(2), ZZ(5)]], (3, 3), ZZ) + assert A.det() == ZZ(0) + + Ans = DomainMatrix([[QQ(1), QQ(2)]], (1, 2), QQ) + raises(NonSquareMatrixError, lambda: Ans.det()) + + A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + assert A.det() == QQ(-2) + + +def test_DomainMatrix_lu(): + A = DomainMatrix([], (0, 0), QQ) + assert A.lu() == (A, A, []) + + A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + L = DomainMatrix([[QQ(1), QQ(0)], [QQ(3), QQ(1)]], (2, 2), QQ) + U = DomainMatrix([[QQ(1), QQ(2)], [QQ(0), QQ(-2)]], (2, 2), QQ) + swaps = [] + assert A.lu() == (L, U, swaps) + + A = DomainMatrix([[QQ(0), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + L = DomainMatrix([[QQ(1), QQ(0)], [QQ(0), QQ(1)]], (2, 2), QQ) + U = DomainMatrix([[QQ(3), QQ(4)], [QQ(0), QQ(2)]], (2, 2), QQ) + swaps = [(0, 1)] + assert A.lu() == (L, U, swaps) + + A = DomainMatrix([[QQ(1), QQ(2)], [QQ(2), QQ(4)]], (2, 2), QQ) + L = DomainMatrix([[QQ(1), QQ(0)], [QQ(2), QQ(1)]], (2, 2), QQ) + U = DomainMatrix([[QQ(1), QQ(2)], [QQ(0), QQ(0)]], (2, 2), QQ) + swaps = [] + assert A.lu() == (L, U, swaps) + + A = DomainMatrix([[QQ(0), QQ(2)], [QQ(0), QQ(4)]], (2, 2), QQ) + L = DomainMatrix([[QQ(1), QQ(0)], [QQ(0), QQ(1)]], (2, 2), QQ) + U = DomainMatrix([[QQ(0), QQ(2)], [QQ(0), QQ(4)]], (2, 2), QQ) + swaps = [] + assert A.lu() == (L, U, swaps) + + A = DomainMatrix([[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)]], (2, 3), QQ) + L = DomainMatrix([[QQ(1), QQ(0)], [QQ(4), QQ(1)]], (2, 2), QQ) + U = DomainMatrix([[QQ(1), QQ(2), QQ(3)], [QQ(0), QQ(-3), QQ(-6)]], (2, 3), QQ) + swaps = [] + assert A.lu() == (L, U, swaps) + + A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]], (3, 2), QQ) + L = DomainMatrix([ + [QQ(1), QQ(0), QQ(0)], + [QQ(3), QQ(1), QQ(0)], + [QQ(5), QQ(2), QQ(1)]], (3, 3), QQ) + U = DomainMatrix([[QQ(1), QQ(2)], [QQ(0), QQ(-2)], [QQ(0), QQ(0)]], (3, 2), QQ) + swaps = [] + assert A.lu() == (L, U, swaps) + + A = [[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 1, 2]] + L = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 1, 1]] + U = [[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]] + to_dom = lambda rows, dom: [[dom(e) for e in row] for row in rows] + A = DomainMatrix(to_dom(A, QQ), (4, 4), QQ) + L = DomainMatrix(to_dom(L, QQ), (4, 4), QQ) + U = DomainMatrix(to_dom(U, QQ), (4, 4), QQ) + assert A.lu() == (L, U, []) + + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + raises(ValueError, lambda: A.lu()) + + +def test_DomainMatrix_lu_solve(): + # Base case + A = b = x = DomainMatrix([], (0, 0), QQ) + assert A.lu_solve(b) == x + + # Basic example + A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + b = DomainMatrix([[QQ(1)], [QQ(2)]], (2, 1), QQ) + x = DomainMatrix([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) + assert A.lu_solve(b) == x + + # Example with swaps + A = DomainMatrix([[QQ(0), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + b = DomainMatrix([[QQ(1)], [QQ(2)]], (2, 1), QQ) + x = DomainMatrix([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) + assert A.lu_solve(b) == x + + # Non-invertible + A = DomainMatrix([[QQ(1), QQ(2)], [QQ(2), QQ(4)]], (2, 2), QQ) + b = DomainMatrix([[QQ(1)], [QQ(2)]], (2, 1), QQ) + raises(NonInvertibleMatrixError, lambda: A.lu_solve(b)) + + # Overdetermined, consistent + A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]], (3, 2), QQ) + b = DomainMatrix([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ) + x = DomainMatrix([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) + assert A.lu_solve(b) == x + + # Overdetermined, inconsistent + A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]], (3, 2), QQ) + b = DomainMatrix([[QQ(1)], [QQ(2)], [QQ(4)]], (3, 1), QQ) + raises(NonInvertibleMatrixError, lambda: A.lu_solve(b)) + + # Underdetermined + A = DomainMatrix([[QQ(1), QQ(2)]], (1, 2), QQ) + b = DomainMatrix([[QQ(1)]], (1, 1), QQ) + raises(NotImplementedError, lambda: A.lu_solve(b)) + + # Non-field + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + b = DomainMatrix([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) + raises(ValueError, lambda: A.lu_solve(b)) + + # Shape mismatch + A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) + b = DomainMatrix([[QQ(1), QQ(2)]], (1, 2), QQ) + raises(ShapeError, lambda: A.lu_solve(b)) + + +def test_DomainMatrix_charpoly(): + A = DomainMatrix([], (0, 0), ZZ) + assert A.charpoly() == [ZZ(1)] + + A = DomainMatrix([[1]], (1, 1), ZZ) + assert A.charpoly() == [ZZ(1), ZZ(-1)] + + A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) + assert A.charpoly() == [ZZ(1), ZZ(-5), ZZ(-2)] + + A = DomainMatrix([[ZZ(1), ZZ(2), ZZ(3)], [ZZ(4), ZZ(5), ZZ(6)], [ZZ(7), ZZ(8), ZZ(9)]], (3, 3), ZZ) + assert A.charpoly() == [ZZ(1), ZZ(-15), ZZ(-18), ZZ(0)] + + Ans = DomainMatrix([[QQ(1), QQ(2)]], (1, 2), QQ) + raises(NonSquareMatrixError, lambda: Ans.charpoly()) + + +def test_DomainMatrix_eye(): + A = DomainMatrix.eye(3, QQ) + assert A.rep == DDM.eye(3, QQ) + assert A.shape == (3, 3) + assert A.domain == QQ diff --git a/sympy/polys/tests/test_domainmatrix.py b/sympy/polys/tests/test_domainmatrix.py deleted file mode 100644 index ca17bc4c4462..000000000000 --- a/sympy/polys/tests/test_domainmatrix.py +++ /dev/null @@ -1,1085 +0,0 @@ -from sympy.core.compatibility import HAS_GMPY -from sympy.core.numbers import Rational -from sympy.functions import sqrt -from sympy.matrices.common import (NonInvertibleMatrixError, - NonSquareMatrixError, ShapeError) -from sympy.matrices.dense import Matrix -from sympy.polys import ZZ, QQ -from sympy.polys.domainmatrix import ( - DDM, - DDMBadInputError, DDMDomainError, DDMShapeError, - DomainMatrix, - ddm_iadd, ddm_isub, ddm_ineg, ddm_imatmul, ddm_imul, - ddm_irref, ddm_idet, ddm_iinv, ddm_ilu, ddm_ilu_split, ddm_ilu_solve, - ddm_berk, - ) - -from sympy.testing.pytest import raises - - -def test_DDM_init(): - items = [[ZZ(0), ZZ(1), ZZ(2)], [ZZ(3), ZZ(4), ZZ(5)]] - shape = (2, 3) - ddm = DDM(items, shape, ZZ) - assert ddm.shape == shape - assert ddm.rows == 2 - assert ddm.cols == 3 - assert ddm.domain == ZZ - - raises(DDMBadInputError, lambda: DDM([[ZZ(2), ZZ(3)]], (2, 2), ZZ)) - raises(DDMBadInputError, lambda: DDM([[ZZ(1)], [ZZ(2), ZZ(3)]], (2, 2), ZZ)) - - -def test_DDM_getsetitem(): - ddm = DDM([[ZZ(2), ZZ(3)], [ZZ(4), ZZ(5)]], (2, 2), ZZ) - - assert ddm[0][0] == ZZ(2) - assert ddm[0][1] == ZZ(3) - assert ddm[1][0] == ZZ(4) - assert ddm[1][1] == ZZ(5) - - raises(IndexError, lambda: ddm[2][0]) - raises(IndexError, lambda: ddm[0][2]) - - ddm[0][0] = ZZ(-1) - assert ddm[0][0] == ZZ(-1) - - -def test_DDM_str(): - ddm = DDM([[ZZ(0), ZZ(1)], [ZZ(2), ZZ(3)]], (2, 2), ZZ) - if HAS_GMPY: - assert str(ddm) == 'DDM([[mpz(0), mpz(1)], [mpz(2), mpz(3)]], (2, 2), ZZ)' - else: - assert str(ddm) == 'DDM([[0, 1], [2, 3]], (2, 2), ZZ)' - - -def test_DDM_eq(): - items = [[ZZ(0), ZZ(1)], [ZZ(2), ZZ(3)]] - ddm1 = DDM(items, (2, 2), ZZ) - ddm2 = DDM(items, (2, 2), ZZ) - - assert (ddm1 == ddm1) is True - assert (ddm1 == items) is False - assert (items == ddm1) is False - assert (ddm1 == ddm2) is True - assert (ddm2 == ddm1) is True - - assert (ddm1 != ddm1) is False - assert (ddm1 != items) is True - assert (items != ddm1) is True - assert (ddm1 != ddm2) is False - assert (ddm2 != ddm1) is False - - ddm3 = DDM([[ZZ(0), ZZ(1)], [ZZ(3), ZZ(3)]], (2, 2), ZZ) - ddm3 = DDM(items, (2, 2), QQ) - - assert (ddm1 == ddm3) is False - assert (ddm3 == ddm1) is False - assert (ddm1 != ddm3) is True - assert (ddm3 != ddm1) is True - - -def test_DDM_zeros(): - ddmz = DDM.zeros((3, 4), QQ) - assert list(ddmz) == [[QQ(0)] * 4] * 3 - assert ddmz.shape == (3, 4) - assert ddmz.domain == QQ - - -def test_DDM_eye(): - ddmz = DDM.eye(3, QQ) - f = lambda i, j: QQ(1) if i == j else QQ(0) - assert list(ddmz) == [[f(i, j) for i in range(3)] for j in range(3)] - assert ddmz.shape == (3, 3) - assert ddmz.domain == QQ - - -def test_DDM_copy(): - ddm1 = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) - ddm2 = ddm1.copy() - assert (ddm1 == ddm2) is True - ddm1[0][0] = QQ(-1) - assert (ddm1 == ddm2) is False - ddm2[0][0] = QQ(-1) - assert (ddm1 == ddm2) is True - - -def test_DDM_add(): - A = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) - B = DDM([[ZZ(3)], [ZZ(4)]], (2, 1), ZZ) - C = DDM([[ZZ(4)], [ZZ(6)]], (2, 1), ZZ) - AQ = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) - assert A + B == A.add(B) == C - - raises(DDMShapeError, lambda: A + DDM([[ZZ(5)]], (1, 1), ZZ)) - raises(TypeError, lambda: A + ZZ(1)) - raises(TypeError, lambda: ZZ(1) + A) - raises(DDMDomainError, lambda: A + AQ) - raises(DDMDomainError, lambda: AQ + A) - - -def test_DDM_sub(): - A = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) - B = DDM([[ZZ(3)], [ZZ(4)]], (2, 1), ZZ) - C = DDM([[ZZ(-2)], [ZZ(-2)]], (2, 1), ZZ) - AQ = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) - D = DDM([[ZZ(5)]], (1, 1), ZZ) - assert A - B == A.sub(B) == C - - raises(TypeError, lambda: A - ZZ(1)) - raises(TypeError, lambda: ZZ(1) - A) - raises(DDMShapeError, lambda: A - D) - raises(DDMShapeError, lambda: D - A) - raises(DDMShapeError, lambda: A.sub(D)) - raises(DDMShapeError, lambda: D.sub(A)) - raises(DDMDomainError, lambda: A - AQ) - raises(DDMDomainError, lambda: AQ - A) - raises(DDMDomainError, lambda: A.sub(AQ)) - raises(DDMDomainError, lambda: AQ.sub(A)) - - -def test_DDM_neg(): - A = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) - An = DDM([[ZZ(-1)], [ZZ(-2)]], (2, 1), ZZ) - assert -A == A.neg() == An - assert -An == An.neg() == A - - -def test_DDM_mul(): - A = DDM([[ZZ(1)]], (1, 1), ZZ) - raises(TypeError, lambda: [[1]] * A) - raises(TypeError, lambda: A * [[1]]) - - -def test_DDM_matmul(): - A = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) - B = DDM([[ZZ(3), ZZ(4)]], (1, 2), ZZ) - AB = DDM([[ZZ(3), ZZ(4)], [ZZ(6), ZZ(8)]], (2, 2), ZZ) - BA = DDM([[ZZ(11)]], (1, 1), ZZ) - - assert A @ B == A.matmul(B) == AB - assert B @ A == B.matmul(A) == BA - - raises(TypeError, lambda: A @ 1) - raises(TypeError, lambda: A @ [[3, 4]]) - - Bq = DDM([[QQ(3), QQ(4)]], (1, 2), QQ) - - raises(DDMDomainError, lambda: A @ Bq) - raises(DDMDomainError, lambda: Bq @ A) - - C = DDM([[ZZ(1)]], (1, 1), ZZ) - - assert A @ C == A.matmul(C) == A - - raises(DDMShapeError, lambda: C @ A) - raises(DDMShapeError, lambda: C.matmul(A)) - - Z04 = DDM([], (0, 4), ZZ) - Z40 = DDM([[]]*4, (4, 0), ZZ) - Z50 = DDM([[]]*5, (5, 0), ZZ) - Z05 = DDM([], (0, 5), ZZ) - Z45 = DDM([[0] * 5] * 4, (4, 5), ZZ) - Z54 = DDM([[0] * 4] * 5, (5, 4), ZZ) - Z00 = DDM([], (0, 0), ZZ) - - assert Z04 @ Z45 == Z04.matmul(Z45) == Z05 - assert Z45 @ Z50 == Z45.matmul(Z50) == Z40 - assert Z00 @ Z04 == Z00.matmul(Z04) == Z04 - assert Z50 @ Z00 == Z50.matmul(Z00) == Z50 - assert Z00 @ Z00 == Z00.matmul(Z00) == Z00 - assert Z50 @ Z04 == Z50.matmul(Z04) == Z54 - - raises(DDMShapeError, lambda: Z05 @ Z40) - raises(DDMShapeError, lambda: Z05.matmul(Z40)) - - -def test_DDM_rref(): - - A = DDM([], (0, 4), QQ) - assert A.rref() == (A, []) - - A = DDM([[QQ(0), QQ(1)], [QQ(1), QQ(1)]], (2, 2), QQ) - Ar = DDM([[QQ(1), QQ(0)], [QQ(0), QQ(1)]], (2, 2), QQ) - pivots = [0, 1] - assert A.rref() == (Ar, pivots) - - A = DDM([[QQ(1), QQ(2), QQ(1)], [QQ(3), QQ(4), QQ(1)]], (2, 3), QQ) - Ar = DDM([[QQ(1), QQ(0), QQ(-1)], [QQ(0), QQ(1), QQ(1)]], (2, 3), QQ) - pivots = [0, 1] - assert A.rref() == (Ar, pivots) - - A = DDM([[QQ(3), QQ(4), QQ(1)], [QQ(1), QQ(2), QQ(1)]], (2, 3), QQ) - Ar = DDM([[QQ(1), QQ(0), QQ(-1)], [QQ(0), QQ(1), QQ(1)]], (2, 3), QQ) - pivots = [0, 1] - assert A.rref() == (Ar, pivots) - - A = DDM([[QQ(1), QQ(0)], [QQ(1), QQ(3)], [QQ(0), QQ(1)]], (3, 2), QQ) - Ar = DDM([[QQ(1), QQ(0)], [QQ(0), QQ(1)], [QQ(0), QQ(0)]], (3, 2), QQ) - pivots = [0, 1] - assert A.rref() == (Ar, pivots) - - A = DDM([[QQ(1), QQ(0), QQ(1)], [QQ(3), QQ(0), QQ(1)]], (2, 3), QQ) - Ar = DDM([[QQ(1), QQ(0), QQ(0)], [QQ(0), QQ(0), QQ(1)]], (2, 3), QQ) - pivots = [0, 2] - assert A.rref() == (Ar, pivots) - - -def test_DDM_det(): - # 0x0 case - A = DDM([], (0, 0), ZZ) - assert A.det() == ZZ(1) - - # 1x1 case - A = DDM([[ZZ(2)]], (1, 1), ZZ) - assert A.det() == ZZ(2) - - # 2x2 case - A = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - assert A.det() == ZZ(-2) - - # 3x3 with swap - A = DDM([[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(2), ZZ(5)]], (3, 3), ZZ) - assert A.det() == ZZ(0) - - # 2x2 QQ case - A = DDM([[QQ(1, 2), QQ(1, 2)], [QQ(1, 3), QQ(1, 4)]], (2, 2), QQ) - assert A.det() == QQ(-1, 24) - - # Nonsquare error - A = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) - raises(DDMShapeError, lambda: A.det()) - - # Nonsquare error with empty matrix - A = DDM([], (0, 1), ZZ) - raises(DDMShapeError, lambda: A.det()) - - -def test_DDM_inv(): - A = DDM([[QQ(1, 1), QQ(2, 1)], [QQ(3, 1), QQ(4, 1)]], (2, 2), QQ) - Ainv = DDM([[QQ(-2, 1), QQ(1, 1)], [QQ(3, 2), QQ(-1, 2)]], (2, 2), QQ) - assert A.inv() == Ainv - - A = DDM([[QQ(1), QQ(2)]], (1, 2), QQ) - raises(DDMShapeError, lambda: A.inv()) - - A = DDM([[ZZ(2)]], (1, 1), ZZ) - raises(ValueError, lambda: A.inv()) - - A = DDM([], (0, 0), QQ) - assert A.inv() == A - - A = DDM([[QQ(1), QQ(2)], [QQ(2), QQ(4)]], (2, 2), QQ) - raises(NonInvertibleMatrixError, lambda: A.inv()) - - -def test_DDM_lu(): - A = DDM([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - L, U, swaps = A.lu() - assert L == DDM([[QQ(1), QQ(0)], [QQ(3), QQ(1)]], (2, 2), QQ) - assert U == DDM([[QQ(1), QQ(2)], [QQ(0), QQ(-2)]], (2, 2), QQ) - assert swaps == [] - - A = [[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 1, 2]] - Lexp = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 1, 1]] - Uexp = [[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]] - to_dom = lambda rows, dom: [[dom(e) for e in row] for row in rows] - A = DDM(to_dom(A, QQ), (4, 4), QQ) - Lexp = DDM(to_dom(Lexp, QQ), (4, 4), QQ) - Uexp = DDM(to_dom(Uexp, QQ), (4, 4), QQ) - L, U, swaps = A.lu() - assert L == Lexp - assert U == Uexp - assert swaps == [] - - -def test_DDM_lu_solve(): - # Basic example - A = DDM([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) - x = DDM([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) - assert A.lu_solve(b) == x - - # Example with swaps - A = DDM([[QQ(0), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - assert A.lu_solve(b) == x - - # Overdetermined, consistent - A = DDM([[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]], (3, 2), QQ) - b = DDM([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ) - assert A.lu_solve(b) == x - - # Overdetermined, inconsistent - b = DDM([[QQ(1)], [QQ(2)], [QQ(4)]], (3, 1), QQ) - raises(NonInvertibleMatrixError, lambda: A.lu_solve(b)) - - # Square, noninvertible - A = DDM([[QQ(1), QQ(2)], [QQ(1), QQ(2)]], (2, 2), QQ) - b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) - raises(NonInvertibleMatrixError, lambda: A.lu_solve(b)) - - # Underdetermined - A = DDM([[QQ(1), QQ(2)]], (1, 2), QQ) - b = DDM([[QQ(3)]], (1, 1), QQ) - raises(NotImplementedError, lambda: A.lu_solve(b)) - - # Domain mismatch - bz = DDM([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) - raises(DDMDomainError, lambda: A.lu_solve(bz)) - - # Shape mismatch - b3 = DDM([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ) - raises(DDMShapeError, lambda: A.lu_solve(b3)) - - -def test_DDM_charpoly(): - A = DDM([], (0, 0), ZZ) - assert A.charpoly() == [ZZ(1)] - - A = DDM([ - [ZZ(1), ZZ(2), ZZ(3)], - [ZZ(4), ZZ(5), ZZ(6)], - [ZZ(7), ZZ(8), ZZ(9)]], (3, 3), ZZ) - Avec = [ZZ(1), ZZ(-15), ZZ(-18), ZZ(0)] - assert A.charpoly() == Avec - - A = DDM([[ZZ(1), ZZ(2)]], (1, 2), ZZ) - raises(DDMShapeError, lambda: A.charpoly()) - - -def test_ddm_iadd(): - a = [[1, 2], [3, 4]] - b = [[5, 6], [7, 8]] - ddm_iadd(a, b) - assert a == [[6, 8], [10, 12]] - - -def test_ddm_isub(): - a = [[1, 2], [3, 4]] - b = [[5, 6], [7, 8]] - ddm_isub(a, b) - assert a == [[-4, -4], [-4, -4]] - - -def test_ddm_ineg(): - a = [[1, 2], [3, 4]] - ddm_ineg(a) - assert a == [[-1, -2], [-3, -4]] - - -def test_ddm_matmul(): - a = [[1, 2], [3, 4]] - ddm_imul(a, 2) - assert a == [[2, 4], [6, 8]] - - a = [[1, 2], [3, 4]] - ddm_imul(a, 0) - assert a == [[0, 0], [0, 0]] - - -def test_ddm_imatmul(): - a = [[1, 2, 3], [4, 5, 6]] - b = [[1, 2], [3, 4], [5, 6]] - - c1 = [[0, 0], [0, 0]] - ddm_imatmul(c1, a, b) - assert c1 == [[22, 28], [49, 64]] - - c2 = [[0, 0, 0], [0, 0, 0], [0, 0, 0]] - ddm_imatmul(c2, b, a) - assert c2 == [[9, 12, 15], [19, 26, 33], [29, 40, 51]] - - b3 = [[1], [2], [3]] - c3 = [[0], [0]] - ddm_imatmul(c3, a, b3) - assert c3 == [[14], [32]] - - -def test_ddm_irref(): - # Empty matrix - A = [] - Ar = [] - pivots = [] - assert ddm_irref(A) == pivots - assert A == Ar - - # Standard square case - A = [[QQ(0), QQ(1)], [QQ(1), QQ(1)]] - Ar = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]] - pivots = [0, 1] - assert ddm_irref(A) == pivots - assert A == Ar - - # m < n case - A = [[QQ(1), QQ(2), QQ(1)], [QQ(3), QQ(4), QQ(1)]] - Ar = [[QQ(1), QQ(0), QQ(-1)], [QQ(0), QQ(1), QQ(1)]] - pivots = [0, 1] - assert ddm_irref(A) == pivots - assert A == Ar - - # same m < n but reversed - A = [[QQ(3), QQ(4), QQ(1)], [QQ(1), QQ(2), QQ(1)]] - Ar = [[QQ(1), QQ(0), QQ(-1)], [QQ(0), QQ(1), QQ(1)]] - pivots = [0, 1] - assert ddm_irref(A) == pivots - assert A == Ar - - # m > n case - A = [[QQ(1), QQ(0)], [QQ(1), QQ(3)], [QQ(0), QQ(1)]] - Ar = [[QQ(1), QQ(0)], [QQ(0), QQ(1)], [QQ(0), QQ(0)]] - pivots = [0, 1] - assert ddm_irref(A) == pivots - assert A == Ar - - # Example with missing pivot - A = [[QQ(1), QQ(0), QQ(1)], [QQ(3), QQ(0), QQ(1)]] - Ar = [[QQ(1), QQ(0), QQ(0)], [QQ(0), QQ(0), QQ(1)]] - pivots = [0, 2] - assert ddm_irref(A) == pivots - assert A == Ar - - # Example with missing pivot and no replacement - A = [[QQ(0), QQ(1)], [QQ(0), QQ(2)], [QQ(1), QQ(0)]] - Ar = [[QQ(1), QQ(0)], [QQ(0), QQ(1)], [QQ(0), QQ(0)]] - pivots = [0, 1] - assert ddm_irref(A) == pivots - assert A == Ar - - -def test_ddm_idet(): - A = [] - assert ddm_idet(A, ZZ) == ZZ(1) - - A = [[ZZ(2)]] - assert ddm_idet(A, ZZ) == ZZ(2) - - A = [[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]] - assert ddm_idet(A, ZZ) == ZZ(-2) - - A = [[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(3), ZZ(5)]] - assert ddm_idet(A, ZZ) == ZZ(-1) - - A = [[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(2), ZZ(5)]] - assert ddm_idet(A, ZZ) == ZZ(0) - - A = [[QQ(1, 2), QQ(1, 2)], [QQ(1, 3), QQ(1, 4)]] - assert ddm_idet(A, QQ) == QQ(-1, 24) - - -def test_ddm_inv(): - A = [] - Ainv = [] - ddm_iinv(Ainv, A, QQ) - assert Ainv == A - - A = [] - Ainv = [] - raises(ValueError, lambda: ddm_iinv(Ainv, A, ZZ)) - - A = [[QQ(1), QQ(2)]] - Ainv = [[QQ(0), QQ(0)]] - raises(NonSquareMatrixError, lambda: ddm_iinv(Ainv, A, QQ)) - - A = [[QQ(1, 1), QQ(2, 1)], [QQ(3, 1), QQ(4, 1)]] - Ainv = [[QQ(0), QQ(0)], [QQ(0), QQ(0)]] - Ainv_expected = [[QQ(-2, 1), QQ(1, 1)], [QQ(3, 2), QQ(-1, 2)]] - ddm_iinv(Ainv, A, QQ) - assert Ainv == Ainv_expected - - A = [[QQ(1, 1), QQ(2, 1)], [QQ(2, 1), QQ(4, 1)]] - Ainv = [[QQ(0), QQ(0)], [QQ(0), QQ(0)]] - raises(NonInvertibleMatrixError, lambda: ddm_iinv(Ainv, A, QQ)) - - -def test_ddm_ilu(): - A = [] - Alu = [] - swaps = ddm_ilu(A) - assert A == Alu - assert swaps == [] - - A = [[]] - Alu = [[]] - swaps = ddm_ilu(A) - assert A == Alu - assert swaps == [] - - A = [[QQ(1), QQ(2)], [QQ(3), QQ(4)]] - Alu = [[QQ(1), QQ(2)], [QQ(3), QQ(-2)]] - swaps = ddm_ilu(A) - assert A == Alu - assert swaps == [] - - A = [[QQ(0), QQ(2)], [QQ(3), QQ(4)]] - Alu = [[QQ(3), QQ(4)], [QQ(0), QQ(2)]] - swaps = ddm_ilu(A) - assert A == Alu - assert swaps == [(0, 1)] - - A = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)], [QQ(7), QQ(8), QQ(9)]] - Alu = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(-3), QQ(-6)], [QQ(7), QQ(2), QQ(0)]] - swaps = ddm_ilu(A) - assert A == Alu - assert swaps == [] - - A = [[QQ(0), QQ(1), QQ(2)], [QQ(0), QQ(1), QQ(3)], [QQ(1), QQ(1), QQ(2)]] - Alu = [[QQ(1), QQ(1), QQ(2)], [QQ(0), QQ(1), QQ(3)], [QQ(0), QQ(1), QQ(-1)]] - swaps = ddm_ilu(A) - assert A == Alu - assert swaps == [(0, 2)] - - A = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)]] - Alu = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(-3), QQ(-6)]] - swaps = ddm_ilu(A) - assert A == Alu - assert swaps == [] - - A = [[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]] - Alu = [[QQ(1), QQ(2)], [QQ(3), QQ(-2)], [QQ(5), QQ(2)]] - swaps = ddm_ilu(A) - assert A == Alu - assert swaps == [] - - -def test_ddm_ilu_split(): - U = [] - L = [] - Uexp = [] - Lexp = [] - swaps = ddm_ilu_split(L, U, QQ) - assert U == Uexp - assert L == Lexp - assert swaps == [] - - U = [[]] - L = [[QQ(1)]] - Uexp = [[]] - Lexp = [[QQ(1)]] - swaps = ddm_ilu_split(L, U, QQ) - assert U == Uexp - assert L == Lexp - assert swaps == [] - - U = [[QQ(1), QQ(2)], [QQ(3), QQ(4)]] - L = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]] - Uexp = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)]] - Lexp = [[QQ(1), QQ(0)], [QQ(3), QQ(1)]] - swaps = ddm_ilu_split(L, U, QQ) - assert U == Uexp - assert L == Lexp - assert swaps == [] - - U = [[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)]] - L = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]] - Uexp = [[QQ(1), QQ(2), QQ(3)], [QQ(0), QQ(-3), QQ(-6)]] - Lexp = [[QQ(1), QQ(0)], [QQ(4), QQ(1)]] - swaps = ddm_ilu_split(L, U, QQ) - assert U == Uexp - assert L == Lexp - assert swaps == [] - - U = [[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]] - L = [[QQ(1), QQ(0), QQ(0)], [QQ(0), QQ(1), QQ(0)], [QQ(0), QQ(0), QQ(1)]] - Uexp = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)], [QQ(0), QQ(0)]] - Lexp = [[QQ(1), QQ(0), QQ(0)], [QQ(3), QQ(1), QQ(0)], [QQ(5), QQ(2), QQ(1)]] - swaps = ddm_ilu_split(L, U, QQ) - assert U == Uexp - assert L == Lexp - assert swaps == [] - - -def test_ddm_ilu_solve(): - # Basic example - # A = [[QQ(1), QQ(2)], [QQ(3), QQ(4)]] - U = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)]] - L = [[QQ(1), QQ(0)], [QQ(3), QQ(1)]] - swaps = [] - b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) - x = DDM([[QQ(0)], [QQ(0)]], (2, 1), QQ) - xexp = DDM([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) - ddm_ilu_solve(x, L, U, swaps, b) - assert x == xexp - - # Example with swaps - # A = [[QQ(0), QQ(2)], [QQ(3), QQ(4)]] - U = [[QQ(3), QQ(4)], [QQ(0), QQ(2)]] - L = [[QQ(1), QQ(0)], [QQ(0), QQ(1)]] - swaps = [(0, 1)] - b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) - x = DDM([[QQ(0)], [QQ(0)]], (2, 1), QQ) - xexp = DDM([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) - ddm_ilu_solve(x, L, U, swaps, b) - assert x == xexp - - # Overdetermined, consistent - # A = DDM([[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]], (3, 2), QQ) - U = [[QQ(1), QQ(2)], [QQ(0), QQ(-2)], [QQ(0), QQ(0)]] - L = [[QQ(1), QQ(0), QQ(0)], [QQ(3), QQ(1), QQ(0)], [QQ(5), QQ(2), QQ(1)]] - swaps = [] - b = DDM([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ) - x = DDM([[QQ(0)], [QQ(0)]], (2, 1), QQ) - xexp = DDM([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) - ddm_ilu_solve(x, L, U, swaps, b) - assert x == xexp - - # Overdetermined, inconsistent - b = DDM([[QQ(1)], [QQ(2)], [QQ(4)]], (3, 1), QQ) - raises(NonInvertibleMatrixError, lambda: ddm_ilu_solve(x, L, U, swaps, b)) - - # Square, noninvertible - # A = DDM([[QQ(1), QQ(2)], [QQ(1), QQ(2)]], (2, 2), QQ) - U = [[QQ(1), QQ(2)], [QQ(0), QQ(0)]] - L = [[QQ(1), QQ(0)], [QQ(1), QQ(1)]] - swaps = [] - b = DDM([[QQ(1)], [QQ(2)]], (2, 1), QQ) - raises(NonInvertibleMatrixError, lambda: ddm_ilu_solve(x, L, U, swaps, b)) - - # Underdetermined - # A = DDM([[QQ(1), QQ(2)]], (1, 2), QQ) - U = [[QQ(1), QQ(2)]] - L = [[QQ(1)]] - swaps = [] - b = DDM([[QQ(3)]], (1, 1), QQ) - raises(NotImplementedError, lambda: ddm_ilu_solve(x, L, U, swaps, b)) - - # Shape mismatch - b3 = DDM([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ) - raises(DDMShapeError, lambda: ddm_ilu_solve(x, L, U, swaps, b3)) - - # Empty shape mismatch - U = [[QQ(1)]] - L = [[QQ(1)]] - swaps = [] - x = [[QQ(1)]] - b = [] - raises(DDMShapeError, lambda: ddm_ilu_solve(x, L, U, swaps, b)) - - # Empty system - U = [] - L = [] - swaps = [] - b = [] - x = [] - ddm_ilu_solve(x, L, U, swaps, b) - assert x == [] - - -def test_ddm_charpoly(): - A = [] - assert ddm_berk(A, ZZ) == [[ZZ(1)]] - - A = [[ZZ(1), ZZ(2), ZZ(3)], [ZZ(4), ZZ(5), ZZ(6)], [ZZ(7), ZZ(8), ZZ(9)]] - Avec = [[ZZ(1)], [ZZ(-15)], [ZZ(-18)], [ZZ(0)]] - assert ddm_berk(A, ZZ) == Avec - - A = DDM([[ZZ(1), ZZ(2)]], (1, 2), ZZ) - raises(DDMShapeError, lambda: ddm_berk(A, ZZ)) - - -def test_DomainMatrix_init(): - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - assert A.rep == DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - assert A.shape == (2, 2) - assert A.domain == ZZ - - raises(DDMBadInputError, lambda: DomainMatrix([[ZZ(1), ZZ(2)]], (2, 2), ZZ)) - - -def test_DomainMatrix_from_ddm(): - ddm = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - A = DomainMatrix.from_ddm(ddm) - assert A.rep == ddm - assert A.shape == (2, 2) - assert A.domain == ZZ - - -def test_DomainMatrix_from_list_sympy(): - ddm = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - A = DomainMatrix.from_list_sympy(2, 2, [[1, 2], [3, 4]]) - assert A.rep == ddm - assert A.shape == (2, 2) - assert A.domain == ZZ - - K = QQ.algebraic_field(sqrt(2)) - ddm = DDM( - [[K.convert(1 + sqrt(2)), K.convert(2 + sqrt(2))], - [K.convert(3 + sqrt(2)), K.convert(4 + sqrt(2))]], - (2, 2), - K - ) - A = DomainMatrix.from_list_sympy( - 2, 2, [[1 + sqrt(2), 2 + sqrt(2)], [3 + sqrt(2), 4 + sqrt(2)]], - extension=True) - assert A.rep == ddm - assert A.shape == (2, 2) - assert A.domain == K - - -def test_DomainMatrix_from_Matrix(): - ddm = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - A = DomainMatrix.from_Matrix(Matrix([[1, 2], [3, 4]])) - assert A.rep == ddm - assert A.shape == (2, 2) - assert A.domain == ZZ - - K = QQ.algebraic_field(sqrt(2)) - ddm = DDM( - [[K.convert(1 + sqrt(2)), K.convert(2 + sqrt(2))], - [K.convert(3 + sqrt(2)), K.convert(4 + sqrt(2))]], - (2, 2), - K - ) - A = DomainMatrix.from_Matrix( - Matrix([[1 + sqrt(2), 2 + sqrt(2)], [3 + sqrt(2), 4 + sqrt(2)]]), - extension=True) - assert A.rep == ddm - assert A.shape == (2, 2) - assert A.domain == K - - -def test_DomainMatrix_eq(): - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - assert A == A - B = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(1)]], (2, 2), ZZ) - assert A != B - C = [[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]] - assert A != C - - -def test_DomainMatrix_get_domain(): - K, items = DomainMatrix.get_domain([1, 2, 3, 4]) - assert items == [ZZ(1), ZZ(2), ZZ(3), ZZ(4)] - assert K == ZZ - - K, items = DomainMatrix.get_domain([1, 2, 3, Rational(1, 2)]) - assert items == [QQ(1), QQ(2), QQ(3), QQ(1, 2)] - assert K == QQ - - -def test_DomainMatrix_convert_to(): - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - Aq = A.convert_to(QQ) - assert Aq == DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - - -def test_DomainMatrix_to_field(): - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - Aq = A.to_field() - assert Aq == DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - - -def test_DomainMatrix_unify(): - Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - Aq = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - assert Az.unify(Az) == (Az, Az) - assert Az.unify(Aq) == (Aq, Aq) - assert Aq.unify(Az) == (Aq, Aq) - assert Aq.unify(Aq) == (Aq, Aq) - - -def test_DomainMatrix_to_Matrix(): - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - assert A.to_Matrix() == Matrix([[1, 2], [3, 4]]) - - -def test_DomainMatrix_repr(): - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - assert repr(A) == 'DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ)' - - -def test_DomainMatrix_add(): - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - B = DomainMatrix([[ZZ(2), ZZ(4)], [ZZ(6), ZZ(8)]], (2, 2), ZZ) - assert A + A == A.add(A) == B - - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - L = [[2, 3], [3, 4]] - raises(TypeError, lambda: A + L) - raises(TypeError, lambda: L + A) - - A1 = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - A2 = DomainMatrix([[ZZ(1), ZZ(2)]], (1, 2), ZZ) - raises(ShapeError, lambda: A1 + A2) - raises(ShapeError, lambda: A2 + A1) - raises(ShapeError, lambda: A1.add(A2)) - raises(ShapeError, lambda: A2.add(A1)) - - Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - Aq = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - raises(ValueError, lambda: Az + Aq) - raises(ValueError, lambda: Aq + Az) - raises(ValueError, lambda: Az.add(Aq)) - raises(ValueError, lambda: Aq.add(Az)) - - -def test_DomainMatrix_sub(): - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - B = DomainMatrix([[ZZ(0), ZZ(0)], [ZZ(0), ZZ(0)]], (2, 2), ZZ) - assert A - A == A.sub(A) == B - - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - L = [[2, 3], [3, 4]] - raises(TypeError, lambda: A - L) - raises(TypeError, lambda: L - A) - - A1 = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - A2 = DomainMatrix([[ZZ(1), ZZ(2)]], (1, 2), ZZ) - raises(ShapeError, lambda: A1 - A2) - raises(ShapeError, lambda: A2 - A1) - raises(ShapeError, lambda: A1.sub(A2)) - raises(ShapeError, lambda: A2.sub(A1)) - - Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - Aq = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - raises(ValueError, lambda: Az - Aq) - raises(ValueError, lambda: Aq - Az) - raises(ValueError, lambda: Az.sub(Aq)) - raises(ValueError, lambda: Aq.sub(Az)) - - -def test_DomainMatrix_neg(): - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - Aneg = DomainMatrix([[ZZ(-1), ZZ(-2)], [ZZ(-3), ZZ(-4)]], (2, 2), ZZ) - assert -A == A.neg() == Aneg - - -def test_DomainMatrix_mul(): - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - A2 = DomainMatrix([[ZZ(7), ZZ(10)], [ZZ(15), ZZ(22)]], (2, 2), ZZ) - assert A*A == A.matmul(A) == A2 - - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - L = [[1, 2], [3, 4]] - raises(TypeError, lambda: A * L) - raises(TypeError, lambda: L * A) - - Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - Aq = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - raises(DDMDomainError, lambda: Az * Aq) - raises(DDMDomainError, lambda: Aq * Az) - raises(DDMDomainError, lambda: Az.matmul(Aq)) - raises(DDMDomainError, lambda: Aq.matmul(Az)) - - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - AA = DomainMatrix([[ZZ(2), ZZ(4)], [ZZ(6), ZZ(8)]], (2, 2), ZZ) - x = ZZ(2) - assert A * x == x * A == A.mul(x) == AA - - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - AA = DomainMatrix([[ZZ(0), ZZ(0)], [ZZ(0), ZZ(0)]], (2, 2), ZZ) - x = ZZ(0) - assert A * x == x * A == A.mul(x) == AA - - -def test_DomainMatrix_pow(): - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - A2 = DomainMatrix([[ZZ(7), ZZ(10)], [ZZ(15), ZZ(22)]], (2, 2), ZZ) - A3 = DomainMatrix([[ZZ(37), ZZ(54)], [ZZ(81), ZZ(118)]], (2, 2), ZZ) - eye = DomainMatrix([[ZZ(1), ZZ(0)], [ZZ(0), ZZ(1)]], (2, 2), ZZ) - assert A**0 == A.pow(0) == eye - assert A**1 == A.pow(1) == A - assert A**2 == A.pow(2) == A2 - assert A**3 == A.pow(3) == A3 - - raises(TypeError, lambda: A ** Rational(1, 2)) - raises(NotImplementedError, lambda: A ** -1) - raises(NotImplementedError, lambda: A.pow(-1)) - - -def test_DomainMatrix_rref(): - A = DomainMatrix([], (0, 1), QQ) - assert A.rref() == (A, ()) - - A = DomainMatrix([[QQ(1)]], (1, 1), QQ) - assert A.rref() == (A, (0,)) - - A = DomainMatrix([[QQ(0)]], (1, 1), QQ) - assert A.rref() == (A, ()) - - A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - Ar, pivots = A.rref() - assert Ar == DomainMatrix([[QQ(1), QQ(0)], [QQ(0), QQ(1)]], (2, 2), QQ) - assert pivots == (0, 1) - - A = DomainMatrix([[QQ(0), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - Ar, pivots = A.rref() - assert Ar == DomainMatrix([[QQ(1), QQ(0)], [QQ(0), QQ(1)]], (2, 2), QQ) - assert pivots == (0, 1) - - A = DomainMatrix([[QQ(0), QQ(2)], [QQ(0), QQ(4)]], (2, 2), QQ) - Ar, pivots = A.rref() - assert Ar == DomainMatrix([[QQ(0), QQ(1)], [QQ(0), QQ(0)]], (2, 2), QQ) - assert pivots == (1,) - - Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - raises(ValueError, lambda: Az.rref()) - - -def test_DomainMatrix_inv(): - A = DomainMatrix([], (0, 0), QQ) - assert A.inv() == A - - A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - Ainv = DomainMatrix([[QQ(-2), QQ(1)], [QQ(3, 2), QQ(-1, 2)]], (2, 2), QQ) - assert A.inv() == Ainv - - Az = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - raises(ValueError, lambda: Az.inv()) - - Ans = DomainMatrix([[QQ(1), QQ(2)]], (1, 2), QQ) - raises(NonSquareMatrixError, lambda: Ans.inv()) - - Aninv = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(6)]], (2, 2), QQ) - raises(NonInvertibleMatrixError, lambda: Aninv.inv()) - - -def test_DomainMatrix_det(): - A = DomainMatrix([], (0, 0), ZZ) - assert A.det() == 1 - - A = DomainMatrix([[1]], (1, 1), ZZ) - assert A.det() == 1 - - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - assert A.det() == ZZ(-2) - - A = DomainMatrix([[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(3), ZZ(5)]], (3, 3), ZZ) - assert A.det() == ZZ(-1) - - A = DomainMatrix([[ZZ(1), ZZ(2), ZZ(3)], [ZZ(1), ZZ(2), ZZ(4)], [ZZ(1), ZZ(2), ZZ(5)]], (3, 3), ZZ) - assert A.det() == ZZ(0) - - Ans = DomainMatrix([[QQ(1), QQ(2)]], (1, 2), QQ) - raises(NonSquareMatrixError, lambda: Ans.det()) - - A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - assert A.det() == QQ(-2) - - -def test_DomainMatrix_lu(): - A = DomainMatrix([], (0, 0), QQ) - assert A.lu() == (A, A, []) - - A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - L = DomainMatrix([[QQ(1), QQ(0)], [QQ(3), QQ(1)]], (2, 2), QQ) - U = DomainMatrix([[QQ(1), QQ(2)], [QQ(0), QQ(-2)]], (2, 2), QQ) - swaps = [] - assert A.lu() == (L, U, swaps) - - A = DomainMatrix([[QQ(0), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - L = DomainMatrix([[QQ(1), QQ(0)], [QQ(0), QQ(1)]], (2, 2), QQ) - U = DomainMatrix([[QQ(3), QQ(4)], [QQ(0), QQ(2)]], (2, 2), QQ) - swaps = [(0, 1)] - assert A.lu() == (L, U, swaps) - - A = DomainMatrix([[QQ(1), QQ(2)], [QQ(2), QQ(4)]], (2, 2), QQ) - L = DomainMatrix([[QQ(1), QQ(0)], [QQ(2), QQ(1)]], (2, 2), QQ) - U = DomainMatrix([[QQ(1), QQ(2)], [QQ(0), QQ(0)]], (2, 2), QQ) - swaps = [] - assert A.lu() == (L, U, swaps) - - A = DomainMatrix([[QQ(0), QQ(2)], [QQ(0), QQ(4)]], (2, 2), QQ) - L = DomainMatrix([[QQ(1), QQ(0)], [QQ(0), QQ(1)]], (2, 2), QQ) - U = DomainMatrix([[QQ(0), QQ(2)], [QQ(0), QQ(4)]], (2, 2), QQ) - swaps = [] - assert A.lu() == (L, U, swaps) - - A = DomainMatrix([[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)]], (2, 3), QQ) - L = DomainMatrix([[QQ(1), QQ(0)], [QQ(4), QQ(1)]], (2, 2), QQ) - U = DomainMatrix([[QQ(1), QQ(2), QQ(3)], [QQ(0), QQ(-3), QQ(-6)]], (2, 3), QQ) - swaps = [] - assert A.lu() == (L, U, swaps) - - A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]], (3, 2), QQ) - L = DomainMatrix([ - [QQ(1), QQ(0), QQ(0)], - [QQ(3), QQ(1), QQ(0)], - [QQ(5), QQ(2), QQ(1)]], (3, 3), QQ) - U = DomainMatrix([[QQ(1), QQ(2)], [QQ(0), QQ(-2)], [QQ(0), QQ(0)]], (3, 2), QQ) - swaps = [] - assert A.lu() == (L, U, swaps) - - A = [[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 1, 2]] - L = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 1, 1]] - U = [[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]] - to_dom = lambda rows, dom: [[dom(e) for e in row] for row in rows] - A = DomainMatrix(to_dom(A, QQ), (4, 4), QQ) - L = DomainMatrix(to_dom(L, QQ), (4, 4), QQ) - U = DomainMatrix(to_dom(U, QQ), (4, 4), QQ) - assert A.lu() == (L, U, []) - - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - raises(ValueError, lambda: A.lu()) - - -def test_DomainMatrix_lu_solve(): - # Base case - A = b = x = DomainMatrix([], (0, 0), QQ) - assert A.lu_solve(b) == x - - # Basic example - A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - b = DomainMatrix([[QQ(1)], [QQ(2)]], (2, 1), QQ) - x = DomainMatrix([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) - assert A.lu_solve(b) == x - - # Example with swaps - A = DomainMatrix([[QQ(0), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - b = DomainMatrix([[QQ(1)], [QQ(2)]], (2, 1), QQ) - x = DomainMatrix([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) - assert A.lu_solve(b) == x - - # Non-invertible - A = DomainMatrix([[QQ(1), QQ(2)], [QQ(2), QQ(4)]], (2, 2), QQ) - b = DomainMatrix([[QQ(1)], [QQ(2)]], (2, 1), QQ) - raises(NonInvertibleMatrixError, lambda: A.lu_solve(b)) - - # Overdetermined, consistent - A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]], (3, 2), QQ) - b = DomainMatrix([[QQ(1)], [QQ(2)], [QQ(3)]], (3, 1), QQ) - x = DomainMatrix([[QQ(0)], [QQ(1, 2)]], (2, 1), QQ) - assert A.lu_solve(b) == x - - # Overdetermined, inconsistent - A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)], [QQ(5), QQ(6)]], (3, 2), QQ) - b = DomainMatrix([[QQ(1)], [QQ(2)], [QQ(4)]], (3, 1), QQ) - raises(NonInvertibleMatrixError, lambda: A.lu_solve(b)) - - # Underdetermined - A = DomainMatrix([[QQ(1), QQ(2)]], (1, 2), QQ) - b = DomainMatrix([[QQ(1)]], (1, 1), QQ) - raises(NotImplementedError, lambda: A.lu_solve(b)) - - # Non-field - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - b = DomainMatrix([[ZZ(1)], [ZZ(2)]], (2, 1), ZZ) - raises(ValueError, lambda: A.lu_solve(b)) - - # Shape mismatch - A = DomainMatrix([[QQ(1), QQ(2)], [QQ(3), QQ(4)]], (2, 2), QQ) - b = DomainMatrix([[QQ(1), QQ(2)]], (1, 2), QQ) - raises(ShapeError, lambda: A.lu_solve(b)) - - -def test_DomainMatrix_charpoly(): - A = DomainMatrix([], (0, 0), ZZ) - assert A.charpoly() == [ZZ(1)] - - A = DomainMatrix([[1]], (1, 1), ZZ) - assert A.charpoly() == [ZZ(1), ZZ(-1)] - - A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ) - assert A.charpoly() == [ZZ(1), ZZ(-5), ZZ(-2)] - - A = DomainMatrix([[ZZ(1), ZZ(2), ZZ(3)], [ZZ(4), ZZ(5), ZZ(6)], [ZZ(7), ZZ(8), ZZ(9)]], (3, 3), ZZ) - assert A.charpoly() == [ZZ(1), ZZ(-15), ZZ(-18), ZZ(0)] - - Ans = DomainMatrix([[QQ(1), QQ(2)]], (1, 2), QQ) - raises(NonSquareMatrixError, lambda: Ans.charpoly()) - - -def test_DomainMatrix_eye(): - A = DomainMatrix.eye(3, QQ) - assert A.rep == DDM.eye(3, QQ) - assert A.shape == (3, 3) - assert A.domain == QQ From 4c2bc06b4fb024ca4c85e0e0a5ca2b81b410d8be Mon Sep 17 00:00:00 2001 From: Oscar Benjamin Date: Thu, 7 Jan 2021 19:53:32 +0000 Subject: [PATCH 2/6] doc(polys): add module docstrings in sympy.polys.matrices Add basic docstrings for each of the modules in sympy.polys.matrices just to give some explanation of what each module is for and the classes that they define. --- sympy/polys/matrices/__init__.py | 10 +++++ sympy/polys/matrices/ddm.py | 63 ++++++++++++++++++++++++++++ sympy/polys/matrices/dense.py | 36 ++++++++++++++++ sympy/polys/matrices/domainmatrix.py | 11 +++++ sympy/polys/matrices/eigen.py | 5 +++ sympy/polys/matrices/exceptions.py | 11 +++++ 6 files changed, 136 insertions(+) diff --git a/sympy/polys/matrices/__init__.py b/sympy/polys/matrices/__init__.py index 1d58d4174344..3099d7efa9a9 100644 --- a/sympy/polys/matrices/__init__.py +++ b/sympy/polys/matrices/__init__.py @@ -1,3 +1,13 @@ +""" + +sympy.polys.matrices package. + +The main export from this package is the DomainMatrix class which is a +lower-level implementation of matrices based on the polys Domains. This +implentation is typically a lot faster than sympy's standard Matrix class but +is a work in progress and is still experimental. + +""" from .domainmatrix import DomainMatrix __all__ = [ diff --git a/sympy/polys/matrices/ddm.py b/sympy/polys/matrices/ddm.py index b0652fea4237..bca4e0bdc150 100644 --- a/sympy/polys/matrices/ddm.py +++ b/sympy/polys/matrices/ddm.py @@ -1,3 +1,66 @@ +""" + +Module for the DDM class. + +The DDM class is an internal representation used by DomainMatrix. The letters +DDM stand for Dense Domain Matrix. A DDM instance represents a matrix using +elements from a polynomial Domain (e.g. ZZ, QQ, ...) in a dense-matrix +representation. + +Basic usage: + + >>> from sympy import ZZ, QQ + >>> from sympy.polys.matrices.ddm import DDM + >>> A = DDM([[ZZ(0), ZZ(1)], [ZZ(-1), ZZ(0)]], (2, 2), ZZ) + >>> A.shape + (2, 2) + >>> A + [[0, 1], [-1, 0]] + >>> type(A) + + >>> A @ A + [[-1, 0], [0, -1]] + +The ddm_* functions are designed to operate on DDM as well as on an ordinary +list of lists: + + >>> from sympy.polys.matrices.dense import ddm_idet + >>> ddm_idet(A, QQ) + 1 + >>> ddm_idet([[0, 1], [-1, 0]], QQ) + 1 + >>> A + [[-1, 0], [0, 1]] + +Note that ddm_idet modifies the input matrix in-place. It is recommended to +use the DDM.det method as a friendlier interface to this instead which takes +care of copying the matrix: + + >>> B = DDM([[ZZ(0), ZZ(1)], [ZZ(-1), ZZ(0)]], (2, 2), ZZ) + >>> B.det() + 1 + +Normally DDM would not be used directly and is just part of the internal +representation of DomainMatrix which adds further functionality including e.g. +unifying domains. + +The dense format used by DDM is a list of lists of elements e.g. the 2x2 +identity matrix is like [[1, 0], [0, 1]]. The DDM class itself is a subclass +of list and its list items are plain lists. Elements are accessed as e.g. +ddm[i][j] where ddm[i] gives the ith row and ddm[i][j] gets the element in the +jth column of that row. Subclassing list makes e.g. iteration and indexing +very efficient. We do not override __getitem__ because it would lose that +benefit. + +The core routines are implemented by the ddm_* functions defined in dense.py. +Those functions are intended to be able to operate on a raw list-of-lists +representation of matrices with most functions operating in-place. The DDM +class takes care of copying etc and also stores a Domain object associated +with its elements. This makes it possible to implement things like A + B with +domain checking and also shape checking so that the list of lists +representation is friendlier. + +""" from .exceptions import DDMBadInputError, DDMShapeError, DDMDomainError from .dense import ( diff --git a/sympy/polys/matrices/dense.py b/sympy/polys/matrices/dense.py index 01d5dfc42fb5..f04b5c81fdd1 100644 --- a/sympy/polys/matrices/dense.py +++ b/sympy/polys/matrices/dense.py @@ -1,3 +1,39 @@ +""" + +Module for the ddm_* routines for operating on a matrix in list of lists +matrix representation. + +These routines are used internally by the DDM class which also provides a +friendlier interface for them. The idea here is to implement core matrix +routines in a way that can be applied to any simple list representation +without the need to use any particular matrix class. For example we can +compute the RREF of a matrix like: + + >>> from sympy.polys.matrices.dense import ddm_irref + >>> M = [[1, 2, 3], [4, 5, 6]] + >>> pivots = ddm_irref(M) + >>> M + [[1.0, 0.0, -1.0], [0, 1.0, 2.0]] + +These are lower-level routines that work mostly in place.The routines at this +level should not need to know what the domain of the elements is but should +ideally document what operations they will use and what functions they need to +be provided with. + +The next-level up is the DDM class which uses these routines but wraps them up +with an interface that handles copying etc and keeps track of the Domain of +the elements of the matrix: + + >>> from sympy.polys.domains import QQ + >>> from sympy.polys.matrices.ddm import DDM + >>> M = DDM([[QQ(1), QQ(2), QQ(3)], [QQ(4), QQ(5), QQ(6)]], (2, 3), QQ) + >>> M + [[1, 2, 3], [4, 5, 6]] + >>> Mrref, pivots = M.rref() + >>> Mrref + [[1, 0, -1], [0, 1, 2]] + +""" from operator import mul from .exceptions import ( diff --git a/sympy/polys/matrices/domainmatrix.py b/sympy/polys/matrices/domainmatrix.py index 808fdd861919..f473807512be 100644 --- a/sympy/polys/matrices/domainmatrix.py +++ b/sympy/polys/matrices/domainmatrix.py @@ -1,3 +1,14 @@ +""" + +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 sympy.core.sympify import _sympify from ..constructor import construct_domain diff --git a/sympy/polys/matrices/eigen.py b/sympy/polys/matrices/eigen.py index 6cbbd4fbf26e..b8f13c0e0e12 100644 --- a/sympy/polys/matrices/eigen.py +++ b/sympy/polys/matrices/eigen.py @@ -1,3 +1,8 @@ +""" + +Routines for computing eigenvectors with DomainMatrix. + +""" from sympy.core.symbol import Dummy from ..agca.extensions import FiniteExtension diff --git a/sympy/polys/matrices/exceptions.py b/sympy/polys/matrices/exceptions.py index dc13dccb5d67..ce22a3b50bd8 100644 --- a/sympy/polys/matrices/exceptions.py +++ b/sympy/polys/matrices/exceptions.py @@ -1,3 +1,14 @@ +""" + +Module to define exceptions to be used in sympy.polys.matrices modules and +classes. + +Ideally all exceptions raised in these modules would be defined and documented +here and not e.g. imported from matrices. Also ideally generic exceptions like +ValueError/TypeError would not be raised anywhere. + +""" + from sympy.matrices.common import (NonInvertibleMatrixError, NonSquareMatrixError, ShapeError) From da196f0c4980d51645929dda0b99579526e1c427 Mon Sep 17 00:00:00 2001 From: Oscar Benjamin Date: Thu, 7 Jan 2021 19:57:20 +0000 Subject: [PATCH 3/6] docs(polys): add docstring to redundant domainmatrix.py --- sympy/polys/domainmatrix.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/sympy/polys/domainmatrix.py b/sympy/polys/domainmatrix.py index 2414264bdcf6..c0ccaaa4cb96 100644 --- a/sympy/polys/domainmatrix.py +++ b/sympy/polys/domainmatrix.py @@ -1,6 +1,12 @@ -# -# Stub to expose DomainMatrix which has now moved to -# sympy.polys.matrices package. It should really be imported from there. -# +""" +Stub module to expose DomainMatrix which has now moved to +sympy.polys.matrices package. It should now be imported as: + + >>> from sympy.polys.matrices import DomainMatrix + +This module might be removed in future. +""" from sympy.polys.matrices.domainmatrix import DomainMatrix + +__all__ = ['DomainMatrix'] From b6af2f3be70b26aa04761657886c0787ca3b45f7 Mon Sep 17 00:00:00 2001 From: Oscar Benjamin Date: Thu, 7 Jan 2021 20:03:39 +0000 Subject: [PATCH 4/6] maint(polys): remove redundant newlines --- sympy/polys/matrices/tests/test_dense.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/sympy/polys/matrices/tests/test_dense.py b/sympy/polys/matrices/tests/test_dense.py index 16f21b622e5d..5061052d0d85 100644 --- a/sympy/polys/matrices/tests/test_dense.py +++ b/sympy/polys/matrices/tests/test_dense.py @@ -337,5 +337,3 @@ def test_ddm_charpoly(): A = DDM([[ZZ(1), ZZ(2)]], (1, 2), ZZ) raises(DDMShapeError, lambda: ddm_berk(A, ZZ)) - - From 429c2fbb2f583ee70164d9d2d7bb087ade395bcc Mon Sep 17 00:00:00 2001 From: Oscar Benjamin Date: Thu, 7 Jan 2021 20:15:34 +0000 Subject: [PATCH 5/6] maint(polys): remove redundant newlines --- sympy/polys/matrices/tests/test_ddm.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/sympy/polys/matrices/tests/test_ddm.py b/sympy/polys/matrices/tests/test_ddm.py index f2700aa357d1..b8c75638241d 100644 --- a/sympy/polys/matrices/tests/test_ddm.py +++ b/sympy/polys/matrices/tests/test_ddm.py @@ -337,5 +337,3 @@ def test_DDM_charpoly(): A = DDM([[ZZ(1), ZZ(2)]], (1, 2), ZZ) raises(DDMShapeError, lambda: A.charpoly()) - - From ac9c802db304e9a500c2244409e49f34ba77cc7e Mon Sep 17 00:00:00 2001 From: Oscar Benjamin Date: Fri, 8 Jan 2021 18:37:33 +0000 Subject: [PATCH 6/6] doc(polys): add docstring in sympy.polys.matrices Fix typo in docstring. --- sympy/polys/matrices/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sympy/polys/matrices/__init__.py b/sympy/polys/matrices/__init__.py index 3099d7efa9a9..704fa9829469 100644 --- a/sympy/polys/matrices/__init__.py +++ b/sympy/polys/matrices/__init__.py @@ -4,8 +4,8 @@ The main export from this package is the DomainMatrix class which is a lower-level implementation of matrices based on the polys Domains. This -implentation is typically a lot faster than sympy's standard Matrix class but -is a work in progress and is still experimental. +implementation is typically a lot faster than sympy's standard Matrix class +but is a work in progress and is still experimental. """ from .domainmatrix import DomainMatrix