Another implementation of Groebner bases for the exterior algebra.

src/sage/algebras/clifford_algebra.py

@@ -30,6 +30,7 @@
 from sage.modules.with_basis.morphism import ModuleMorphismByLinearity
 from sage.categories.poor_man_map import PoorManMap
 from sage.rings.integer_ring import ZZ
+from sage.rings.noncommutative_ideals import Ideal_nc
 from sage.modules.free_module import FreeModule, FreeModule_generic
 from sage.matrix.constructor import Matrix
 from sage.matrix.args import MatrixArgs
@@ -2238,6 +2239,16 @@ def f4_reduce(P, G):
 
         return G
 
+    def _ideal_class_(self, n=0):
+        """
+        Return the class that is used to implement ideals of ``self``.
+
+        TESTS::
+
+            sage: E.<x,y,z> = ExteriorAlgebra(QQ)
+            sage: E._ideal_class_()
+        """
+        return ExteriorAlgebraIdeal
 
     class Element(CliffordAlgebraElement):
         """
@@ -2315,6 +2326,41 @@ def _mul_(self, other):
 
             return self.__class__(P, d)
 
+        def reduce(self, I, left=True):
+            r"""
+            Reduce ``self`` with respect to the elements in ``I``.
+
+            INPUT:
+
+            - ``I`` -- a list of exterior algebra elements or an ideal
+            - ``side`` -- the side, ignored if ``I`` is an ideal
+            """
+            if isinstance(I, ExteriorAlgebraIdeal):
+                I = I.groebner_basis()
+                left = (I.side() == "left")
+
+            E = self.parent()
+            f = self
+            from sage.algebras.exterior_algebra_cython import leading_support
+            for g in I:
+                lm = leading_support(g)
+                reduction = True
+                while reduction:
+                    supp = f.support()
+                    reduction = False
+                    for s in supp:
+                        if lm <= s:
+                            reduction = True
+                            mon = E.monomial(s - lm)
+                            if left:
+                                gp = mon * g
+                                f = f - f[s] / gp[s] * gp
+                            else:
+                                gp = g * mon
+                                f = f - f[s] / gp[s] * gp
+                            break
+            return f
+
         def interior_product(self, x):
             r"""
             Return the interior product (also known as antiderivation) of
@@ -3154,3 +3200,44 @@ def chain_complex(self, R=None):
             basis = next_basis
 
         return ChainComplex(data, degree=1)
+
+class ExteriorAlgebraIdeal(Ideal_nc):
+    """
+    An ideal of the exterior algebra.
+    """
+    def reduce(self, f):
+        """
+        Reduce ``f`` modulo ``self``.
+        """
+        return f.reduce(self.groebner_basis())
+
+    @cached_method
+    def groebner_basis(self):
+        r"""
+        Return the reduced Gröbner basis of ``self``.
+
+        EXAMPLES:
+
+        We compute an example that was checked against Macaulay2::
+
+            sage: E.<a,b,c,d,e> = ExteriorAlgebra(QQ)
+            sage: rels = [c*d*e - b*d*e + b*c*e - b*c*d,
+            ....:         c*d*e - a*d*e + a*c*e - a*c*d,
+            ....:         b*d*e - a*d*e + a*b*e - a*b*d,
+            ....:         b*c*e - a*c*e + a*b*e - a*b*c,
+            ....:         b*c*d - a*c*d + a*b*d - a*b*c]
+            sage: I = E.ideal(rels)
+            sage: I.groebner_basis()
+            (-b*c*d + b*c*e - b*d*e + c*d*e,
+             -a*c*d + a*c*e - a*d*e + c*d*e,
+             -a*b*d + a*b*e - a*d*e + b*d*e,
+             -a*b*c + a*b*e - a*c*e + b*c*e)
+        """
+        from sage.algebras.exterior_algebra_cython import compute_groebner
+        side = 2
+        if self.side() == "left":
+            side = 0
+        elif self.side() == "right":
+            side = 1
+        return compute_groebner(self, side)
+

src/sage/algebras/exterior_algebra_cython.pxd

@@ -0,0 +1,20 @@
+"""
+Exterior algebras backend
+"""
+
+from sage.data_structures.bitset cimport FrozenBitset
+from sage.rings.integer cimport Integer
+
+cdef inline Integer bitset_to_int(FrozenBitset X)
+cdef inline FrozenBitset int_to_bitset(Integer n)
+cdef inline unsigned long degree(FrozenBitset X)
+cdef inline FrozenBitset leading_supp(f)
+cpdef tuple get_leading_supports(tuple I)
+
+# Grobner basis functions
+cdef inline build_monomial(supp, E)
+cdef inline partial_S_poly(f, g, E, int side)
+cdef inline set preprocessing(list P, list G, E, int side)
+cdef inline list reduction(list P, list G, E, int side)
+#cpdef tuple compute_groebner(tuple I, int side)
+

src/sage/algebras/exterior_algebra_cython.pyx

@@ -0,0 +1,245 @@
+"""
+Exterior algebras backend
+
+This contains the backend implementations in Cython for the exterior algebra.
+"""
+
+from sage.libs.gmp.mpz cimport mpz_sizeinbase, mpz_setbit, mpz_tstbit, mpz_cmp_si, mpz_sgn
+from sage.data_structures.bitset_base cimport bitset_t, bitset_init, bitset_first, bitset_next, bitset_set_to
+from sage.structure.parent cimport Parent
+
+cdef inline Integer bitset_to_int(FrozenBitset X):
+    """
+    Convert ``X`` to an :class:`Integer`.
+    """
+    cdef Integer ret = Integer(0)
+    cdef long elt = bitset_first(X._bitset)
+    while elt >= 0:
+        mpz_setbit(ret.value, elt)
+        elt = bitset_next(X._bitset, elt + 1)
+    return ret
+
+cdef inline FrozenBitset int_to_bitset(Integer n):
+    """
+    Convert a nonnegative integer ``n`` to a :class:`FrozenBitset`.
+    """
+    cdef unsigned long i
+
+    if mpz_sgn(n.value) == 0:
+        return FrozenBitset()
+
+    cdef FrozenBitset ret = <FrozenBitset> FrozenBitset()
+    cdef size_t s = mpz_sizeinbase(n.value, 2)
+    bitset_init(ret._bitset, s)
+    for i in range(s):
+        bitset_set_to(ret._bitset, i, mpz_tstbit(n.value, i))
+    return ret
+
+
+cdef inline unsigned long degree(FrozenBitset X):
+    """
+    Compute the degree of ``X``.
+    """
+    cdef unsigned long ret = 0
+    cdef long elt = bitset_first(X._bitset)
+    while elt >= 0:
+        ret += 1
+        elt = bitset_next(X._bitset, elt + 1)
+    return ret
+
+
+#TODO: Bring the exterior algebra elements as a cdef class and make f know its type!
+cdef inline FrozenBitset leading_supp(f):
+    """
+    Return the leading support of the exterior algebra element ``f``.
+    """
+    cdef dict mc = <dict> f._monomial_coefficients
+    return int_to_bitset(min(bitset_to_int(k) for k in mc))
+
+def leading_support(f):
+    return leading_supp(f)
+
+cpdef tuple get_leading_supports(tuple I):
+    """
+    Return the leading supports of the elements in ``I``.
+
+    Used for testing mostly
+
+    INPUT:
+
+    - ``I`` -- a tuple of elements of an exterior algebra
+    """
+    # We filter out any elements that are 0
+    return tuple(set([leading_supp(f) for f in I if f._monomial_coefficients]))
+
+
+cdef inline build_monomial(E, supp):
+    """
+    Helper function for the fastest way to build a monomial.
+    """
+    return E.element_class(E, {supp: (<Parent> E)._base.one()})
+
+cdef inline partial_S_poly(f, g, E, int side):
+    """
+    Compute one half of the `S`-polynomial for ``f`` and ``g``.
+
+    This computes:
+
+    .. MATH::
+
+        LCM(LM(f), LM(g)) / LT(f) \cdot f.
+    """
+    cdef FrozenBitset lmf = leading_supp(f)
+    cdef FrozenBitset lmg = leading_supp(g)
+    cdef FrozenBitset D = <FrozenBitset> lmg.difference(lmf)
+    if side == 0:
+        ret = build_monomial(E, D) * f
+        return (~ret[lmf._union(lmg)]) * ret
+    ret = f * build_monomial(E, D)
+    return ret * (~ret[lmf._union(lmg)])
+
+cdef inline set preprocessing(list P, list G, E, int side):
+    """
+    Perform the preprocessing step.
+    """
+    #print("Start preprocessing:", P)
+    cdef set L = set(partial_S_poly(f0, f1, E, side) for f0,f1 in P)
+    L.update(partial_S_poly(f1, f0, E, side) for f0,f1 in P)
+    if side == 2:
+        # in partial_S_poly, side == 2 gets treated like right (== 1)
+        L.update(partial_S_poly(f0, f1, E, 0) for f0,f1 in P)
+        L.update(partial_S_poly(f1, f0, E, 0) for f0,f1 in P)
+
+    cdef set done = set(leading_supp(f) for f in L)
+    cdef set monL = set()
+    for f in L:
+        monL.update(f.support())
+    monL.difference_update(done)
+
+    while monL:
+        m = int_to_bitset(max(bitset_to_int(k) for k in monL))
+        done.add(m)
+        monL.remove(m)
+        for g in G:
+            lm = leading_supp(g)
+            if lm <= m:
+                f = build_monomial(E, m.difference(lm)) * g
+                if f in L:
+                    break
+                monL.update(set(f.support()) - done)
+                L.add(f)
+                break
+    #print("preprocessing:", L)
+    return L
+
+cdef inline list reduction(list P, list G, E, int side):
+    """
+    Perform the reduction of ``P`` mod ``G`` in ``E``.
+    """
+    cdef set L = preprocessing(P, G, E, side)
+    cdef Py_ssize_t i
+    from sage.matrix.constructor import matrix
+    M = matrix({(i, bitset_to_int(<FrozenBitset> m)): c for i,f in enumerate(L) for m,c in f._monomial_coefficients.items()},
+               sparse=True)
+    M.echelonize()  # Do this in place
+    lead_supports = set(leading_supp(f) for f in L)
+    return [E.element_class(E, {int_to_bitset(Integer(j)): c for j,c in M[i].iteritems()})
+            for i,p in enumerate(M.pivots())
+            if int_to_bitset(Integer(p)) not in lead_supports]
+
+def compute_groebner(I, side):
+    """
+    Compute the reduced ``side`` Gröbner basis for the ideal ``I``.
+
+    INPUT:
+
+    - ``I`` -- the ideal
+    - ``side`` -- integer; the side of the ideal: ``0`` for left, ``1`` for
+      right, and ``2`` for two-sided
+    """
+    E = I.ring()
+    cdef FrozenBitset p0, p1
+    cdef unsigned long deg
+    cdef Py_ssize_t i, j, k
+
+    cdef list G = [f for f in I.gens() if f]  # Remove 0s TODO: We should make this unnecessary here
+    cdef Py_ssize_t n = len(G)
+    cdef dict P = {}
+    cdef list Gp
+
+    # for ideals generated by homogeneous (wrt Z_2-grading) polynomials, we can just consider it as a left ideal
+    # TODO: We can reduce the number of S-poly computations for Z_2-graded homogeneous
+    #   ideals by throwing out those such that LCM(LM(f), LM(g)) == LM(f) * LM(g).
+    if all(f.is_super_homogeneous() for f in G):
+        side = 0
+
+    for i in range(n):
+        f0 = G[i]
+        p0 = leading_supp(f0)
+        for j in range(i+1, n):
+            f1 = G[j]
+            p1 = leading_supp(f1)
+            deg = degree(<FrozenBitset> (p0._union(p1)))
+            if deg in P:
+                P[deg].append((f0, f1))
+            else:
+                P[deg] = [(f0, f1)]
+
+    while P:
+        #print("Cur G:", G)
+        Pp = P.pop(min(P))  # The selection: lowest lcm degree
+        Gp = reduction(Pp, G, E, side)
+        #print("Reduction yielded:", Gp)
+        G.extend(Gp)
+        for j in range(n, len(G)):
+            f1 = G[j]
+            p1 = leading_supp(f1)
+            for i in range(j):
+                f0 = G[i]
+                p0 = leading_supp(f0)
+                deg = degree(<FrozenBitset> (p0._union(p1)))
+                if deg in P:
+                    P[deg].append((f0, f1))
+                else:
+                    P[deg] = [(f0, f1)]
+        n = len(G)
+
+    #print(G)
+
+    # Now that we have a Gröbner basis, we make this into a reduced Gröbner basis
+    cdef set pairs = set((i, j) for i in range(n) for j in range(n) if i != j)
+    cdef list supp
+    cdef bint did_reduction
+    cdef FrozenBitset lm, s
+    while pairs:
+        i,j = pairs.pop()
+        # We perform the classical reduction algorithm here on each pair
+        # TODO: Make this faster by using the previous technique
+        f = G[i]
+        g = G[j]
+        lm = leading_supp(g)
+        did_reduction = True
+        while did_reduction:
+            supp = f.support()
+            did_reduction = False
+            for s in supp:
+                if lm <= s:
+                    did_reduction = True
+                    mon = E.monomial(s - lm)
+                    if side == 0:
+                        gp = mon * g
+                        f = f - f[s] / gp[s] * gp
+                    else:
+                        gp = g * mon
+                        f = f - f[s] / gp[s] * gp
+                    break
+        if G[i] != f:
+            G[i] = f
+            #print("reduction:", G)
+            if not f:
+                pairs.difference_update((k, i) for k in range(n))
+            else:
+                pairs.update((k, i) for k in range(n) if k != i)
+
+    return tuple([f for f in G if f])
+