1"""
2
3Module for the DomainMatrix class.
4
5A DomainMatrix represents a matrix with elements that are in a particular
6Domain. Each DomainMatrix internally wraps a DDM which is used for the
7lower-level operations. The idea is that the DomainMatrix class provides the
8convenience routines for converting between Expr and the poly domains as well
9as unifying matrices with different domains.
10
11"""
12from functools import reduce
13from sympy.core.sympify import _sympify
14
15from ..constructor import construct_domain
16
17from .exceptions import (NonSquareMatrixError, ShapeError, DDMShapeError,
18        DDMDomainError, DDMFormatError, DDMBadInputError)
19
20from .ddm import DDM
21
22from .sdm import SDM
23
24from .domainscalar import DomainScalar
25
26from sympy.polys.domains import ZZ, EXRAW
27
28class DomainMatrix:
29    r"""
30    Associate Matrix with :py:class:`~.Domain`
31
32    Explanation
33    ===========
34
35    DomainMatrix uses :py:class:`~.Domain` for its internal representation
36    which makes it more faster for many common operations
37    than current sympy Matrix class, but this advantage makes it not
38    entirely compatible with Matrix.
39    DomainMatrix could be found analogous to numpy arrays with "dtype".
40    In the DomainMatrix, each matrix has a domain such as :ref:`ZZ`
41    or  :ref:`QQ(a)`.
42
43
44    Examples
45    ========
46
47    Creating a DomainMatrix from the existing Matrix class:
48
49    >>> from sympy import Matrix
50    >>> from sympy.polys.matrices import DomainMatrix
51    >>> Matrix1 = Matrix([
52    ...    [1, 2],
53    ...    [3, 4]])
54    >>> A = DomainMatrix.from_Matrix(Matrix1)
55    >>> A
56    DomainMatrix({0: {0: 1, 1: 2}, 1: {0: 3, 1: 4}}, (2, 2), ZZ)
57
58    Driectly forming a DomainMatrix:
59
60    >>> from sympy import ZZ
61    >>> from sympy.polys.matrices import DomainMatrix
62    >>> A = DomainMatrix([
63    ...    [ZZ(1), ZZ(2)],
64    ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
65    >>> A
66    DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ)
67
68    See Also
69    ========
70
71    DDM
72    SDM
73    Domain
74    Poly
75
76    """
77
78    def __new__(cls, rows, shape, domain, *, fmt=None):
79        """
80        Creates a :py:class:`~.DomainMatrix`.
81
82        Parameters
83        ==========
84
85        rows : Represents elements of DomainMatrix as list of lists
86        shape : Represents dimension of DomainMatrix
87        domain : Represents :py:class:`~.Domain` of DomainMatrix
88
89        Raises
90        ======
91
92        TypeError
93            If any of rows, shape and domain are not provided
94
95        """
96        if isinstance(rows, (DDM, SDM)):
97            raise TypeError("Use from_rep to initialise from SDM/DDM")
98        elif isinstance(rows, list):
99            rep = DDM(rows, shape, domain)
100        elif isinstance(rows, dict):
101            rep = SDM(rows, shape, domain)
102        else:
103            msg = "Input should be list-of-lists or dict-of-dicts"
104            raise TypeError(msg)
105
106        if fmt is not None:
107            if fmt == 'sparse':
108                rep = rep.to_sdm()
109            elif fmt == 'dense':
110                rep = rep.to_ddm()
111            else:
112                raise ValueError("fmt should be 'sparse' or 'dense'")
113
114        return cls.from_rep(rep)
115
116    def __getnewargs__(self):
117        rep = self.rep
118        if isinstance(rep, DDM):
119            arg = list(rep)
120        elif isinstance(rep, SDM):
121            arg = dict(rep)
122        else:
123            raise RuntimeError # pragma: no cover
124        return arg, self.shape, self.domain
125
126    def __getitem__(self, key):
127        i, j = key
128        m, n = self.shape
129        if not (isinstance(i, slice) or isinstance(j, slice)):
130            return DomainScalar(self.rep.getitem(i, j), self.domain)
131
132        if not isinstance(i, slice):
133            if not -m <= i < m:
134                raise IndexError("Row index out of range")
135            i = i % m
136            i = slice(i, i+1)
137        if not isinstance(j, slice):
138            if not -n <= j < n:
139                raise IndexError("Column index out of range")
140            j = j % n
141            j = slice(j, j+1)
142
143        return self.from_rep(self.rep.extract_slice(i, j))
144
145    def getitem_sympy(self, i, j):
146        return self.domain.to_sympy(self.rep.getitem(i, j))
147
148    def extract(self, rowslist, colslist):
149        return self.from_rep(self.rep.extract(rowslist, colslist))
150
151    def __setitem__(self, key, value):
152        i, j = key
153        if not self.domain.of_type(value):
154            raise TypeError
155        if isinstance(i, int) and isinstance(j, int):
156            self.rep.setitem(i, j, value)
157        else:
158            raise NotImplementedError
159
160    @classmethod
161    def from_rep(cls, rep):
162        """Create a new DomainMatrix efficiently from DDM/SDM.
163
164        Examples
165        ========
166
167        Create a :py:class:`~.DomainMatrix` with an dense internal
168        representation as :py:class:`~.DDM`:
169
170        >>> from sympy.polys.domains import ZZ
171        >>> from sympy.polys.matrices import DomainMatrix
172        >>> from sympy.polys.matrices.ddm import DDM
173        >>> drep = DDM([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)
174        >>> dM = DomainMatrix.from_rep(drep)
175        >>> dM
176        DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ)
177
178        Create a :py:class:`~.DomainMatrix` with a sparse internal
179        representation as :py:class:`~.SDM`:
180
181        >>> from sympy.polys.matrices import DomainMatrix
182        >>> from sympy.polys.matrices.sdm import SDM
183        >>> from sympy import ZZ
184        >>> drep = SDM({0:{1:ZZ(1)},1:{0:ZZ(2)}}, (2, 2), ZZ)
185        >>> dM = DomainMatrix.from_rep(drep)
186        >>> dM
187        DomainMatrix({0: {1: 1}, 1: {0: 2}}, (2, 2), ZZ)
188
189        Parameters
190        ==========
191
192        rep: SDM or DDM
193            The internal sparse or dense representation of the matrix.
194
195        Returns
196        =======
197
198        DomainMatrix
199            A :py:class:`~.DomainMatrix` wrapping *rep*.
200
201        Notes
202        =====
203
204        This takes ownership of rep as its internal representation. If rep is
205        being mutated elsewhere then a copy should be provided to
206        ``from_rep``. Only minimal verification or checking is done on *rep*
207        as this is supposed to be an efficient internal routine.
208
209        """
210        if not isinstance(rep, (DDM, SDM)):
211            raise TypeError("rep should be of type DDM or SDM")
212        self = super().__new__(cls)
213        self.rep = rep
214        self.shape = rep.shape
215        self.domain = rep.domain
216        return self
217
218    @classmethod
219    def from_list_sympy(cls, nrows, ncols, rows, **kwargs):
220        r"""
221        Convert a list of lists of Expr into a DomainMatrix using construct_domain
222
223        Parameters
224        ==========
225
226        nrows: number of rows
227        ncols: number of columns
228        rows: list of lists
229
230        Returns
231        =======
232
233        DomainMatrix containing elements of rows
234
235        Examples
236        ========
237
238        >>> from sympy.polys.matrices import DomainMatrix
239        >>> from sympy.abc import x, y, z
240        >>> A = DomainMatrix.from_list_sympy(1, 3, [[x, y, z]])
241        >>> A
242        DomainMatrix([[x, y, z]], (1, 3), ZZ[x,y,z])
243
244        See Also
245        ========
246
247        sympy.polys.constructor.construct_domain, from_dict_sympy
248
249        """
250        assert len(rows) == nrows
251        assert all(len(row) == ncols for row in rows)
252
253        items_sympy = [_sympify(item) for row in rows for item in row]
254
255        domain, items_domain = cls.get_domain(items_sympy, **kwargs)
256
257        domain_rows = [[items_domain[ncols*r + c] for c in range(ncols)] for r in range(nrows)]
258
259        return DomainMatrix(domain_rows, (nrows, ncols), domain)
260
261    @classmethod
262    def from_dict_sympy(cls, nrows, ncols, elemsdict, **kwargs):
263        """
264
265        Parameters
266        ==========
267
268        nrows: number of rows
269        ncols: number of cols
270        elemsdict: dict of dicts containing non-zero elements of the DomainMatrix
271
272        Returns
273        =======
274
275        DomainMatrix containing elements of elemsdict
276
277        Examples
278        ========
279
280        >>> from sympy.polys.matrices import DomainMatrix
281        >>> from sympy.abc import x,y,z
282        >>> elemsdict = {0: {0:x}, 1:{1: y}, 2: {2: z}}
283        >>> A = DomainMatrix.from_dict_sympy(3, 3, elemsdict)
284        >>> A
285        DomainMatrix({0: {0: x}, 1: {1: y}, 2: {2: z}}, (3, 3), ZZ[x,y,z])
286
287        See Also
288        ========
289
290        from_list_sympy
291
292        """
293        if not all(0 <= r < nrows for r in elemsdict):
294            raise DDMBadInputError("Row out of range")
295        if not all(0 <= c < ncols for row in elemsdict.values() for c in row):
296            raise DDMBadInputError("Column out of range")
297
298        items_sympy = [_sympify(item) for row in elemsdict.values() for item in row.values()]
299        domain, items_domain = cls.get_domain(items_sympy, **kwargs)
300
301        idx = 0
302        items_dict = {}
303        for i, row in elemsdict.items():
304                items_dict[i] = {}
305                for j in row:
306                    items_dict[i][j] = items_domain[idx]
307                    idx += 1
308
309        return DomainMatrix(items_dict, (nrows, ncols), domain)
310
311    @classmethod
312    def from_Matrix(cls, M, fmt='sparse',**kwargs):
313        r"""
314        Convert Matrix to DomainMatrix
315
316        Parameters
317        ==========
318
319        M: Matrix
320
321        Returns
322        =======
323
324        Returns DomainMatrix with identical elements as M
325
326        Examples
327        ========
328
329        >>> from sympy import Matrix
330        >>> from sympy.polys.matrices import DomainMatrix
331        >>> M = Matrix([
332        ...    [1.0, 3.4],
333        ...    [2.4, 1]])
334        >>> A = DomainMatrix.from_Matrix(M)
335        >>> A
336        DomainMatrix({0: {0: 1.0, 1: 3.4}, 1: {0: 2.4, 1: 1.0}}, (2, 2), RR)
337
338        We can keep internal representation as ddm using fmt='dense'
339        >>> from sympy import Matrix, QQ
340        >>> from sympy.polys.matrices import DomainMatrix
341        >>> A = DomainMatrix.from_Matrix(Matrix([[QQ(1, 2), QQ(3, 4)], [QQ(0, 1), QQ(0, 1)]]), fmt='dense')
342        >>> A.rep
343        [[1/2, 3/4], [0, 0]]
344
345        See Also
346        ========
347
348        Matrix
349
350        """
351        if fmt == 'dense':
352            return cls.from_list_sympy(*M.shape, M.tolist(), **kwargs)
353
354        return cls.from_dict_sympy(*M.shape, M.todod(), **kwargs)
355
356    @classmethod
357    def get_domain(cls, items_sympy, **kwargs):
358        K, items_K = construct_domain(items_sympy, **kwargs)
359        return K, items_K
360
361    def copy(self):
362        return self.from_rep(self.rep.copy())
363
364    def convert_to(self, K):
365        r"""
366        Change the domain of DomainMatrix to desired domain or field
367
368        Parameters
369        ==========
370
371        K : Represents the desired domain or field
372
373        Returns
374        =======
375
376        DomainMatrix
377            DomainMatrix with the desired domain or field
378
379        Examples
380        ========
381
382        >>> from sympy import ZZ, ZZ_I
383        >>> from sympy.polys.matrices import DomainMatrix
384        >>> A = DomainMatrix([
385        ...    [ZZ(1), ZZ(2)],
386        ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
387
388        >>> A.convert_to(ZZ_I)
389        DomainMatrix([[1, 2], [3, 4]], (2, 2), ZZ_I)
390
391        """
392        return self.from_rep(self.rep.convert_to(K))
393
394    def to_sympy(self):
395        return self.convert_to(EXRAW)
396
397    def to_field(self):
398        r"""
399        Returns a DomainMatrix with the appropriate field
400
401        Returns
402        =======
403
404        DomainMatrix
405            DomainMatrix with the appropriate field
406
407        Examples
408        ========
409
410        >>> from sympy import ZZ
411        >>> from sympy.polys.matrices import DomainMatrix
412        >>> A = DomainMatrix([
413        ...    [ZZ(1), ZZ(2)],
414        ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
415
416        >>> A.to_field()
417        DomainMatrix([[1, 2], [3, 4]], (2, 2), QQ)
418
419        """
420        K = self.domain.get_field()
421        return self.convert_to(K)
422
423    def to_sparse(self):
424        """
425        Return a sparse DomainMatrix representation of *self*.
426
427        Examples
428        ========
429
430        >>> from sympy.polys.matrices import DomainMatrix
431        >>> from sympy import QQ
432        >>> A = DomainMatrix([[1, 0],[0, 2]], (2, 2), QQ)
433        >>> A.rep
434        [[1, 0], [0, 2]]
435        >>> B = A.to_sparse()
436        >>> B.rep
437        {0: {0: 1}, 1: {1: 2}}
438        """
439        if self.rep.fmt == 'sparse':
440            return self
441
442        return self.from_rep(SDM.from_ddm(self.rep))
443
444    def to_dense(self):
445        """
446        Return a dense DomainMatrix representation of *self*.
447
448        Examples
449        ========
450
451        >>> from sympy.polys.matrices import DomainMatrix
452        >>> from sympy import QQ
453        >>> A = DomainMatrix({0: {0: 1}, 1: {1: 2}}, (2, 2), QQ)
454        >>> A.rep
455        {0: {0: 1}, 1: {1: 2}}
456        >>> B = A.to_dense()
457        >>> B.rep
458        [[1, 0], [0, 2]]
459
460        """
461        if self.rep.fmt == 'dense':
462            return self
463
464        return self.from_rep(SDM.to_ddm(self.rep))
465
466    @classmethod
467    def _unify_domain(cls, *matrices):
468        """Convert matrices to a common domain"""
469        domains = {matrix.domain for matrix in matrices}
470        if len(domains) == 1:
471            return matrices
472        domain = reduce(lambda x, y: x.unify(y), domains)
473        return tuple(matrix.convert_to(domain) for matrix in matrices)
474
475    @classmethod
476    def _unify_fmt(cls, *matrices, fmt=None):
477        """Convert matrices to the same format.
478
479        If all matrices have the same format, then return unmodified.
480        Otherwise convert both to the preferred format given as *fmt* which
481        should be 'dense' or 'sparse'.
482        """
483        formats = {matrix.rep.fmt for matrix in matrices}
484        if len(formats) == 1:
485            return matrices
486        if fmt == 'sparse':
487            return tuple(matrix.to_sparse() for matrix in matrices)
488        elif fmt == 'dense':
489            return tuple(matrix.to_dense() for matrix in matrices)
490        else:
491            raise ValueError("fmt should be 'sparse' or 'dense'")
492
493    def unify(self, *others, fmt=None):
494        """
495        Unifies the domains and the format of self and other
496        matrices.
497
498        Parameters
499        ==========
500
501        others : DomainMatrix
502
503        fmt: string 'dense', 'sparse' or `None` (default)
504            The preferred format to convert to if self and other are not
505            already in the same format. If `None` or not specified then no
506            conversion if performed.
507
508        Returns
509        =======
510
511        Tuple[DomainMatrix]
512            Matrices with unified domain and format
513
514        Examples
515        ========
516
517        Unify the domain of DomainMatrix that have different domains:
518
519        >>> from sympy import ZZ, QQ
520        >>> from sympy.polys.matrices import DomainMatrix
521        >>> A = DomainMatrix([[ZZ(1), ZZ(2)]], (1, 2), ZZ)
522        >>> B = DomainMatrix([[QQ(1, 2), QQ(2)]], (1, 2), QQ)
523        >>> Aq, Bq = A.unify(B)
524        >>> Aq
525        DomainMatrix([[1, 2]], (1, 2), QQ)
526        >>> Bq
527        DomainMatrix([[1/2, 2]], (1, 2), QQ)
528
529        Unify the format (dense or sparse):
530
531        >>> A = DomainMatrix([[ZZ(1), ZZ(2)]], (1, 2), ZZ)
532        >>> B = DomainMatrix({0:{0: ZZ(1)}}, (2, 2), ZZ)
533        >>> B.rep
534        {0: {0: 1}}
535
536        >>> A2, B2 = A.unify(B, fmt='dense')
537        >>> B2.rep
538        [[1, 0], [0, 0]]
539
540        See Also
541        ========
542
543        convert_to, to_dense, to_sparse
544
545        """
546        matrices = (self,) + others
547        matrices = DomainMatrix._unify_domain(*matrices)
548        if fmt is not None:
549            matrices = DomainMatrix._unify_fmt(*matrices, fmt=fmt)
550        return matrices
551
552    def to_Matrix(self):
553        r"""
554        Convert DomainMatrix to Matrix
555
556        Returns
557        =======
558
559        Matrix
560            MutableDenseMatrix for the DomainMatrix
561
562        Examples
563        ========
564
565        >>> from sympy import ZZ
566        >>> from sympy.polys.matrices import DomainMatrix
567        >>> A = DomainMatrix([
568        ...    [ZZ(1), ZZ(2)],
569        ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
570
571        >>> A.to_Matrix()
572        Matrix([
573            [1, 2],
574            [3, 4]])
575
576        See Also
577        ========
578
579        from_Matrix
580
581        """
582        from sympy.matrices.dense import MutableDenseMatrix
583        elemlist = self.rep.to_list()
584        elements_sympy = [self.domain.to_sympy(e) for row in elemlist for e in row]
585        return MutableDenseMatrix(*self.shape, elements_sympy)
586
587    def to_list(self):
588        return self.rep.to_list()
589
590    def to_list_flat(self):
591        return self.rep.to_list_flat()
592
593    def to_dok(self):
594        return self.rep.to_dok()
595
596    def __repr__(self):
597        return 'DomainMatrix(%s, %r, %r)' % (str(self.rep), self.shape, self.domain)
598
599    def transpose(self):
600        """Matrix transpose of ``self``"""
601        return self.from_rep(self.rep.transpose())
602
603    def flat(self):
604        rows, cols = self.shape
605        return [self[i,j].element for i in range(rows) for j in range(cols)]
606
607    @property
608    def is_zero_matrix(self):
609        return all(self[i, j].element == self.domain.zero for i in range(self.shape[0]) for j in range(self.shape[1]))
610
611    def hstack(A, *B):
612        r"""Horizontally stack the given matrices.
613
614        Parameters
615        ==========
616
617        B: DomainMatrix
618            Matrices to stack horizontally.
619
620        Returns
621        =======
622
623        DomainMatrix
624            DomainMatrix by stacking horizontally.
625
626        Examples
627        ========
628
629        >>> from sympy import ZZ
630        >>> from sympy.polys.matrices import DomainMatrix
631
632        >>> A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)
633        >>> B = DomainMatrix([[ZZ(5), ZZ(6)], [ZZ(7), ZZ(8)]], (2, 2), ZZ)
634        >>> A.hstack(B)
635        DomainMatrix([[1, 2, 5, 6], [3, 4, 7, 8]], (2, 4), ZZ)
636
637        >>> C = DomainMatrix([[ZZ(9), ZZ(10)], [ZZ(11), ZZ(12)]], (2, 2), ZZ)
638        >>> A.hstack(B, C)
639        DomainMatrix([[1, 2, 5, 6, 9, 10], [3, 4, 7, 8, 11, 12]], (2, 6), ZZ)
640
641        See Also
642        ========
643
644        unify
645        """
646        A, *B = A.unify(*B, fmt='dense')
647        return DomainMatrix.from_rep(A.rep.hstack(*(Bk.rep for Bk in B)))
648
649    def vstack(A, *B):
650        r"""Vertically stack the given matrices.
651
652        Parameters
653        ==========
654
655        B: DomainMatrix
656            Matrices to stack vertically.
657
658        Returns
659        =======
660
661        DomainMatrix
662            DomainMatrix by stacking vertically.
663
664        Examples
665        ========
666
667        >>> from sympy import ZZ
668        >>> from sympy.polys.matrices import DomainMatrix
669
670        >>> A = DomainMatrix([[ZZ(1), ZZ(2)], [ZZ(3), ZZ(4)]], (2, 2), ZZ)
671        >>> B = DomainMatrix([[ZZ(5), ZZ(6)], [ZZ(7), ZZ(8)]], (2, 2), ZZ)
672        >>> A.vstack(B)
673        DomainMatrix([[1, 2], [3, 4], [5, 6], [7, 8]], (4, 2), ZZ)
674
675        >>> C = DomainMatrix([[ZZ(9), ZZ(10)], [ZZ(11), ZZ(12)]], (2, 2), ZZ)
676        >>> A.vstack(B, C)
677        DomainMatrix([[1, 2], [3, 4], [5, 6], [7, 8], [9, 10], [11, 12]], (6, 2), ZZ)
678
679        See Also
680        ========
681
682        unify
683        """
684        A, *B = A.unify(*B, fmt='dense')
685        return DomainMatrix.from_rep(A.rep.vstack(*(Bk.rep for Bk in B)))
686
687    def applyfunc(self, func, domain=None):
688        if domain is None:
689            domain = self.domain
690        return self.from_rep(self.rep.applyfunc(func, domain))
691
692    def __add__(A, B):
693        if not isinstance(B, DomainMatrix):
694            return NotImplemented
695        A, B = A.unify(B, fmt='dense')
696        return A.add(B)
697
698    def __sub__(A, B):
699        if not isinstance(B, DomainMatrix):
700            return NotImplemented
701        A, B = A.unify(B, fmt='dense')
702        return A.sub(B)
703
704    def __neg__(A):
705        return A.neg()
706
707    def __mul__(A, B):
708        """A * B"""
709        if isinstance(B, DomainMatrix):
710            A, B = A.unify(B, fmt='dense')
711            return A.matmul(B)
712        elif B in A.domain:
713            return A.scalarmul(B)
714        elif isinstance(B, DomainScalar):
715            A, B = A.unify(B)
716            return A.scalarmul(B.element)
717        else:
718            return NotImplemented
719
720    def __rmul__(A, B):
721        if B in A.domain:
722            return A.rscalarmul(B)
723        elif isinstance(B, DomainScalar):
724            A, B = A.unify(B)
725            return A.rscalarmul(B.element)
726        else:
727            return NotImplemented
728
729    def __pow__(A, n):
730        """A ** n"""
731        if not isinstance(n, int):
732            return NotImplemented
733        return A.pow(n)
734
735    def _check(a, op, b, ashape, bshape):
736        if a.domain != b.domain:
737            msg = "Domain mismatch: %s %s %s" % (a.domain, op, b.domain)
738            raise DDMDomainError(msg)
739        if ashape != bshape:
740            msg = "Shape mismatch: %s %s %s" % (a.shape, op, b.shape)
741            raise DDMShapeError(msg)
742        if a.rep.fmt != b.rep.fmt:
743            msg = "Format mismatch: %s %s %s" % (a.rep.fmt, op, b.rep.fmt)
744            raise DDMFormatError(msg)
745
746    def add(A, B):
747        r"""
748        Adds two DomainMatrix matrices of the same Domain
749
750        Parameters
751        ==========
752
753        A, B: DomainMatrix
754            matrices to add
755
756        Returns
757        =======
758
759        DomainMatrix
760            DomainMatrix after Addition
761
762        Raises
763        ======
764
765        ShapeError
766            If the dimensions of the two DomainMatrix are not equal
767
768        ValueError
769            If the domain of the two DomainMatrix are not same
770
771        Examples
772        ========
773
774        >>> from sympy import ZZ
775        >>> from sympy.polys.matrices import DomainMatrix
776        >>> A = DomainMatrix([
777        ...    [ZZ(1), ZZ(2)],
778        ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
779        >>> B = DomainMatrix([
780        ...    [ZZ(4), ZZ(3)],
781        ...    [ZZ(2), ZZ(1)]], (2, 2), ZZ)
782
783        >>> A.add(B)
784        DomainMatrix([[5, 5], [5, 5]], (2, 2), ZZ)
785
786        See Also
787        ========
788
789        sub, matmul
790
791        """
792        A._check('+', B, A.shape, B.shape)
793        return A.from_rep(A.rep.add(B.rep))
794
795
796    def sub(A, B):
797        r"""
798        Subtracts two DomainMatrix matrices of the same Domain
799
800        Parameters
801        ==========
802
803        A, B: DomainMatrix
804            matrices to substract
805
806        Returns
807        =======
808
809        DomainMatrix
810            DomainMatrix after Substraction
811
812        Raises
813        ======
814
815        ShapeError
816            If the dimensions of the two DomainMatrix are not equal
817
818        ValueError
819            If the domain of the two DomainMatrix are not same
820
821        Examples
822        ========
823
824        >>> from sympy import ZZ
825        >>> from sympy.polys.matrices import DomainMatrix
826        >>> A = DomainMatrix([
827        ...    [ZZ(1), ZZ(2)],
828        ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
829        >>> B = DomainMatrix([
830        ...    [ZZ(4), ZZ(3)],
831        ...    [ZZ(2), ZZ(1)]], (2, 2), ZZ)
832
833        >>> A.sub(B)
834        DomainMatrix([[-3, -1], [1, 3]], (2, 2), ZZ)
835
836        See Also
837        ========
838
839        add, matmul
840
841        """
842        A._check('-', B, A.shape, B.shape)
843        return A.from_rep(A.rep.sub(B.rep))
844
845    def neg(A):
846        r"""
847        Returns the negative of DomainMatrix
848
849        Parameters
850        ==========
851
852        A : Represents a DomainMatrix
853
854        Returns
855        =======
856
857        DomainMatrix
858            DomainMatrix after Negation
859
860        Examples
861        ========
862
863        >>> from sympy import ZZ
864        >>> from sympy.polys.matrices import DomainMatrix
865        >>> A = DomainMatrix([
866        ...    [ZZ(1), ZZ(2)],
867        ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
868
869        >>> A.neg()
870        DomainMatrix([[-1, -2], [-3, -4]], (2, 2), ZZ)
871
872        """
873        return A.from_rep(A.rep.neg())
874
875    def mul(A, b):
876        r"""
877        Performs term by term multiplication for the second DomainMatrix
878        w.r.t first DomainMatrix. Returns a DomainMatrix whose rows are
879        list of DomainMatrix matrices created after term by term multiplication.
880
881        Parameters
882        ==========
883
884        A, B: DomainMatrix
885            matrices to multiply term-wise
886
887        Returns
888        =======
889
890        DomainMatrix
891            DomainMatrix after term by term multiplication
892
893        Examples
894        ========
895
896        >>> from sympy import ZZ
897        >>> from sympy.polys.matrices import DomainMatrix
898        >>> A = DomainMatrix([
899        ...    [ZZ(1), ZZ(2)],
900        ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
901        >>> B = DomainMatrix([
902        ...    [ZZ(1), ZZ(1)],
903        ...    [ZZ(0), ZZ(1)]], (2, 2), ZZ)
904
905        >>> A.mul(B)
906        DomainMatrix([[DomainMatrix([[1, 1], [0, 1]], (2, 2), ZZ),
907        DomainMatrix([[2, 2], [0, 2]], (2, 2), ZZ)],
908        [DomainMatrix([[3, 3], [0, 3]], (2, 2), ZZ),
909        DomainMatrix([[4, 4], [0, 4]], (2, 2), ZZ)]], (2, 2), ZZ)
910
911        See Also
912        ========
913
914        matmul
915
916        """
917        return A.from_rep(A.rep.mul(b))
918
919    def rmul(A, b):
920        return A.from_rep(A.rep.rmul(b))
921
922    def matmul(A, B):
923        r"""
924        Performs matrix multiplication of two DomainMatrix matrices
925
926        Parameters
927        ==========
928
929        A, B: DomainMatrix
930            to multiply
931
932        Returns
933        =======
934
935        DomainMatrix
936            DomainMatrix after multiplication
937
938        Examples
939        ========
940
941        >>> from sympy import ZZ
942        >>> from sympy.polys.matrices import DomainMatrix
943        >>> A = DomainMatrix([
944        ...    [ZZ(1), ZZ(2)],
945        ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
946        >>> B = DomainMatrix([
947        ...    [ZZ(1), ZZ(1)],
948        ...    [ZZ(0), ZZ(1)]], (2, 2), ZZ)
949
950        >>> A.matmul(B)
951        DomainMatrix([[1, 3], [3, 7]], (2, 2), ZZ)
952
953        See Also
954        ========
955
956        mul, pow, add, sub
957
958        """
959
960        A._check('*', B, A.shape[1], B.shape[0])
961        return A.from_rep(A.rep.matmul(B.rep))
962
963    def _scalarmul(A, lamda, reverse):
964        if lamda == A.domain.zero:
965            return DomainMatrix.zeros(A.shape, A.domain)
966        elif lamda == A.domain.one:
967            return A.copy()
968        elif reverse:
969            return A.rmul(lamda)
970        else:
971            return A.mul(lamda)
972
973    def scalarmul(A, lamda):
974        return A._scalarmul(lamda, reverse=False)
975
976    def rscalarmul(A, lamda):
977        return A._scalarmul(lamda, reverse=True)
978
979    def mul_elementwise(A, B):
980        assert A.domain == B.domain
981        return A.from_rep(A.rep.mul_elementwise(B.rep))
982
983    def __truediv__(A, lamda):
984        """ Method for Scalar Divison"""
985        if isinstance(lamda, int):
986            lamda = DomainScalar(ZZ(lamda), ZZ)
987
988        if not isinstance(lamda, DomainScalar):
989            return NotImplemented
990
991        A, lamda = A.to_field().unify(lamda)
992        if lamda.element == lamda.domain.zero:
993            raise ZeroDivisionError
994        if lamda.element == lamda.domain.one:
995            return A.to_field()
996
997        return A.mul(1 / lamda.element)
998
999    def pow(A, n):
1000        r"""
1001        Computes A**n
1002
1003        Parameters
1004        ==========
1005
1006        A : DomainMatrix
1007
1008        n : exponent for A
1009
1010        Returns
1011        =======
1012
1013        DomainMatrix
1014            DomainMatrix on computing A**n
1015
1016        Raises
1017        ======
1018
1019        NotImplementedError
1020            if n is negative.
1021
1022        Examples
1023        ========
1024
1025        >>> from sympy import ZZ
1026        >>> from sympy.polys.matrices import DomainMatrix
1027        >>> A = DomainMatrix([
1028        ...    [ZZ(1), ZZ(1)],
1029        ...    [ZZ(0), ZZ(1)]], (2, 2), ZZ)
1030
1031        >>> A.pow(2)
1032        DomainMatrix([[1, 2], [0, 1]], (2, 2), ZZ)
1033
1034        See Also
1035        ========
1036
1037        matmul
1038
1039        """
1040        nrows, ncols = A.shape
1041        if nrows != ncols:
1042            raise NonSquareMatrixError('Power of a nonsquare matrix')
1043        if n < 0:
1044            raise NotImplementedError('Negative powers')
1045        elif n == 0:
1046            return A.eye(nrows, A.domain)
1047        elif n == 1:
1048            return A
1049        elif n % 2 == 1:
1050            return A * A**(n - 1)
1051        else:
1052            sqrtAn = A ** (n // 2)
1053            return sqrtAn * sqrtAn
1054
1055    def scc(self):
1056        """Compute the strongly connected components of a DomainMatrix
1057
1058        Explanation
1059        ===========
1060
1061        A square matrix can be considered as the adjacency matrix for a
1062        directed graph where the row and column indices are the vertices. In
1063        this graph if there is an edge from vertex ``i`` to vertex ``j`` if
1064        ``M[i, j]`` is nonzero. This routine computes the strongly connected
1065        components of that graph which are subsets of the rows and columns that
1066        are connected by some nonzero element of the matrix. The strongly
1067        connected components are useful because many operations such as the
1068        determinant can be computed by working with the submatrices
1069        corresponding to each component.
1070
1071        Examples
1072        ========
1073
1074        Find the strongly connected components of a matrix:
1075
1076        >>> from sympy import ZZ
1077        >>> from sympy.polys.matrices import DomainMatrix
1078        >>> M = DomainMatrix([[ZZ(1), ZZ(0), ZZ(2)],
1079        ...                   [ZZ(0), ZZ(3), ZZ(0)],
1080        ...                   [ZZ(4), ZZ(6), ZZ(5)]], (3, 3), ZZ)
1081        >>> M.scc()
1082        [[1], [0, 2]]
1083
1084        Compute the determinant from the components:
1085
1086        >>> MM = M.to_Matrix()
1087        >>> MM
1088        Matrix([
1089        [1, 0, 2],
1090        [0, 3, 0],
1091        [4, 6, 5]])
1092        >>> MM[[1], [1]]
1093        Matrix([[3]])
1094        >>> MM[[0, 2], [0, 2]]
1095        Matrix([
1096        [1, 2],
1097        [4, 5]])
1098        >>> MM.det()
1099        -9
1100        >>> MM[[1], [1]].det() * MM[[0, 2], [0, 2]].det()
1101        -9
1102
1103        The components are given in reverse topological order and represent a
1104        permutation of the rows and columns that will bring the matrix into
1105        block lower-triangular form:
1106
1107        >>> MM[[1, 0, 2], [1, 0, 2]]
1108        Matrix([
1109        [3, 0, 0],
1110        [0, 1, 2],
1111        [6, 4, 5]])
1112
1113        Returns
1114        =======
1115
1116        List of lists of integers
1117            Each list represents a strongly connected component.
1118
1119        See also
1120        ========
1121
1122        sympy.matrices.matrices.MatrixBase.strongly_connected_components
1123        sympy.utilities.iterables.strongly_connected_components
1124
1125        """
1126        rows, cols = self.shape
1127        assert rows == cols
1128        return self.rep.scc()
1129
1130    def rref(self):
1131        r"""
1132        Returns reduced-row echelon form and list of pivots for the DomainMatrix
1133
1134        Returns
1135        =======
1136
1137        (DomainMatrix, list)
1138            reduced-row echelon form and list of pivots for the DomainMatrix
1139
1140        Raises
1141        ======
1142
1143        ValueError
1144            If the domain of DomainMatrix not a Field
1145
1146        Examples
1147        ========
1148
1149        >>> from sympy import QQ
1150        >>> from sympy.polys.matrices import DomainMatrix
1151        >>> A = DomainMatrix([
1152        ...     [QQ(2), QQ(-1), QQ(0)],
1153        ...     [QQ(-1), QQ(2), QQ(-1)],
1154        ...     [QQ(0), QQ(0), QQ(2)]], (3, 3), QQ)
1155
1156        >>> rref_matrix, rref_pivots = A.rref()
1157        >>> rref_matrix
1158        DomainMatrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]], (3, 3), QQ)
1159        >>> rref_pivots
1160        (0, 1, 2)
1161
1162        See Also
1163        ========
1164
1165        convert_to, lu
1166
1167        """
1168        if not self.domain.is_Field:
1169            raise ValueError('Not a field')
1170        rref_ddm, pivots = self.rep.rref()
1171        return self.from_rep(rref_ddm), tuple(pivots)
1172
1173    def nullspace(self):
1174        r"""
1175        Returns the Null Space for the DomainMatrix
1176
1177        Returns
1178        =======
1179
1180        DomainMatrix
1181            Null Space of the DomainMatrix
1182
1183        Examples
1184        ========
1185
1186        >>> from sympy import QQ
1187        >>> from sympy.polys.matrices import DomainMatrix
1188        >>> A = DomainMatrix([
1189        ...    [QQ(1), QQ(-1)],
1190        ...    [QQ(2), QQ(-2)]], (2, 2), QQ)
1191        >>> A.nullspace()
1192        DomainMatrix([[1, 1]], (1, 2), QQ)
1193
1194        """
1195        if not self.domain.is_Field:
1196            raise ValueError('Not a field')
1197        return self.from_rep(self.rep.nullspace()[0])
1198
1199    def inv(self):
1200        r"""
1201        Finds the inverse of the DomainMatrix if exists
1202
1203        Returns
1204        =======
1205
1206        DomainMatrix
1207            DomainMatrix after inverse
1208
1209        Raises
1210        ======
1211
1212        ValueError
1213            If the domain of DomainMatrix not a Field
1214
1215        NonSquareMatrixError
1216            If the DomainMatrix is not a not Square DomainMatrix
1217
1218        Examples
1219        ========
1220
1221        >>> from sympy import QQ
1222        >>> from sympy.polys.matrices import DomainMatrix
1223        >>> A = DomainMatrix([
1224        ...     [QQ(2), QQ(-1), QQ(0)],
1225        ...     [QQ(-1), QQ(2), QQ(-1)],
1226        ...     [QQ(0), QQ(0), QQ(2)]], (3, 3), QQ)
1227        >>> A.inv()
1228        DomainMatrix([[2/3, 1/3, 1/6], [1/3, 2/3, 1/3], [0, 0, 1/2]], (3, 3), QQ)
1229
1230        See Also
1231        ========
1232
1233        neg
1234
1235        """
1236        if not self.domain.is_Field:
1237            raise ValueError('Not a field')
1238        m, n = self.shape
1239        if m != n:
1240            raise NonSquareMatrixError
1241        inv = self.rep.inv()
1242        return self.from_rep(inv)
1243
1244    def det(self):
1245        r"""
1246        Returns the determinant of a Square DomainMatrix
1247
1248        Returns
1249        =======
1250
1251        S.Complexes
1252            determinant of Square DomainMatrix
1253
1254        Raises
1255        ======
1256
1257        ValueError
1258            If the domain of DomainMatrix not a Field
1259
1260        Examples
1261        ========
1262
1263        >>> from sympy import ZZ
1264        >>> from sympy.polys.matrices import DomainMatrix
1265        >>> A = DomainMatrix([
1266        ...    [ZZ(1), ZZ(2)],
1267        ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
1268
1269        >>> A.det()
1270        -2
1271
1272        """
1273        m, n = self.shape
1274        if m != n:
1275            raise NonSquareMatrixError
1276        return self.rep.det()
1277
1278    def lu(self):
1279        r"""
1280        Returns Lower and Upper decomposition of the DomainMatrix
1281
1282        Returns
1283        =======
1284
1285        (L, U, exchange)
1286            L, U are Lower and Upper decomposition of the DomainMatrix,
1287            exchange is the list of indices of rows exchanged in the decomposition.
1288
1289        Raises
1290        ======
1291
1292        ValueError
1293            If the domain of DomainMatrix not a Field
1294
1295        Examples
1296        ========
1297
1298        >>> from sympy import QQ
1299        >>> from sympy.polys.matrices import DomainMatrix
1300        >>> A = DomainMatrix([
1301        ...    [QQ(1), QQ(-1)],
1302        ...    [QQ(2), QQ(-2)]], (2, 2), QQ)
1303        >>> A.lu()
1304        (DomainMatrix([[1, 0], [2, 1]], (2, 2), QQ), DomainMatrix([[1, -1], [0, 0]], (2, 2), QQ), [])
1305
1306        See Also
1307        ========
1308
1309        lu_solve
1310
1311        """
1312        if not self.domain.is_Field:
1313            raise ValueError('Not a field')
1314        L, U, swaps = self.rep.lu()
1315        return self.from_rep(L), self.from_rep(U), swaps
1316
1317    def lu_solve(self, rhs):
1318        r"""
1319        Solver for DomainMatrix x in the A*x = B
1320
1321        Parameters
1322        ==========
1323
1324        rhs : DomainMatrix B
1325
1326        Returns
1327        =======
1328
1329        DomainMatrix
1330            x in A*x = B
1331
1332        Raises
1333        ======
1334
1335        ShapeError
1336            If the DomainMatrix A and rhs have different number of rows
1337
1338        ValueError
1339            If the domain of DomainMatrix A not a Field
1340
1341        Examples
1342        ========
1343
1344        >>> from sympy import QQ
1345        >>> from sympy.polys.matrices import DomainMatrix
1346        >>> A = DomainMatrix([
1347        ...    [QQ(1), QQ(2)],
1348        ...    [QQ(3), QQ(4)]], (2, 2), QQ)
1349        >>> B = DomainMatrix([
1350        ...    [QQ(1), QQ(1)],
1351        ...    [QQ(0), QQ(1)]], (2, 2), QQ)
1352
1353        >>> A.lu_solve(B)
1354        DomainMatrix([[-2, -1], [3/2, 1]], (2, 2), QQ)
1355
1356        See Also
1357        ========
1358
1359        lu
1360
1361        """
1362        if self.shape[0] != rhs.shape[0]:
1363            raise ShapeError("Shape")
1364        if not self.domain.is_Field:
1365            raise ValueError('Not a field')
1366        sol = self.rep.lu_solve(rhs.rep)
1367        return self.from_rep(sol)
1368
1369    def _solve(A, b):
1370        # XXX: Not sure about this method or its signature. It is just created
1371        # because it is needed by the holonomic module.
1372        if A.shape[0] != b.shape[0]:
1373            raise ShapeError("Shape")
1374        if A.domain != b.domain or not A.domain.is_Field:
1375            raise ValueError('Not a field')
1376        Aaug = A.hstack(b)
1377        Arref, pivots = Aaug.rref()
1378        particular = Arref.from_rep(Arref.rep.particular())
1379        nullspace_rep, nonpivots = Arref[:,:-1].rep.nullspace()
1380        nullspace = Arref.from_rep(nullspace_rep)
1381        return particular, nullspace
1382
1383    def charpoly(self):
1384        r"""
1385        Returns the coefficients of the characteristic polynomial
1386        of the DomainMatrix. These elements will be domain elements.
1387        The domain of the elements will be same as domain of the DomainMatrix.
1388
1389        Returns
1390        =======
1391
1392        list
1393            coefficients of the characteristic polynomial
1394
1395        Raises
1396        ======
1397
1398        NonSquareMatrixError
1399            If the DomainMatrix is not a not Square DomainMatrix
1400
1401        Examples
1402        ========
1403
1404        >>> from sympy import ZZ
1405        >>> from sympy.polys.matrices import DomainMatrix
1406        >>> A = DomainMatrix([
1407        ...    [ZZ(1), ZZ(2)],
1408        ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
1409
1410        >>> A.charpoly()
1411        [1, -5, -2]
1412
1413        """
1414        m, n = self.shape
1415        if m != n:
1416            raise NonSquareMatrixError("not square")
1417        return self.rep.charpoly()
1418
1419    @classmethod
1420    def eye(cls, shape, domain):
1421        r"""
1422        Return identity matrix of size n
1423
1424        Examples
1425        ========
1426
1427        >>> from sympy.polys.matrices import DomainMatrix
1428        >>> from sympy import QQ
1429        >>> DomainMatrix.eye(3, QQ)
1430        DomainMatrix({0: {0: 1}, 1: {1: 1}, 2: {2: 1}}, (3, 3), QQ)
1431
1432        """
1433        if isinstance(shape, int):
1434            shape = (shape, shape)
1435        return cls.from_rep(SDM.eye(shape, domain))
1436
1437    @classmethod
1438    def diag(cls, diagonal, domain, shape=None):
1439        r"""
1440        Return diagonal matrix with entries from ``diagonal``.
1441
1442        Examples
1443        ========
1444
1445        >>> from sympy.polys.matrices import DomainMatrix
1446        >>> from sympy import ZZ
1447        >>> DomainMatrix.diag([ZZ(5), ZZ(6)], ZZ)
1448        DomainMatrix({0: {0: 5}, 1: {1: 6}}, (2, 2), ZZ)
1449
1450        """
1451        if shape is None:
1452            N = len(diagonal)
1453            shape = (N, N)
1454        return cls.from_rep(SDM.diag(diagonal, domain, shape))
1455
1456    @classmethod
1457    def zeros(cls, shape, domain, *, fmt='sparse'):
1458        """Returns a zero DomainMatrix of size shape, belonging to the specified domain
1459
1460        Examples
1461        ========
1462
1463        >>> from sympy.polys.matrices import DomainMatrix
1464        >>> from sympy import QQ
1465        >>> DomainMatrix.zeros((2, 3), QQ)
1466        DomainMatrix({}, (2, 3), QQ)
1467
1468        """
1469        return cls.from_rep(SDM.zeros(shape, domain))
1470
1471    @classmethod
1472    def ones(cls, shape, domain):
1473        """Returns a zero DomainMatrix of size shape, belonging to the specified domain
1474
1475        Examples
1476        ========
1477
1478        >>> from sympy.polys.matrices import DomainMatrix
1479        >>> from sympy import QQ
1480        >>> DomainMatrix.ones((2,3), QQ)
1481        DomainMatrix([[1, 1, 1], [1, 1, 1]], (2, 3), QQ)
1482
1483        """
1484        return cls.from_rep(DDM.ones(shape, domain))
1485
1486    def __eq__(A, B):
1487        r"""
1488        Checks for two DomainMatrix matrices to be equal or not
1489
1490        Parameters
1491        ==========
1492
1493        A, B: DomainMatrix
1494            to check equality
1495
1496        Returns
1497        =======
1498
1499        Boolean
1500            True for equal, else False
1501
1502        Raises
1503        ======
1504
1505        NotImplementedError
1506            If B is not a DomainMatrix
1507
1508        Examples
1509        ========
1510
1511        >>> from sympy import ZZ
1512        >>> from sympy.polys.matrices import DomainMatrix
1513        >>> A = DomainMatrix([
1514        ...    [ZZ(1), ZZ(2)],
1515        ...    [ZZ(3), ZZ(4)]], (2, 2), ZZ)
1516        >>> B = DomainMatrix([
1517        ...    [ZZ(1), ZZ(1)],
1518        ...    [ZZ(0), ZZ(1)]], (2, 2), ZZ)
1519        >>> A.__eq__(A)
1520        True
1521        >>> A.__eq__(B)
1522        False
1523
1524        """
1525        if not isinstance(A, type(B)):
1526            return NotImplemented
1527        return A.domain == B.domain and A.rep == B.rep
1528
1529    def unify_eq(A, B):
1530        if A.shape != B.shape:
1531            return False
1532        if A.domain != B.domain:
1533            A, B = A.unify(B)
1534        return A == B
1535