groups: fix FpSubgroup's "contains" method

sympy/combinatorics/fp_groups.py

@@ -297,50 +297,125 @@ class FpSubgroup(DefaultPrinting):
     group belongs to the subgroup
 
     '''
-    def __init__(self, G, gens):
+    def __init__(self, G, gens, normal=False):
         super(FpSubgroup,self).__init__()
         self.parent = G
         self.generators = list(set([g for g in gens if g != G.identity]))
         self._min_words = None #for use in __contains__
         self.C = None
+        self.normal = normal
 
     def __contains__(self, g):
-        if self._min_words is None:
-            gens = self.generators[:]
-            gens.extend([e**-1 for e in gens])
-            for w1 in gens:
-                for w2 in gens:
-                    if w2**-1 == w1:
-                        continue
-                    if w2[len(w2)-1]**-1 == w1[0] and w2*w1 not in gens:
-                        gens.append(w2*w1)
-                    if w2[0]**-1 == w1[len(w1)-1] and w1*w2 not in gens:
-                        gens.append(w1*w2)
-            self._min_words = gens
-
-        min_words = self._min_words
-        known = {} #to keep track of words
-
-        def _word_break(w):
-            if len(w) == 0:
-                return True
-            i = 0
-            while i < len(w):
-                i += 1
-                prefix = w.subword(0, i)
-                if prefix not in min_words:
-                    continue
-                rest = w.subword(i, len(w))
-                if rest not in known:
-                    known[rest] = _word_break(rest)
-                if known[rest]:
-                    return True
-            return False
 
-        if _word_break(g):
-            return True
-        elif isinstance(self.parent, FreeGroup):
-            return False
+        if isinstance(self.parent, FreeGroup):
+            if self._min_words is None:
+                # make _min_words - a list of subwords such that
+                # g is in the subgroup if and only if it can be
+                # partitioned into these subwords. Infinite families of
+                # subwords are presented by tuples, e.g. (r, w)
+                # stands fot the family of subwords r*w**n*r**-1
+
+                def _process(w):
+                    # this is to be used before adding new words
+                    # into _min_words; if the word w is not cyclically
+                    # reduced, it will generate an infinite family of
+                    # subwords so should be written as a tuple;
+                    # if it is, w**-1 should be added to the list
+                    # as well
+                    p, r = w.cyclic_reduction(removed=True)
+                    if not r.is_identity:
+                        return [(r, p)]
+                    else:
+                        return [w, w**-1]
+
+                # make the initial list
+                gens = []
+                for w in self.generators:
+                    if self.normal:
+                        w = w.cyclic_reduction()
+                    gens.extend(_process(w))
+
+                for w1 in gens:
+                    for w2 in gens:
+                        # if w1 and w2 are equal or are inverses, continue
+                        if w1 == w2 or (not isinstance(w1, tuple)
+                                                        and w1**-1 == w2):
+                            continue
+
+                        # if the start of one word is the inverse of the
+                        # end of the other, their multuple should be added
+                        # to _min_words because of cancellation
+                        if isinstance(w1, tuple):
+                            # start, end
+                            s1, s2 = w1[0][0], w1[0][0]**-1
+                        else:
+                            s1, s2 = w1[0], w1[len(w1)-1]
+
+                        if isinstance(w2, tuple):
+                            # start, end
+                            r1, r2 = w2[0][0], w2[0][0]**-1
+                        else:
+                            r1, r2 = w2[0], w2[len(w1)-1]
+
+                        # p1 and p2 are w1 and w2 or, in case when
+                        # w1 or w2 is an infinite family, a representative
+                        p1, p2 = w1, w2
+                        if isinstance(w1, tuple):
+                            p1 = w1[0]*w1[1]*w1[0]**-1
+                        if isinstance(w2, tuple):
+                            p2 = w2[0]*w2[1]*w2[0]**-1
+
+                        # add the product of the words to the list is necessary
+                        if r1**-1 == s2 and not (p1*p2).is_identity:
+                            new = _process(p1*p2)
+                            if not new in gens:
+                                gens.extend(new)
+
+                        if r2**-1 == s1 and not (p2*p1).is_identity:
+                            new = _process(p2*p1)
+                            if not new in gens:
+                                gens.extend(new)
+
+                self._min_words = gens
+
+            min_words = self._min_words
+
+            def _is_subword(w):
+                # check if w is a word in _min_words or one of
+                # the infinite families in it
+                w, r = w.cyclic_reduction(removed=True)
+                if r.is_identity or self.normal:
+                    return w in min_words
+                else:
+                    t = [s[1] for s in min_words if isinstance(s, tuple)
+                                                                and s[0] == r]
+                    return [s for s in t if w.power_of(s)] != []
+
+            # store the solution of words for which the result of
+            # _word_break (below) is known
+            known = {}
+
+            def _word_break(w):
+                # check if w can be written as a product of words
+                # in min_words
+                if len(w) == 0:
+                    return True
+                i = 0
+                while i < len(w):
+                    i += 1
+                    prefix = w.subword(0, i)
+                    if not _is_subword(prefix):
+                        continue
+                    rest = w.subword(i, len(w))
+                    if rest not in known:
+                        known[rest] = _word_break(rest)
+                    if known[rest]:
+                        return True
+                return False
+
+            if self.normal:
+                g = g.cyclic_reduction()
+            return _word_break(g)
         else:
             if self.C is None:
                 C = self.parent.coset_enumeration(self.generators)
@@ -369,6 +444,15 @@ def to_FpGroup(self):
             return free_group(', '.join(gen_syms))[0]
         return self.parent.subgroup(C=self.C)
 
+    def __str__(self):
+        if len(self.generators) > 30:
+            str_form = "<fp subgroup with %s generators>" % len(self.generators)
+        else:
+            str_form = "<fp subgroup on the generators %s>" % str(self.generators)
+        return str_form
+
+    __repr__ = __str__
+
 
 ###############################################################################
 #                           LOW INDEX SUBGROUPS                               #

sympy/combinatorics/free_groups.py

@@ -1207,6 +1207,89 @@ def identity_cyclic_reduction(self):
             word = group.dtype(rep)
         return word
 
+    def cyclic_reduction(self, removed=False):
+        """Return a cyclically reduced version of the word. Unlike
+        `identity_cyclic_reduction`, this will not cyclically permute
+        the reduced word - just remove the "unreduced" bits on either
+        side of it. Compare the examples with those of
+        `identity_cyclic_reduction`.
+
+        When `removed` is `True`, return a tuple `(word, r)` where
+        self `r` is such that before the reductin the word was either
+        `r*word*r**-1`.
+
+        Examples
+        ========
+
+        >>> from sympy.combinatorics.free_groups import free_group
+        >>> F, x, y = free_group("x, y")
+        >>> (x**2*y**2*x**-1).cyclic_reduction()
+        x*y**2
+        >>> (x**-3*y**-1*x**5).cyclic_reduction()
+        y**-1*x**2
+        >>> (x**-3*y**-1*x**5).cyclic_reduction(removed=True)
+        (y**-1*x**2, x**-3)
+
+        """
+        word = self.copy()
+        group = self.group
+        g = self.group.identity
+        while not word.is_cyclically_reduced():
+            exp1 = abs(word.exponent_syllable(0))
+            exp2 = abs(word.exponent_syllable(-1))
+            exp = min(exp1, exp2)
+            start = word[0]**abs(exp)
+            end = word[-1]**abs(exp)
+            word = start**-1*word*end**-1
+            g = g*start
+        if removed:
+            return word, g
+        return word
+
+    def power_of(self, other):
+        '''
+        Check if `self == other**n` for some integer n.
+
+        Examples
+        ========
+
+        >>> from sympy.combinatorics.free_groups import free_group
+        >>> F, x, y = free_group("x, y")
+        >>> ((x*y)**2).power_of(x*y)
+        True
+        >>> (x**-3*y**-2*x**3).power_of(x**-3*y*x**3)
+        True
+
+        '''
+        if self.is_identity:
+            return True
+
+        l = len(other)
+        if l == 1:
+            # self has to be a power of one generator
+            gens = self.contains_generators()
+            s = other in gens or other**-1 in gens
+            return len(gens) == 1 and s
+
+        # if self is not cyclically reduced and it is a power of other,
+        # other isn't cyclically reduced and the parts removed during
+        # their reduction must be equal
+        reduced, r1 = self.cyclic_reduction(removed=True)
+        if not r1.is_identity:
+            other, r2 = other.cyclic_reduction(removed=True)
+            if r1 == r2:
+                return reduced.power_of(other)
+            return False
+
+        if len(self) < l or len(self) % l:
+            return False
+
+        prefix = self.subword(0, l)
+        if prefix == other or prefix**-1 == other:
+            rest = self.subword(l, len(self))
+            return rest.power_of(other)
+        return False
+
 
 def letter_form_to_array_form(array_form, group):
     """

sympy/combinatorics/homomorphisms.py

@@ -84,14 +84,14 @@ def _compute_kernel(self):
             raise NotImplementedError(
                 "Kernel computation is not implemented for infinite groups")
         gens = []
-        K = FpSubgroup(G, gens)
+        K = FpSubgroup(G, gens, normal=True)
         i = self.image().order()
         while K.order()*i != G_order:
             r = G.random_element()
             k = r*self.invert(self(r))
             if not k in K:
                 gens.append(k)
-                K = FpSubgroup(G, gens)
+                K = FpSubgroup(G, gens, normal=True)
         return K
 
     def image(self):

sympy/combinatorics/tests/test_fp_groups.py

@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 from sympy import S
 from sympy.combinatorics.fp_groups import (FpGroup, low_index_subgroups,
-                                           reidemeister_presentation)
+                                   reidemeister_presentation, FpSubgroup)
 from sympy.combinatorics.free_groups import free_group
 
 """
@@ -158,3 +158,13 @@ def test_order():
     F, a, b, c = free_group("a, b, c")
     f = FpGroup(F, [a**250, b**2, c*b*c**-1*b, c**4, c**-1*a**-1*c*a, a**-1*b**-1*a*b])
     assert f.order() == 2000
+
+def test_fp_subgroup():
+    F, x, y = free_group("x, y")
+    f = FpGroup(F, [x**4, y**2, x*y*x**-1*y])
+    S = FpSubgroup(f, [x*y])
+    assert (x*y)**-3 in S
+
+    S = FpSubgroup(F, [x**-1*y*x])
+    assert x**-1*y**4*x in S
+    assert x**-1*y**4*x**2 not in S