Reidemeister Schreier algorithm implementation in SymPy

sympy/combinatorics/fp_groups.py

@@ -2,10 +2,12 @@
 from __future__ import print_function, division
 
 from sympy.core.basic import Basic
+from sympy.core import Symbol
 from sympy.printing.defaults import DefaultPrinting
 from sympy.utilities import public
 from sympy.utilities.iterables import flatten
-from sympy.combinatorics.free_group import FreeGroupElement
+from sympy.combinatorics.free_group import FreeGroupElement, free_group
+
 from itertools import chain, product
 from bisect import bisect_left
 
@@ -56,7 +58,7 @@ def __new__(cls, fr_grp, relators):
         obj._free_group = fr_grp
         obj._relators = relators
         obj.generators = obj._generators()
-        obj.dytype = type("FpGroupElement", (FpGroupElement,), {"group": obj})
+        obj.dtype = type("FpGroupElement", (FpGroupElement,), {"group": obj})
         return obj
 
     @property
@@ -328,10 +330,10 @@ def scan_check(self, alpha, word):
     def merge(self, k, lamda, q):
         p = self.p
         phi = self.rep(k)
-        chi = self.rep(lamda)
-        if phi != chi:
-            mu = min(phi, chi)
-            v = max(phi, chi)
+        psi = self.rep(lamda)
+        if phi != psi:
+            mu = min(phi, psi)
+            v = max(phi, psi)
             p[v] = mu
             q.append(v)
 
@@ -1003,4 +1005,263 @@ def first_in_class(C, Y=[]):
     return True
 
 
+# Pg 175 [1]
+def define_schreier_generators(C):
+    y = []
+    gamma = 1
+    f = C.fp_group
+    X = f.generators
+    C.P = [[None]*len(C.A) for i in range(C.n)]
+    for alpha, x in product(C.omega, C.A):
+        beta = C.table[alpha][C.A_dict[x]]
+        if beta == gamma:
+            C.P[alpha][C.A_dict[x]] = "<identity>"
+            C.P[beta][C.A_dict_inv[x]] = "<identity>"
+            gamma += 1
+        elif x in X and C.P[alpha][C.A_dict[x]] is None:
+            y_alpha_x = '%s_%s' % (x, alpha)
+            y.append(y_alpha_x)
+            C.P[alpha][C.A_dict[x]] = y_alpha_x
+    grp_gens = list(free_group(', '.join(y)))
+    C.schreier_free_group = grp_gens.pop(0)
+    C.schreier_generators = grp_gens
+    # replace all elements of P by, free group elements
+    for i, j in product(range(len(C.P)), range(len(C.A))):
+        # if equals "<identity>", replace by identity element
+        if C.P[i][j] == "<identity>":
+            C.P[i][j] = C.schreier_free_group.identity
+        elif isinstance(C.P[i][j], str):
+            r = C.schreier_generators[y.index(C.P[i][j])]
+            C.P[i][j] = r
+            beta = C.table[i][j]
+            C.P[beta][j + 1] = r**-1
+
+
+def reidemeister_relators(C):
+    R = C.fp_group.relators()
+    rels = set([rewrite(C, coset, word) for word in R for coset in range(C.n)])
+    identity = C.schreier_free_group.identity
+    order_1_gens = set([i for i in rels if len(i) == 1])
+    # remove all the order 1 generators from relators
+    rels.difference_update(order_1_gens)
+    order_2_gens = set([i for i in rels if len(i) == 2])
+    # replace order 1 generators by identity element in reidemeister relators
+    rels = [w.eliminate_word(x, identity) for w in rels for x in order_1_gens]
+    C.schreier_generators = [i for i in C.schreier_generators if i not in order_1_gens]
+    # remove cyclic conjugate elements from relators
+    i = 0
+    while i < len(rels):
+        j = i + 1
+        while j < len(rels):
+            if rels[i].is_cyclic_conjugate(rels[j]):
+                del rels[j]
+            else:
+                j += 1
+        i += 1
+    C.reidemeister_relators = rels
+
+
+def rewrite(C, alpha, w):
+    """
+    Parameters
+    ----------
+
+    C: CosetTable
+    α: A live coset
+    w: A word in `A*`
+
+    Returns
+    -------
+
+    ρ(τ(α), w)
+
+    Examples
+    ========
+
+    >>> from sympy.combinatorics.fp_groups import FpGroup, CosetTable, define_schreier_generators, rewrite
+    >>> from sympy.combinatorics.free_group import free_group
+    >>> F, x, y = free_group("x ,y")
+    >>> f = FpGroup(F, [x**2, y**3, (x*y)**6])
+    >>> C = CosetTable(f, [])
+    >>> C.table = [[1, 1, 2, 3], [0, 0, 4, 5], [4, 4, 3, 0], [5, 5, 0, 2], [2, 2, 5, 1], [3, 3, 1, 4]]
+    >>> C.p = [0, 1, 2, 3, 4, 5]
+    >>> define_schreier_generators(C)
+    >>> rewrite(C, 0, (x*y)**6)
+    x_4*y_2*x_3*x_1*x_2*y_4*x_5
+
+    """
+    v = C.schreier_free_group.identity
+    for i in range(len(w)):
+        x_i = w[i]
+        v = v*C.P[alpha][C.A_dict[x_i]]
+        alpha = C.table[alpha][C.A_dict[x_i]]
+    return v
+
+
+# routines for modified Todd Coxeter procedure
+# Section 5.3.2 from [1]
+def modified_define(C, alpha, x):
+    """
+    see also
+    ========
+    CosetTable.define
+
+    """
+    C.define(alpha, x)
+    C.P[alpha][C.A_dict[x]] = "<identity>"
+    C.p_p[beta] = "<identity>"
+
+def modified_scan(C, alpha, word, y):
+    """"
+    See also
+    ========
+    CosetTable.scan
+
+    """
+    table = C.table
+    f = alpha
+    # f_p and b_p are words with coset numbers "f" and "b".
+    # f_p when we scan forward from beginning, while b_p when we scan backward
+    # from end.
+    f_p = "<identity>"
+    b_p = y
+    i = 0
+    r = len(word)
+    b = alpha
+    j = r - 1
+    while i <= j and table[f][C.A_dict[word[i]]] is not None:
+        f_p = f_p*C.P[f][C.A_dict[word[i]]]
+        i += 1
+    if i > j:
+        if f != b:
+            modified_coincidence(C, f, alpha, f_p**-1*y)
+        return
+    while j >= i and table[b][C.A_dict_inv[word[j]]] is not None:
+        b_p = b_p*C.P[b][C.A_dict_inv[word[i]]]
+        b = table[b][C.A_dict_inv[word[i]]]
+        j -= 1
+    if j < i:
+        modified_coincidence(C, f, b, f_p**-1*b_p)
+    elif j == i:
+        # deduction making
+        table[f][C.A_dict[word[i]]] = b; table[b][C.A_dict_inv[word[i]]] = f
+        C.P[f][C.A_dict[word[i]]] = f_p**-1*b_p
+        C.P[b][C.A_dict_inv[word[i]]] = b_p**-1*f_p
+    # otherwise scan is incomplete and yields no information
+
+def modified_rep(C, k):
+    p = C.p
+    p_p = C.p_p
+    lamda = k
+    rho = p[lamda]
+    # `s` is used to trace back the compression path
+    s = [None]*len(p)
+    while rho != lamda:
+        s[rho] = lamda
+        lamda = rho
+        rho = p[lamda]
+    rho = s[lamda]
+    while rho != k:
+        mu = rho
+        rho = s[mu]
+        p[rho] = lamda
+        p_p[rho] = p_p[rho]*p_p[mu]
+    return lamda
+
+def modified_merge(C, k, lamda, w, q):
+    p = C.p
+    p_p = C.p_p
+    phi = C.rep(k)
+    psi = C.rep(lamda)
+    # can't use `max/min` function as in `CosetTable.merge`
+    if phi > psi:
+        p[phi] = psi
+        p_p[phi] = p_p[k]**-1*w*p_p[lamda]
+        q.append(psi)
+    elif psi > phi:
+        p[psi] = phi
+        p_p[psi] = p_p[lamda]**-1*w**-1*p_p[k]
+        q.append(psi)
+
+def modified_coincidence(C, alpha, beta, w):
+    A_dict = C.A_dict
+    A_dict_inv = C.A_dict_inv
+    table = C.table
+    p_p = C.p_p
+    l = 0
+    q = []
+    modified_merge(C, alpha, beta, w, q)
+    while len(q) > 0:
+        gamma = q.pop(0)
+        for x in A_dict:
+            delta = table[gamma][A_dict[x]]
+            if delta is not None:
+                table[delta][A_dict_inv[x]] = None
+                mu = C.rep(gamma)
+                nu = C.rep(delta)
+                if table[mu][A_dict[x]] is not None:
+                    v = p_p[delta]**-1*C.P[gamma][A_dict[x]]**-1*p_p[gamma]*C.P[mu][A_dict[x]]
+                    modified_merge(C, nu, table[mu][A_dict[x]], v, q)
+                elif table[nu][A_dict_inv[x]] is not None:
+                    v = p_p[gamma]**-1*C.P[gamma][A_dict[x]]*p_p[delta]*C.P[mu][A_dict_inv[x]]
+                else:
+                    table[mu][A_dict[x]] = nu; table[nu][A_dict_inv[x]] = mu
+                    v = p_p[gamma]**-1*C.P[gamma][A_dict[x]]*p_p[delta]
+                    C.P[mu][A_dict[x]] = v
+                    C.P[nu][A_dict_inv[x]] = v**-1
+
+
+def elimination_technique_1(C):
+    rels = list(C.reidemeister_relators)
+    # the shorter relators are examined first so that generators selected for
+    # elimination will have shorter strings as equivalent
+    rels.sort(reverse=True)
+    gens = list(C.schreier_generators)
+    # TODO: redundant generator can also be present as inverse in relator
+    redundant_gens = {}
+    # examine each relator in relator list for any generator occuring exactly
+    # once
+    for i in range(len(rels) -1, -1, -1):
+        j = 0
+        while j < len(gens):
+            if rels[i].generator_exponent_sum(gens[j]) == 1:
+                gen_index = rels[i].index(gens[j])
+                bk = rels[i].subword(gen_index + 1, len(rels[i]))
+                fw = rels[i].subword(0, gen_index)
+                redundant_gens[gens[j]] = (bk*fw)**-1
+                del rels[i]; del gens[j]
+            j += 1
+    # eliminate the redundant generator from remaing relators
+    for i in range(len(rels)):
+        w = rels[i]
+        for gen in redundant_gens:
+            rels[i] = rels[i].eliminate_word(gen, redundant_gens[gen])
+    return rels
+
+
+def reidemeister_presentation(fp_grp, H):
+    """
+    fp_group: A finitely presented group, an instance of FpGroup
+    H: A subgroup whose presentation is to be found, given as a list
+    of words in generators of `fp_grp`
+
+    Examples
+    ========
+
+    >>> from sympy.combinatorics.free_group import free_group, FpGroup
+    >>> F, x, y = free_group("x, y")
+    >>> f = FpGroup(F, [x**3, y**5, (x*y)**2])
+    >>> H = [x*y, x**-1*y**-1*x*y*x]
+    >>> reidemeister_relators(f, H)
+    [x_4**3, y_3**-3, y_2**-1*y_3**-1*y_2**-1*y_3**-1]
+
+    """
+    C = coset_enumeration_r(fp_grp, H)
+    C.compress(); C.standardize()
+    define_schreier_generators(C)
+    reidemeister_relators(C)
+    rels = elimination_technique_1(C)
+    return rels
+
+
 FpGroupElement = FreeGroupElement