normalforms.py 4.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156
  1. '''Functions returning normal forms of matrices'''
  2. from sympy.polys.domains.integerring import ZZ
  3. from sympy.polys.polytools import Poly
  4. from sympy.polys.matrices import DomainMatrix
  5. from sympy.polys.matrices.normalforms import (
  6. smith_normal_form as _snf,
  7. is_smith_normal_form as _is_snf,
  8. smith_normal_decomp as _snd,
  9. invariant_factors as _invf,
  10. hermite_normal_form as _hnf,
  11. )
  12. def _to_domain(m, domain=None):
  13. """Convert Matrix to DomainMatrix"""
  14. # XXX: deprecated support for RawMatrix:
  15. ring = getattr(m, "ring", None)
  16. m = m.applyfunc(lambda e: e.as_expr() if isinstance(e, Poly) else e)
  17. dM = DomainMatrix.from_Matrix(m)
  18. domain = domain or ring
  19. if domain is not None:
  20. dM = dM.convert_to(domain)
  21. return dM
  22. def smith_normal_form(m, domain=None):
  23. '''
  24. Return the Smith Normal Form of a matrix `m` over the ring `domain`.
  25. This will only work if the ring is a principal ideal domain.
  26. Examples
  27. ========
  28. >>> from sympy import Matrix, ZZ
  29. >>> from sympy.matrices.normalforms import smith_normal_form
  30. >>> m = Matrix([[12, 6, 4], [3, 9, 6], [2, 16, 14]])
  31. >>> print(smith_normal_form(m, domain=ZZ))
  32. Matrix([[1, 0, 0], [0, 10, 0], [0, 0, 30]])
  33. '''
  34. dM = _to_domain(m, domain)
  35. return _snf(dM).to_Matrix()
  36. def is_smith_normal_form(m, domain=None):
  37. '''
  38. Checks that the matrix is in Smith Normal Form
  39. '''
  40. dM = _to_domain(m, domain)
  41. return _is_snf(dM)
  42. def smith_normal_decomp(m, domain=None):
  43. '''
  44. Return the Smith Normal Decomposition of a matrix `m` over the ring
  45. `domain`. This will only work if the ring is a principal ideal domain.
  46. Examples
  47. ========
  48. >>> from sympy import Matrix, ZZ
  49. >>> from sympy.matrices.normalforms import smith_normal_decomp
  50. >>> m = Matrix([[12, 6, 4], [3, 9, 6], [2, 16, 14]])
  51. >>> a, s, t = smith_normal_decomp(m, domain=ZZ)
  52. >>> assert a == s * m * t
  53. '''
  54. dM = _to_domain(m, domain)
  55. a, s, t = _snd(dM)
  56. return a.to_Matrix(), s.to_Matrix(), t.to_Matrix()
  57. def invariant_factors(m, domain=None):
  58. '''
  59. Return the tuple of abelian invariants for a matrix `m`
  60. (as in the Smith-Normal form)
  61. References
  62. ==========
  63. .. [1] https://en.wikipedia.org/wiki/Smith_normal_form#Algorithm
  64. .. [2] https://web.archive.org/web/20200331143852/https://sierra.nmsu.edu/morandi/notes/SmithNormalForm.pdf
  65. '''
  66. dM = _to_domain(m, domain)
  67. factors = _invf(dM)
  68. factors = tuple(dM.domain.to_sympy(f) for f in factors)
  69. # XXX: deprecated.
  70. if hasattr(m, "ring"):
  71. if m.ring.is_PolynomialRing:
  72. K = m.ring
  73. to_poly = lambda f: Poly(f, K.symbols, domain=K.domain)
  74. factors = tuple(to_poly(f) for f in factors)
  75. return factors
  76. def hermite_normal_form(A, *, D=None, check_rank=False):
  77. r"""
  78. Compute the Hermite Normal Form of a Matrix *A* of integers.
  79. Examples
  80. ========
  81. >>> from sympy import Matrix
  82. >>> from sympy.matrices.normalforms import hermite_normal_form
  83. >>> m = Matrix([[12, 6, 4], [3, 9, 6], [2, 16, 14]])
  84. >>> print(hermite_normal_form(m))
  85. Matrix([[10, 0, 2], [0, 15, 3], [0, 0, 2]])
  86. Parameters
  87. ==========
  88. A : $m \times n$ ``Matrix`` of integers.
  89. D : int, optional
  90. Let $W$ be the HNF of *A*. If known in advance, a positive integer *D*
  91. being any multiple of $\det(W)$ may be provided. In this case, if *A*
  92. also has rank $m$, then we may use an alternative algorithm that works
  93. mod *D* in order to prevent coefficient explosion.
  94. check_rank : boolean, optional (default=False)
  95. The basic assumption is that, if you pass a value for *D*, then
  96. you already believe that *A* has rank $m$, so we do not waste time
  97. checking it for you. If you do want this to be checked (and the
  98. ordinary, non-modulo *D* algorithm to be used if the check fails), then
  99. set *check_rank* to ``True``.
  100. Returns
  101. =======
  102. ``Matrix``
  103. The HNF of matrix *A*.
  104. Raises
  105. ======
  106. DMDomainError
  107. If the domain of the matrix is not :ref:`ZZ`.
  108. DMShapeError
  109. If the mod *D* algorithm is used but the matrix has more rows than
  110. columns.
  111. References
  112. ==========
  113. .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
  114. (See Algorithms 2.4.5 and 2.4.8.)
  115. """
  116. # Accept any of Python int, SymPy Integer, and ZZ itself:
  117. if D is not None and not ZZ.of_type(D):
  118. D = ZZ(int(D))
  119. return _hnf(A._rep, D=D, check_rank=check_rank).to_Matrix()