matrices and groups: Smith Normal Form and an infinity test for fp groups

sympy/combinatorics/fp_groups.py

@@ -135,15 +135,39 @@ def order(self, strategy="relator_based"):
         2
 
         """
+        from sympy.core import S
         if self._order != None:
             return self._order
         if self._coset_table != None:
             self._order = len(self._coset_table.table)
         else:
-            self._coset_table = self.coset_enumeration([], strategy)
-            self._order = len(self._coset_table.table)
+            if self._is_infinite():
+                self._order = S.Infinity
+            else:
+                self._coset_table = self.coset_enumeration([], strategy)
+                self._order = len(self._coset_table.table)
         return self._order
 
+    def _is_infinite(self):
+        '''
+        Test if the group is infinite. Return `True` if the test succeeds
+        and `None` otherwise
+
+        '''
+        # Abelianisation test: check is the abelianisation is infinite
+        abelian_rels = []
+        from sympy.polys.solvers import RawMatrix as Matrix
+        from sympy.polys.domains import ZZ
+        from sympy.matrices.normalforms import smith_normal_invariants
+        for rel in self.relators():
+            abelian_rels.append([rel.exponent_sum(g) for g in self.generators])
+        m = Matrix(abelian_rels)
+        setattr(m, "ring", ZZ)
+        if 0 in smith_normal_invariants(m):
+            return True
+        else:
+            return None
+
     def index(self, H, strategy="relator_based"):
         """
         Returns the index of subgroup ``H`` in group ``self``.

sympy/combinatorics/free_groups.py

@@ -585,7 +585,8 @@ def order(self):
             return S.Infinity
 
     def commutator(self, other):
-        """Returns the commutator of `self` and `x`: ``~x*~self*x*self``
+        """
+        Return the commutator of `self` and `x`: ``~x*~self*x*self``
         """
         group = self.group
         if not isinstance(other, group.dtype):

sympy/combinatorics/tests/test_fp_groups.py

@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+from sympy import S
 from sympy.combinatorics.fp_groups import (FpGroup, CosetTable, low_index_subgroups,
         coset_enumeration_r, coset_enumeration_c, reidemeister_presentation)
 from sympy.combinatorics.free_groups import free_group
@@ -848,3 +849,8 @@ def test_subgroup_presentations():
         "a_0**5*b_1**-2*a_0*b_1**2*c_3**-1*b_1**-2*c_3**-1*b_1*a_0**5*b_1**-2*a_0*b_1**2*c_3**-1*b_1**-2*c_3**-1*b_1*a_0**5*b_1**-2*a_0*b_1**2*c_3**-1*b_1**-2*c_3**-1*b_1))"
         )
     assert str(reidemeister_presentation(f, H)) == k
+
+def test_is_infinite():
+    F, x, y = free_group("x, y")
+    f = FpGroup(F, [x*y*x**-1*y**-1, y**2])
+    assert f.order() == S.Infinity

sympy/matrices/normalforms.py

@@ -0,0 +1,145 @@
+from __future__ import print_function, division
+from sympy.matrices import diag
+
+'''
+Functions returning normal forms of matrices
+
+'''
+def smith_normal_form(m, domain = None):
+    '''
+    Return the Smith Normal Form of a matrix `m` over the ring `domain`.
+    This will only work if the ring is a principal ideal domain.
+
+    Examples
+    ========
+
+    >>> from sympy.polys.solvers import RawMatrix as Matrix
+    >>> from sympy.polys.domains import ZZ
+    >>> from sympy.matrices.normalforms import smith_normal_form
+    >>> m = Matrix([[12, 6, 4], [3, 9, 6], [2, 16, 14]])
+    >>> setattr(m, "ring", ZZ)
+    >>> print(smith_normal_form(m))
+    Matrix([[1, 0, 0], [0, 10, 0], [0, 0, -30]])
+
+    '''
+    invs = smith_normal_invariants(m, domain=domain)
+    smf = diag(*invs)
+    n = len(invs)
+    if m.rows > n:
+        smf = smf.row_insert(m.rows, zeros(m.rows-n, m.cols))
+    elif m.cols > n:
+        smf = smf.col_insert(m.cols, zeros(m.rows, m.cols-n))
+    return smf
+
+def smith_normal_invariants(m, domain = None):
+    '''
+    Return the tuple of abelian invariants for a matrix `m`
+    (as in the Smith-Normal form)
+
+    References
+    ==========
+
+    [1] https://en.wikipedia.org/wiki/Smith_normal_form#Algorithm
+    [2] http://sierra.nmsu.edu/morandi/notes/SmithNormalForm.pdf
+
+    '''
+    if not domain:
+        if not (hasattr(m, "ring") and m.ring.is_PID):
+            raise ValueError(
+                "The matrix entries must be over a principal ideal domain")
+        else:
+            domain = m.ring
+
+    if len(m) == 0:
+        return ()
+
+    m = m[:, :]
+
+    def add_rows(m, i, j, a, b, c, d):
+        # replace m[i, :] by a*m[i, :] + b*m[j, :]
+        # and m[j, :] by c*m[i, :] + d*m[j, :]
+        for k in range(m.cols):
+            e = m[i, k]
+            m[i, k] = a*e + b*m[j, k]
+            m[j, k] = c*e + d*m[j, k]
+
+    def add_columns(m, i, j, a, b, c, d):
+        # replace m[:, i] by a*m[:, i] + b*m[:, j]
+        # and m[:, j] by c*m[:, i] + d*m[:, j]
+        for k in range(m.rows):
+            e = m[k, i]
+            m[k, i] = a*e + b*m[k, j]
+            m[k, j] = c*e + d*m[k, j]
+
+    def clear_column(m):
+        # make m[1:, 0] zero by row and column operations
+        if m[0,0] == 0:
+            return m
+        pivot = m[0, 0]
+        for j in range(1, m.rows):
+            if m[j, 0] == 0:
+                continue
+            d, r = domain.div(m[j,0], pivot)
+            if r == 0:
+                add_rows(m, 0, j, 1, 0, -d, 1)
+            else:
+                a, b, g = domain.gcdex(pivot, m[j,0])
+                d_0 = domain.div(m[j, 0], g)[0]
+                d_j = domain.div(pivot, g)[0]
+                add_rows(m, 0, j, a, b, d_0, -d_j)
+                pivot = g
+        return m
+
+    def clear_row(m):
+        # make m[0, 1:] zero by row and column operations
+        if m[0] == 0:
+            return m
+        pivot = m[0, 0]
+        for j in range(1, m.cols):
+            if m[0, j] == 0:
+                continue
+            d, r = domain.div(m[0, j], pivot)
+            if r == 0:
+                add_columns(m, 0, j, 1, 0, -d, 1)
+            else:
+                a, b, g = domain.gcdex(pivot, m[0, j])
+                d_0 = domain.div(m[0, j], g)[0]
+                d_j = domain.div(pivot, g)[0]
+                add_columns(m, 0, j, a, b, d_0, -d_j)
+                pivot = g
+        return m
+
+    # permute the rows and columns until m[0,0] is non-zero if possible
+    ind = [i for i in range(m.rows) if m[i,0] != 0]
+    if ind:
+        m = m.permute_rows([[0, ind[0]]])
+    else:
+        ind = [j for j in range(m.cols) if m[0,j] != 0]
+        if ind:
+            m = m.permute_cols([[0, ind[0]]])
+
+    # make the first row and column except m[0,0] zero
+    while (any([m[0,i] != 0 for i in range(1,m.cols)]) or
+           any([m[i,0] != 0 for i in range(1,m.rows)])):
+        m = clear_column(m)
+        m = clear_row(m)
+
+    if 1 in m.shape:
+        invs = ()
+    else:
+        invs = smith_normal_invariants(m[1:,1:], domain=domain)
+
+    if m[0,0]:
+        result = [m[0,0]]
+        result.extend(invs)
+        # in case m[0] doesn't divide the invariants of the rest of the matrix
+        for i in range(len(result)-1):
+            if result[i] and domain.div(result[i+1], result[i])[1] != 0:
+                g = domain.gcd(result[i+1], result[i])
+                result[i+1] = domain.div(result[i], g)[0]*result[i+1]
+                result[i] = g
+            else:
+                break
+    else:
+        result = invs + (m[0,0],)
+    return tuple(result)

sympy/matrices/tests/test_normalforms.py

@@ -0,0 +1,20 @@
+from __future__ import print_function, division
+
+from sympy import Symbol, Poly
+from sympy.polys.solvers import RawMatrix as Matrix
+from sympy.matrices.normalforms import smith_normal_invariants, smith_normal_form
+from sympy.polys.domains import ZZ, QQ
+
+def test_smith_normal():
+    m = Matrix([[12, 6, 4,8],[3,9,6,12],[2,16,14,28],[20,10,10,20]])
+    setattr(m, 'ring', ZZ)
+    smf = Matrix([[1, 0, 0, 0], [0, 10, 0, 0], [0, 0, -30, 0], [0, 0, 0, 0]])
+    assert smith_normal_form(m) == smf
+
+    x = Symbol('x')
+    m = Matrix([[Poly(x-1), Poly(1, x),Poly(-1,x)],
+                [0, Poly(x), Poly(-1,x)],
+                [Poly(0,x),Poly(-1,x),Poly(x)]])
+    setattr(m, 'ring', QQ[x])
+    invs = (Poly(1, x), Poly(x - 1), Poly(x**2 - 1))
+    assert smith_normal_invariants(m) == invs

sympy/polys/domains/domain.py

@@ -43,6 +43,7 @@ class Domain(object):
 
     is_Simple = False
     is_Composite = False
+    is_PID = False
 
     has_CharacteristicZero = False
 

sympy/polys/domains/field.py

@@ -11,6 +11,7 @@ class Field(Ring):
     """Represents a field domain. """
 
     is_Field = True
+    is_PID = True
 
     def get_ring(self):
         """Returns a ring associated with ``self``. """

sympy/polys/domains/integerring.py

@@ -18,6 +18,7 @@ class IntegerRing(Ring, CharacteristicZero, SimpleDomain):
 
     is_IntegerRing = is_ZZ = True
     is_Numerical = True
+    is_PID = True
 
     has_assoc_Ring = True
     has_assoc_Field = True

sympy/polys/domains/polynomialring.py

@@ -33,6 +33,11 @@ def __init__(self, domain_or_ring, symbols=None, order=None):
         self.symbols = ring.symbols
         self.domain = ring.domain
 
+
+        if symbols:
+            if ring.domain.is_Field and ring.domain.is_Exact and len(symbols)==1:
+                self.is_PID = True
+
         # TODO: remove this
         self.dom = self.domain
 

sympy/polys/domains/realfield.py

@@ -21,6 +21,7 @@ class RealField(Field, CharacteristicZero, SimpleDomain):
 
     is_Exact = False
     is_Numerical = True
+    is_PID = False
 
     has_assoc_Ring = False
     has_assoc_Field = True

sympy/polys/solvers.py

@@ -7,6 +7,10 @@
 class RawMatrix(Matrix):
     _sympify = staticmethod(lambda x: x)
 
+    def is_zero():
+        from sympy.matrices import MatrixShaping
+        return MatrixShaping.is_zero(self)
+
 def eqs_to_matrix(eqs, ring):
     """Transform from equations to matrix form. """
     xs = ring.gens