From d0ab1a1613d9152c7a518a17b405e3b9b615e0f4 Mon Sep 17 00:00:00 2001 From: "Trevor K. Karn" Date: Tue, 21 Jun 2022 12:31:12 -0500 Subject: [PATCH 01/27] Add _basis_index_function as a method, fix some doctests, and fix some type issues (tuple -> FrozenBitset) --- src/sage/algebras/clifford_algebra.py | 61 +++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index dfdafc49a06..6aa6b541724 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -885,9 +885,15 @@ def _element_constructor_(self, x): if x in self.free_module(): R = self.base_ring() if x.parent().base_ring() is R: +<<<<<<< HEAD return self.element_class(self, {FrozenBitset((i, )): c for i, c in x.items()}) # if the base ring is different, attempt to coerce it into R return self.element_class(self, {FrozenBitset((i, )): R(c) for i, c in x.items() if R(c) != R.zero()}) +======= + return self.element_class(self, {FrozenBitset((i,)): c for i,c in x.items()}) + # if the base ring is different, attempt to coerce it into R + return self.element_class(self, {FrozenBitset((i,)): R(c) for i,c in x.items() if R(c) != R.zero()}) +>>>>>>> a7cd8bf99e (Add _basis_index_function as a method, fix some doctests, and fix some type issues (tuple -> FrozenBitset)) if (isinstance(x, CliffordAlgebraElement) and self.has_coerce_map_from(x.parent())): @@ -898,10 +904,55 @@ def _element_constructor_(self, x): R = self.base_ring() return self.element_class(self, {FrozenBitset((i,)): R.one() for i in x}) +<<<<<<< HEAD try: return super(CliffordAlgebra, self)._element_constructor_(x) except TypeError: raise TypeError(f'do not know how to make {x=} an element of self') +======= + if isinstance(x, tuple): + R = self.base_ring() + return self.element_class(self, {FrozenBitset((i,)): R.one() for i in x}) + + return super(CliffordAlgebra, self)._element_constructor_(x) +>>>>>>> a7cd8bf99e (Add _basis_index_function as a method, fix some doctests, and fix some type issues (tuple -> FrozenBitset)) + + def _basis_index_function(self, x): + """ + Given an integer indexing the basis, return the correct + bitset. + + For backwards compatibility, tuples are also accepted. + + EXAMPLES:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl = CliffordAlgebra(Q) + sage: Cl._basis_index_function(7) + 111 + sage: Cl._basis_index_function(5) + 101 + sage: Cl._basis_index_function(4) + 001 + + sage: Cl._basis_index_function((0, 1, 2)) + 111 + sage: Cl._basis_index_function((0, 2)) + 101 + sage: Cl._basis_index_function((2,)) + 001 + """ + Q = self._quadratic_form + format_style = f"0{Q.dim()}b" + + # if the input is a tuple, assume that it has + # entries in {0, ..., 2**Q.dim()-1} + if isinstance(x, tuple): + return FrozenBitset(x, capacity = Q.dim()) + + # slice the output of format in order to make conventions + # of format and FrozenBitset agree. + return FrozenBitset(format(x, format_style)[::-1], capacity=Q.dim()) def gen(self, i): """ @@ -918,7 +969,11 @@ def gen(self, i): sage: [Cl.gen(i) for i in range(3)] [x, y, z] """ +<<<<<<< HEAD return self._from_dict({FrozenBitset((i, )): self.base_ring().one()}, remove_zeros=False) +======= + return self._from_dict({FrozenBitset((i,)): self.base_ring().one()}, remove_zeros=False) +>>>>>>> a7cd8bf99e (Add _basis_index_function as a method, fix some doctests, and fix some type issues (tuple -> FrozenBitset)) def algebra_generators(self): """ @@ -1247,8 +1302,14 @@ def lift_module_morphism(self, m, names=None): Cl = CliffordAlgebra(Q, names) n = self._quadratic_form.dim() +<<<<<<< HEAD f = lambda x: self.prod(self._from_dict({FrozenBitset((j, )): m[j, i] for j in range(n)}, remove_zeros=True) for i in x) +======= + f = lambda x: self.prod(self._from_dict( {FrozenBitset((j,)): m[j,i] for j in range(n)}, + remove_zeros=True ) + for i in x) +>>>>>>> a7cd8bf99e (Add _basis_index_function as a method, fix some doctests, and fix some type issues (tuple -> FrozenBitset)) cat = AlgebrasWithBasis(self.category().base_ring()).Super().FiniteDimensional() return Cl.module_morphism(on_basis=f, codomain=self, category=cat) From 74d731f2c50e632789d372344e758e35042cda53 Mon Sep 17 00:00:00 2001 From: "Trevor K. Karn" Date: Thu, 23 Jun 2022 15:39:06 -0500 Subject: [PATCH 02/27] Add index class and first draft of exterior __mul__ --- src/sage/algebras/clifford_algebra.py | 86 +-------------------------- 1 file changed, 2 insertions(+), 84 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 6aa6b541724..612237c6e85 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -324,65 +324,6 @@ def conjugate(self): clifford_conjugate = conjugate - # TODO: This is a general function which should be moved to a - # superalgebras category when one is implemented. - def supercommutator(self, x): - r""" - Return the supercommutator of ``self`` and ``x``. - - Let `A` be a superalgebra. The *supercommutator* of homogeneous - elements `x, y \in A` is defined by - - .. MATH:: - - [x, y\} = x y - (-1)^{|x| |y|} y x - - and extended to all elements by linearity. - - EXAMPLES:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: a = x*y - z - sage: b = x - y + y*z - sage: a.supercommutator(b) - -5*x*y + 8*x*z - 2*y*z - 6*x + 12*y - 5*z - sage: a.supercommutator(Cl.one()) - 0 - sage: Cl.one().supercommutator(a) - 0 - sage: Cl.zero().supercommutator(a) - 0 - sage: a.supercommutator(Cl.zero()) - 0 - - sage: Q = QuadraticForm(ZZ, 2, [-1,1,-3]) - sage: Cl. = CliffordAlgebra(Q) - sage: [a.supercommutator(b) for a in Cl.basis() for b in Cl.basis()] - [0, 0, 0, 0, 0, -2, 1, -x - 2*y, 0, 1, - -6, 6*x + y, 0, x + 2*y, -6*x - y, 0] - sage: [a*b-b*a for a in Cl.basis() for b in Cl.basis()] - [0, 0, 0, 0, 0, 0, 2*x*y - 1, -x - 2*y, 0, - -2*x*y + 1, 0, 6*x + y, 0, x + 2*y, -6*x - y, 0] - - Exterior algebras inherit from Clifford algebras, so - supercommutators work as well. We verify the exterior algebra - is supercommutative:: - - sage: E. = ExteriorAlgebra(QQ) - sage: all(b1.supercommutator(b2) == 0 - ....: for b1 in E.basis() for b2 in E.basis()) - True - """ - P = self.parent() - ret = P.zero() - for ms, cs in self: - for mx, cx in x: - ret += P.term(ms, cs) * P.term(mx, cx) - s = (-1)**(P.degree_on_basis(ms) * P.degree_on_basis(mx)) - ret -= s * P.term(mx, cx) * P.term(ms, cs) - return ret - class CliffordAlgebraIndices(Parent): r""" @@ -400,7 +341,6 @@ def __call__(self, el): sage: E = ExteriorAlgebra(QQ, 7) sage: B = E.basis() """ - if not isinstance(el, Element): return self._element_constructor_(el) else: @@ -885,15 +825,9 @@ def _element_constructor_(self, x): if x in self.free_module(): R = self.base_ring() if x.parent().base_ring() is R: -<<<<<<< HEAD return self.element_class(self, {FrozenBitset((i, )): c for i, c in x.items()}) # if the base ring is different, attempt to coerce it into R return self.element_class(self, {FrozenBitset((i, )): R(c) for i, c in x.items() if R(c) != R.zero()}) -======= - return self.element_class(self, {FrozenBitset((i,)): c for i,c in x.items()}) - # if the base ring is different, attempt to coerce it into R - return self.element_class(self, {FrozenBitset((i,)): R(c) for i,c in x.items() if R(c) != R.zero()}) ->>>>>>> a7cd8bf99e (Add _basis_index_function as a method, fix some doctests, and fix some type issues (tuple -> FrozenBitset)) if (isinstance(x, CliffordAlgebraElement) and self.has_coerce_map_from(x.parent())): @@ -904,18 +838,10 @@ def _element_constructor_(self, x): R = self.base_ring() return self.element_class(self, {FrozenBitset((i,)): R.one() for i in x}) -<<<<<<< HEAD try: return super(CliffordAlgebra, self)._element_constructor_(x) except TypeError: raise TypeError(f'do not know how to make {x=} an element of self') -======= - if isinstance(x, tuple): - R = self.base_ring() - return self.element_class(self, {FrozenBitset((i,)): R.one() for i in x}) - - return super(CliffordAlgebra, self)._element_constructor_(x) ->>>>>>> a7cd8bf99e (Add _basis_index_function as a method, fix some doctests, and fix some type issues (tuple -> FrozenBitset)) def _basis_index_function(self, x): """ @@ -969,11 +895,7 @@ def gen(self, i): sage: [Cl.gen(i) for i in range(3)] [x, y, z] """ -<<<<<<< HEAD return self._from_dict({FrozenBitset((i, )): self.base_ring().one()}, remove_zeros=False) -======= - return self._from_dict({FrozenBitset((i,)): self.base_ring().one()}, remove_zeros=False) ->>>>>>> a7cd8bf99e (Add _basis_index_function as a method, fix some doctests, and fix some type issues (tuple -> FrozenBitset)) def algebra_generators(self): """ @@ -1302,14 +1224,8 @@ def lift_module_morphism(self, m, names=None): Cl = CliffordAlgebra(Q, names) n = self._quadratic_form.dim() -<<<<<<< HEAD f = lambda x: self.prod(self._from_dict({FrozenBitset((j, )): m[j, i] for j in range(n)}, remove_zeros=True) for i in x) -======= - f = lambda x: self.prod(self._from_dict( {FrozenBitset((j,)): m[j,i] for j in range(n)}, - remove_zeros=True ) - for i in x) ->>>>>>> a7cd8bf99e (Add _basis_index_function as a method, fix some doctests, and fix some type issues (tuple -> FrozenBitset)) cat = AlgebrasWithBasis(self.category().base_ring()).Super().FiniteDimensional() return Cl.module_morphism(on_basis=f, codomain=self, category=cat) @@ -1394,6 +1310,7 @@ def lift_isometry(self, m, names=None): Cl = CliffordAlgebra(Q, names) n = Q.dim() + f = lambda x: Cl.prod(Cl._from_dict({FrozenBitset((j, )): m[j, i] for j in range(n)}, remove_zeros=True) for i in x) cat = AlgebrasWithBasis(self.category().base_ring()).Super().FiniteDimensional() @@ -2252,6 +2169,7 @@ def _mul_(self, other): d = {} n = P.ngens() +<<<<<<< HEAD for ml, cl in self: # ml for "monomial on the left" for mr, cr in other: # mr for "monomial on the right" if ml.intersection(mr): From 4be416f78c285a86972e160d5bbb0958e805bb08 Mon Sep 17 00:00:00 2001 From: "Trevor K. Karn" Date: Tue, 5 Jul 2022 13:26:10 -0700 Subject: [PATCH 03/27] Add doctest --- src/sage/algebras/clifford_algebra.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 612237c6e85..49e6b8ccc97 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -2169,7 +2169,6 @@ def _mul_(self, other): d = {} n = P.ngens() -<<<<<<< HEAD for ml, cl in self: # ml for "monomial on the left" for mr, cr in other: # mr for "monomial on the right" if ml.intersection(mr): @@ -2186,8 +2185,9 @@ def _mul_(self, other): num_cross = 0 # keep track of the number of signs tot_cross = 0 for i in ml: + num_cross_new = 0 while i > j: - num_cross += 1 + num_cross_new += 1 try: j = next(it) except StopIteration: From 65984c5d2253024d6a5566ab4ce90194332f4299 Mon Sep 17 00:00:00 2001 From: "Trevor K. Karn" Date: Sat, 9 Jul 2022 07:23:25 -0500 Subject: [PATCH 04/27] Initial commit of f4 --- src/sage/algebras/clifford_algebra.py | 113 ++++++++++++++++++++++++++ 1 file changed, 113 insertions(+) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 49e6b8ccc97..df96ebfa4fb 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -2126,6 +2126,119 @@ def lifted_form(x, y): codomain=self.base_ring(), name="Bilinear Form") + def exterior_bitset_f4(self, I): + r""" + Return a Groebner basis for an ideal `I` of the exterior algebra. + + EXAMPLES:: + + sage: 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: exterior_bitset_f4(I) + (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) + + The example above was computed first using M2: + + E = QQ[a..e, SkewCommutative => true] + I = ideal( 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) + groebnerBasis(I) + + returns: + o3 = | bcd-bce+bde-cde acd-ace+ade-cde abd-abe+ade-bde abc-abe+ace-bce | + """ + + # following https://pi.math.cornell.edu/~djp282/documents/math6140-2017a.pdf + + def get_pairs(pairs): + # this implements the `select` function of Peifer + + # change pairs in here so I don't have to do it after calling get_pairs. + return pairs.pop() # temporarily do Buchbergers + + def f4_sel(P): + # this is the function on page 13 of the original F4 paper. + + return P + + # P is a list of pairs + # TODO: implement the better choice. + + def symbolic_preprocessing(P, G): + # the first/second terms of the S polynomial might + # have to be adjusted a la Stokes' paper + left = set() # the first term of the S polynomial + right = set() # the second term of the S polynomial + L = left.union(right) + done = set(f.monomials()[0] for f in L) + + from itertools import chain + mon_L = set(chain.from_iterable(f.monomials() for f in L)) + + while done != mon_L: + m = sorted(mon_L.difference(done), key = lambda x: (-len(x.support_of_term()), tuple(x.support_of_term())))[0] + done = done.add(m) + for g in G: + if divides(g.monomials()[0], m): + L.add(m * g/g.monomials()[0]) + break + return L + + def f4_reduce(P, G): + # given a current set of pairs P and a + # current basis G, return G' of new basis + + L = symbolic_preprocessing(P, G) + lm_L = set(f.monomials()[0] for f in L) + + d = dict() + for i, f in enumerate(F): + d.update(((i,hash(m)),c) for m,c in f._monomial_coefficients.items()) + M = MatrixSpace(E.base_ring(), len(L), 2^self.ngens(), sparse=True) + poly_matrix = M(d).rref() + + Lprime = set(E._from_dict(dict((FrozenBitset(format(k,'b')[::-1]),v) for k,v in row.dict().items())) for row in poly_matrix) + + Gprime = set() + + for f in Lprime: + if f.monomials()[0] in lm_L: + continue + Gprime.add(f) + + return Gprime + + F = I.gens() + G = set(F) + k = I.ngens() + + from itertools import combinations + pairs = set(combinations(range(k), 2)) # this is Peifer's P + + while pairs: + P = f4_sel(pairs) # this is different from Buchbergers which would be pairs.pop() + Gtemp = f4_reduce(P, G) + pairs.difference_update(P) + + for h in Gtemp: + G.add(h) + k += 1 + pairs.update((i,k) for i in range(k)) + + return G + + class Element(CliffordAlgebraElement): """ An element of an exterior algebra. From 7a79608e8b39fd52dc0dee2f563d03f6c06dd43f Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Tue, 12 Jul 2022 18:37:53 +0900 Subject: [PATCH 05/27] Another implementation of Groebner bases for the exterior algebra. --- src/sage/algebras/clifford_algebra.py | 87 +++++++ src/sage/algebras/exterior_algebra_cython.pxd | 20 ++ src/sage/algebras/exterior_algebra_cython.pyx | 245 ++++++++++++++++++ 3 files changed, 352 insertions(+) create mode 100644 src/sage/algebras/exterior_algebra_cython.pxd create mode 100644 src/sage/algebras/exterior_algebra_cython.pyx diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index df96ebfa4fb..7a72cd010c4 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/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. = 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. = 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) + diff --git a/src/sage/algebras/exterior_algebra_cython.pxd b/src/sage/algebras/exterior_algebra_cython.pxd new file mode 100644 index 00000000000..66c33994f23 --- /dev/null +++ b/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) + diff --git a/src/sage/algebras/exterior_algebra_cython.pyx b/src/sage/algebras/exterior_algebra_cython.pyx new file mode 100644 index 00000000000..e3c4dc64888 --- /dev/null +++ b/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() + 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 = 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: ( 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 = 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( 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( (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( (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]) + From c8a3766ef1e6b5d1eebcbd052772b2b61aca013e Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Tue, 12 Jul 2022 19:01:16 +0900 Subject: [PATCH 06/27] Fixing some little details. --- src/sage/algebras/clifford_algebra.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 7a72cd010c4..caed257405c 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -2336,8 +2336,8 @@ def reduce(self, I, left=True): - ``side`` -- the side, ignored if ``I`` is an ideal """ if isinstance(I, ExteriorAlgebraIdeal): - I = I.groebner_basis() left = (I.side() == "left") + I = I.groebner_basis() E = self.parent() f = self @@ -3209,7 +3209,7 @@ def reduce(self, f): """ Reduce ``f`` modulo ``self``. """ - return f.reduce(self.groebner_basis()) + return f.reduce(self) @cached_method def groebner_basis(self): From b62d5782d78a613b7b0664a895c00fb80298f14a Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Wed, 13 Jul 2022 14:56:23 +0900 Subject: [PATCH 07/27] Cythonizing the element classes. --- src/sage/algebras/clifford_algebra.py | 301 +----------------- src/sage/algebras/exterior_algebra_cython.pxd | 17 +- 2 files changed, 22 insertions(+), 296 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index caed257405c..91d4efc10e7 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -22,8 +22,8 @@ from sage.structure.parent import Parent from sage.structure.element import Element from sage.data_structures.bitset import Bitset, FrozenBitset -from copy import copy +from sage.algebras.clifford_algebra_element import CliffordAlgebraElement, ExteriorAlgebraElement from sage.categories.algebras_with_basis import AlgebrasWithBasis from sage.categories.hopf_algebras_with_basis import HopfAlgebrasWithBasis from sage.categories.finite_enumerated_sets import FiniteEnumeratedSets @@ -38,294 +38,12 @@ from sage.combinat.free_module import CombinatorialFreeModule from sage.combinat.subset import Subsets from sage.quadratic_forms.quadratic_form import QuadraticForm -from sage.algebras.weyl_algebra import repr_from_monomials from sage.misc.inherit_comparison import InheritComparisonClasscallMetaclass from sage.typeset.ascii_art import ascii_art from sage.typeset.unicode_art import unicode_art import unicodedata -class CliffordAlgebraElement(CombinatorialFreeModule.Element): - """ - An element in a Clifford algebra. - - TESTS:: - - sage: Q = QuadraticForm(ZZ, 3, [1, 2, 3, 4, 5, 6]) - sage: Cl. = CliffordAlgebra(Q) - sage: elt = ((x^3-z)*x + y)^2 - sage: TestSuite(elt).run() - """ - def _repr_(self): - """ - Return a string representation of ``self``. - - TESTS:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: ((x^3-z)*x + y)^2 - -2*x*y*z - x*z + 5*x - 4*y + 2*z + 2 - sage: Cl.zero() - 0 - """ - return repr_from_monomials(self.list(), self.parent()._repr_term) - - def _latex_(self): - r""" - Return a `\LaTeX` representation of ``self``. - - TESTS:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: latex( ((x^3-z)*x + y)^2 ) - -2 x y z - x z + 5 x - 4 y + 2 z + 2 - sage: Cl. = CliffordAlgebra(Q) - sage: latex( (x1 - x2)*x0 + 5*x0*x1*x2 ) - 5 x_{0} x_{1} x_{2} - x_{0} x_{1} + x_{0} x_{2} - 1 - """ - return repr_from_monomials(self.list(), - self.parent()._latex_term, - True) - - def _mul_(self, other): - """ - Return ``self`` multiplied by ``other``. - - INPUT: - - - ``other`` -- element of the same Clifford algebra as ``self`` - - EXAMPLES:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: (x^3 - z*y)*x*(y*z + x*y*z) - x*y*z + y*z - 24*x + 12*y + 2*z - 24 - sage: y*x - -x*y + 2 - sage: z*x - -x*z + 3 - sage: z*z - 6 - sage: x*0 - 0 - sage: 0*x - 0 - """ - Q = self.parent()._quadratic_form - zero = self.parent().base_ring().zero() - d = {} - - for ml, cl in self: - # Distribute the current term ``cl`` * ``ml`` over ``other``. - cur = copy(other._monomial_coefficients) # The current distribution of the term - for i in reversed(ml): - # Distribute the current factor ``e[i]`` (the ``i``-th - # element of the standard basis). - next = {} - # At the end of the following for-loop, ``next`` will be - # the dictionary describing the element - # ``e[i]`` * (the element described by the dictionary ``cur``) - # (where ``e[i]`` is the ``i``-th standard basis vector). - for mr, cr in cur.items(): - - # Commute the factor as necessary until we are in order - for j in mr: - if i <= j: - break - # Add the additional term from the commutation - # get a non-frozen bitset to manipulate - t = Bitset(mr) # a mutable copy - t.discard(j) - t = FrozenBitset(t) - next[t] = next.get(t, zero) + cr * Q[i, j] - # Note: ``Q[i,j] == Q(e[i]+e[j]) - Q(e[i]) - Q(e[j])`` for - # ``i != j``, where ``e[k]`` is the ``k``-th standard - # basis vector. - cr = -cr - if next[t] == zero: - del next[t] - - # Check to see if we have a squared term or not - mr = Bitset(mr) # temporarily mutable - if i in mr: - mr.discard(i) - cr *= Q[i, i] - # Note: ``Q[i,i] == Q(e[i])`` where ``e[i]`` is the - # ``i``-th standard basis vector. - else: - # mr is implicitly sorted - mr.add(i) - mr = FrozenBitset(mr) # refreeze it - next[mr] = next.get(mr, zero) + cr - if next[mr] == zero: - del next[mr] - cur = next - - # Add the distributed terms to the total - for index, coeff in cur.items(): - d[index] = d.get(index, zero) + cl * coeff - if d[index] == zero: - del d[index] - - return self.__class__(self.parent(), d) - - def list(self): - """ - Return the list of monomials and their coefficients in ``self`` - (as a list of `2`-tuples, each of which has the form - ``(monomial, coefficient)``). - - EXAMPLES:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: elt = 5*x + y - sage: elt.list() - [(1, 5), (01, 1)] - """ - return sorted(self._monomial_coefficients.items(), key=lambda m: (-len(m[0]), list(m[0]))) - - def support(self): - """ - Return the support of ``self``. - - This is the list of all monomials which appear with nonzero - coefficient in ``self``. - - EXAMPLES:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: elt = 5*x + y - sage: elt.support() - [1, 01] - """ - return sorted(self._monomial_coefficients.keys(), key=lambda x: (-len(x), list(x))) - - def reflection(self): - r""" - Return the image of the reflection automorphism on ``self``. - - The *reflection automorphism* of a Clifford algebra is defined - as the linear endomorphism of this algebra which maps - - .. MATH:: - - x_1 \wedge x_2 \wedge \cdots \wedge x_m \mapsto - (-1)^m x_1 \wedge x_2 \wedge \cdots \wedge x_m. - - It is an algebra automorphism of the Clifford algebra. - - :meth:`degree_negation` is an alias for :meth:`reflection`. - - EXAMPLES:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: elt = 5*x + y + x*z - sage: r = elt.reflection(); r - x*z - 5*x - y - sage: r.reflection() == elt - True - - TESTS: - - We check that the reflection is an involution:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: all(x.reflection().reflection() == x for x in Cl.basis()) - True - """ - return self.__class__(self.parent(), {m: (-1)**len(m) * c for m, c in self}) - - degree_negation = reflection - - def transpose(self): - r""" - Return the transpose of ``self``. - - The transpose is an anti-algebra involution of a Clifford algebra - and is defined (using linearity) by - - .. MATH:: - - x_1 \wedge x_2 \wedge \cdots \wedge x_m \mapsto - x_m \wedge \cdots \wedge x_2 \wedge x_1. - - EXAMPLES:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: elt = 5*x + y + x*z - sage: t = elt.transpose(); t - -x*z + 5*x + y + 3 - sage: t.transpose() == elt - True - sage: Cl.one().transpose() - 1 - - TESTS: - - We check that the transpose is an involution:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: all(x.transpose().transpose() == x for x in Cl.basis()) - True - - Zero is sent to zero:: - - sage: Cl.zero().transpose() == Cl.zero() - True - """ - P = self.parent() - if not self._monomial_coefficients: - return P.zero() - g = P.gens() - return P.sum(c * P.prod(g[i] for i in reversed(m)) for m, c in self) - - def conjugate(self): - r""" - Return the Clifford conjugate of ``self``. - - The Clifford conjugate of an element `x` of a Clifford algebra is - defined as - - .. MATH:: - - \bar{x} := \alpha(x^t) = \alpha(x)^t - - where `\alpha` denotes the :meth:`reflection ` - automorphism and `t` the :meth:`transposition `. - - EXAMPLES:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: elt = 5*x + y + x*z - sage: c = elt.conjugate(); c - -x*z - 5*x - y + 3 - sage: c.conjugate() == elt - True - - TESTS: - - We check that the conjugate is an involution:: - - sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) - sage: Cl. = CliffordAlgebra(Q) - sage: all(x.conjugate().conjugate() == x for x in Cl.basis()) - True - """ - return self.reflection().transpose() - - clifford_conjugate = conjugate - - class CliffordAlgebraIndices(Parent): r""" A facade parent for the indices of Clifford algebra. @@ -367,7 +85,7 @@ def __init__(self, Qdim): self._nbits = Qdim self._cardinality = 2**Qdim # the if statement here is in case Qdim is 0. - self._maximal_set = FrozenBitset('1'*self._nbits) if self._nbits else FrozenBitset('0') + self._maximal_set = FrozenBitset('1'*self._nbits) if self._nbits else FrozenBitset() category = FiniteEnumeratedSets().Facade() Parent.__init__(self, category=category, facade=True) @@ -446,7 +164,7 @@ def __iter__(self): """ import itertools n = self._nbits - yield FrozenBitset('0') + yield FrozenBitset() k = 1 while k <= n: for C in itertools.combinations(range(n), k): @@ -953,7 +671,7 @@ def one_basis(self): sage: Cl.one_basis() 0 """ - return FrozenBitset('0') + return FrozenBitset() def is_commutative(self): """ @@ -1908,7 +1626,11 @@ def coproduct_on_basis(self, a): one = self.base_ring().one() L = unshuffle_iterator(tuple(a), one) return self.tensor_square()._from_dict( +<<<<<<< HEAD {tuple(FrozenBitset(e) if e else FrozenBitset('0') for e in t): c for t, c in L if c}, +======= + {tuple(FrozenBitset(e) if e else FrozenBitset() for e in t): c for t,c in L if c}, +>>>>>>> b0f66e328e (Cythonizing the element classes.) coerce=False, remove_zeros=False) @@ -2247,9 +1969,11 @@ def _ideal_class_(self, n=0): sage: E. = ExteriorAlgebra(QQ) sage: E._ideal_class_() + """ return ExteriorAlgebraIdeal +<<<<<<< HEAD class Element(CliffordAlgebraElement): """ An element of an exterior algebra. @@ -2534,6 +2258,9 @@ def scalar(self, other): -2*x*y + 5*x*z + y*z """ return (self.transpose() * other).constant_coefficient() +======= + Element = ExteriorAlgebraElement +>>>>>>> b0f66e328e (Cythonizing the element classes.) ##################################################################### # Differentials @@ -2798,7 +2525,7 @@ def _on_basis(self, m): sage: E. = ExteriorAlgebra(QQ) sage: par = E.boundary({(0,1): z, (1,2): x, (2,0): y}) - sage: par._on_basis(FrozenBitset('0')) + sage: par._on_basis(FrozenBitset()) 0 sage: par._on_basis((0,)) 0 diff --git a/src/sage/algebras/exterior_algebra_cython.pxd b/src/sage/algebras/exterior_algebra_cython.pxd index 66c33994f23..022197c2911 100644 --- a/src/sage/algebras/exterior_algebra_cython.pxd +++ b/src/sage/algebras/exterior_algebra_cython.pxd @@ -5,16 +5,15 @@ 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) +cdef Integer bitset_to_int(FrozenBitset X) +cdef FrozenBitset int_to_bitset(Integer n) +cdef unsigned long degree(FrozenBitset X) +cdef 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) +cdef build_monomial(supp, E) +cdef partial_S_poly(f, g, E, int side) +cdef set preprocessing(list P, list G, E, int side) +cdef list reduction(list P, list G, E, int side) From d91a556f756dc7248382478ec27c5bd3986b9ba4 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Wed, 13 Jul 2022 16:09:10 +0900 Subject: [PATCH 08/27] Adding more type declarations and renaming the Groebner file. --- src/sage/algebras/clifford_algebra.py | 2 +- src/sage/algebras/exterior_algebra_cython.pxd | 19 ------ .../algebras/exterior_algebra_groebner.pxd | 21 +++++++ ...thon.pyx => exterior_algebra_groebner.pyx} | 62 ++++++++++++------- 4 files changed, 61 insertions(+), 43 deletions(-) delete mode 100644 src/sage/algebras/exterior_algebra_cython.pxd create mode 100644 src/sage/algebras/exterior_algebra_groebner.pxd rename src/sage/algebras/{exterior_algebra_cython.pyx => exterior_algebra_groebner.pyx} (77%) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 91d4efc10e7..38508538ef9 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -2960,7 +2960,7 @@ def groebner_basis(self): -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 + from sage.algebras.exterior_algebra_groebner import compute_groebner side = 2 if self.side() == "left": side = 0 diff --git a/src/sage/algebras/exterior_algebra_cython.pxd b/src/sage/algebras/exterior_algebra_cython.pxd deleted file mode 100644 index 022197c2911..00000000000 --- a/src/sage/algebras/exterior_algebra_cython.pxd +++ /dev/null @@ -1,19 +0,0 @@ -""" -Exterior algebras backend -""" - -from sage.data_structures.bitset cimport FrozenBitset -from sage.rings.integer cimport Integer - -cdef Integer bitset_to_int(FrozenBitset X) -cdef FrozenBitset int_to_bitset(Integer n) -cdef unsigned long degree(FrozenBitset X) -cdef FrozenBitset leading_supp(f) -cpdef tuple get_leading_supports(tuple I) - -# Grobner basis functions -cdef build_monomial(supp, E) -cdef partial_S_poly(f, g, E, int side) -cdef set preprocessing(list P, list G, E, int side) -cdef list reduction(list P, list G, E, int side) - diff --git a/src/sage/algebras/exterior_algebra_groebner.pxd b/src/sage/algebras/exterior_algebra_groebner.pxd new file mode 100644 index 00000000000..60916f57e0f --- /dev/null +++ b/src/sage/algebras/exterior_algebra_groebner.pxd @@ -0,0 +1,21 @@ +""" +Exterior algebras Gröbner bases +""" + +from sage.data_structures.bitset cimport FrozenBitset +from sage.rings.integer cimport Integer +from sage.algebras.clifford_algebra_element cimport CliffordAlgebraElement +from sage.structure.parent cimport Parent + +cdef Integer bitset_to_int(FrozenBitset X) +cdef FrozenBitset int_to_bitset(Integer n) +cdef long degree(FrozenBitset X) +cdef FrozenBitset leading_supp(CliffordAlgebraElement f) +cpdef tuple get_leading_supports(tuple I) + +# Grobner basis functions +cdef build_monomial(Parent E, FrozenBitset supp) +cdef partial_S_poly(CliffordAlgebraElement f, CliffordAlgebraElement g, Parent E, int side) +cdef set preprocessing(list P, list G, Parent E, int side) +cdef list reduction(list P, list G, Parent E, int side) + diff --git a/src/sage/algebras/exterior_algebra_cython.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx similarity index 77% rename from src/sage/algebras/exterior_algebra_cython.pyx rename to src/sage/algebras/exterior_algebra_groebner.pyx index e3c4dc64888..d0d9ecda1ad 100644 --- a/src/sage/algebras/exterior_algebra_cython.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -1,11 +1,28 @@ """ -Exterior algebras backend +Exterior algebras Gröbner bases -This contains the backend implementations in Cython for the exterior algebra. +This contains the backend implementations in Cython for the Gröbner bases +of exterior algebra. + +AUTHORS: + +- Trevor Karn, Travis Scrimshaw (July 2022): Initial implementation """ +#***************************************************************************** +# Copyright (C) 2022 Trevor Karn +# (C) 2022 Travis Scrimshaw +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# http://www.gnu.org/licenses/ +#***************************************************************************** + 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.data_structures.bitset_base cimport (bitset_t, bitset_init, bitset_first, + bitset_next, bitset_set_to, bitset_len) from sage.structure.parent cimport Parent cdef inline Integer bitset_to_int(FrozenBitset X): @@ -23,7 +40,7 @@ cdef inline FrozenBitset int_to_bitset(Integer n): """ Convert a nonnegative integer ``n`` to a :class:`FrozenBitset`. """ - cdef unsigned long i + cdef Py_ssize_t i if mpz_sgn(n.value) == 0: return FrozenBitset() @@ -36,20 +53,15 @@ cdef inline FrozenBitset int_to_bitset(Integer n): return ret -cdef inline unsigned long degree(FrozenBitset X): +cdef inline 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 + return bitset_len(X._bitset) #TODO: Bring the exterior algebra elements as a cdef class and make f know its type! -cdef inline FrozenBitset leading_supp(f): +cdef inline FrozenBitset leading_supp(CliffordAlgebraElement f): """ Return the leading support of the exterior algebra element ``f``. """ @@ -70,16 +82,17 @@ cpdef tuple get_leading_supports(tuple I): - ``I`` -- a tuple of elements of an exterior algebra """ # We filter out any elements that are 0 + cdef CliffordAlgebraElement f return tuple(set([leading_supp(f) for f in I if f._monomial_coefficients])) -cdef inline build_monomial(E, supp): +cdef inline build_monomial(Parent E, FrozenBitset supp): """ Helper function for the fastest way to build a monomial. """ - return E.element_class(E, {supp: ( E)._base.one()}) + return E.element_class(E, {supp: E._base.one()}) -cdef inline partial_S_poly(f, g, E, int side): +cdef inline partial_S_poly(CliffordAlgebraElement f, CliffordAlgebraElement g, Parent E, int side): """ Compute one half of the `S`-polynomial for ``f`` and ``g``. @@ -98,11 +111,12 @@ cdef inline partial_S_poly(f, g, E, int side): ret = f * build_monomial(E, D) return ret * (~ret[lmf._union(lmg)]) -cdef inline set preprocessing(list P, list G, E, int side): +cdef inline set preprocessing(list P, list G, Parent E, int side): """ Perform the preprocessing step. """ #print("Start preprocessing:", P) + cdef CliffordAlgebraElement f 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: @@ -123,26 +137,28 @@ cdef inline set preprocessing(list P, list G, E, int side): for g in G: lm = leading_supp(g) if lm <= m: - f = build_monomial(E, m.difference(lm)) * g + f = build_monomial(E, m.difference(lm)) * g if f in L: break - monL.update(set(f.support()) - done) + monL.update(set(f._monomial_coefficients) - done) L.add(f) break #print("preprocessing:", L) return L -cdef inline list reduction(list P, list G, E, int side): +cdef inline list reduction(list P, list G, Parent 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( m)): c for i,f in enumerate(L) for m,c in f._monomial_coefficients.items()}, + M = matrix({(i, bitset_to_int( 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) + 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] @@ -157,9 +173,9 @@ def compute_groebner(I, side): - ``side`` -- integer; the side of the ideal: ``0`` for left, ``1`` for right, and ``2`` for two-sided """ - E = I.ring() + cdef Parent E = I.ring() cdef FrozenBitset p0, p1 - cdef unsigned long deg + cdef 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 From 3a9d790c122cf90a3efd9bf86a0db04d2d3d3395 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Wed, 13 Jul 2022 16:38:32 +0900 Subject: [PATCH 09/27] Initial version of GB strategy file for exterior algebra. --- src/sage/algebras/clifford_algebra.py | 9 +- .../algebras/exterior_algebra_groebner.pxd | 26 +- .../algebras/exterior_algebra_groebner.pyx | 434 +++++++++--------- 3 files changed, 245 insertions(+), 224 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 38508538ef9..bb3387bc88f 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -2960,11 +2960,6 @@ def groebner_basis(self): -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_groebner import compute_groebner - side = 2 - if self.side() == "left": - side = 0 - elif self.side() == "right": - side = 1 - return compute_groebner(self, side) + from sage.algebras.exterior_algebra_groebner import GroebnerStrategy + return GroebnerStrategy(self).compute_groebner() diff --git a/src/sage/algebras/exterior_algebra_groebner.pxd b/src/sage/algebras/exterior_algebra_groebner.pxd index 60916f57e0f..a35f22ae36a 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pxd +++ b/src/sage/algebras/exterior_algebra_groebner.pxd @@ -6,16 +6,26 @@ from sage.data_structures.bitset cimport FrozenBitset from sage.rings.integer cimport Integer from sage.algebras.clifford_algebra_element cimport CliffordAlgebraElement from sage.structure.parent cimport Parent +from sage.structure.element cimport MonoidElement -cdef Integer bitset_to_int(FrozenBitset X) -cdef FrozenBitset int_to_bitset(Integer n) cdef long degree(FrozenBitset X) -cdef FrozenBitset leading_supp(CliffordAlgebraElement f) -cpdef tuple get_leading_supports(tuple I) +cdef CliffordAlgebraElement build_monomial(Parent E, FrozenBitset supp) # Grobner basis functions -cdef build_monomial(Parent E, FrozenBitset supp) -cdef partial_S_poly(CliffordAlgebraElement f, CliffordAlgebraElement g, Parent E, int side) -cdef set preprocessing(list P, list G, Parent E, int side) -cdef list reduction(list P, list G, Parent E, int side) +cdef class GroebnerStrategy: + cdef Parent E # the exterior algebra + cdef int side + cdef MonoidElement ideal + + cdef inline FrozenBitset leading_supp(self, CliffordAlgebraElement f) + cdef inline partial_S_poly_left(self, CliffordAlgebraElement f, CliffordAlgebraElement g) + cdef inline partial_S_poly_right(self, CliffordAlgebraElement f, CliffordAlgebraElement g) + cdef set preprocessing(self, list P, list G) + cdef list reduction(self, list P, list G) + + # These are the methods that determine the ordering of the monomials. + # Override these for other orderings. + # TODO: Make them abstract methods that must be implemented in subclasses + cdef inline Integer bitset_to_int(self, FrozenBitset X) + cdef inline FrozenBitset int_to_bitset(self, Integer n) diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index d0d9ecda1ad..df13f4f30b4 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -25,34 +25,6 @@ from sage.data_structures.bitset_base cimport (bitset_t, bitset_init, bitset_fir bitset_next, bitset_set_to, bitset_len) 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 Py_ssize_t i - - if mpz_sgn(n.value) == 0: - return FrozenBitset() - - cdef FrozenBitset ret = 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 long degree(FrozenBitset X): """ Compute the degree of ``X``. @@ -60,202 +32,246 @@ cdef inline long degree(FrozenBitset X): return bitset_len(X._bitset) -#TODO: Bring the exterior algebra elements as a cdef class and make f know its type! -cdef inline FrozenBitset leading_supp(CliffordAlgebraElement f): - """ - Return the leading support of the exterior algebra element ``f``. - """ - cdef dict mc = 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 - cdef CliffordAlgebraElement f - return tuple(set([leading_supp(f) for f in I if f._monomial_coefficients])) - - -cdef inline build_monomial(Parent E, FrozenBitset supp): +cdef inline CliffordAlgebraElement build_monomial(Parent E, FrozenBitset supp): """ Helper function for the fastest way to build a monomial. """ - return E.element_class(E, {supp: E._base.one()}) + return E.element_class(E, {supp: E._base.one()}) -cdef inline partial_S_poly(CliffordAlgebraElement f, CliffordAlgebraElement g, Parent E, int side): +cdef class GroebnerStrategy: """ - Compute one half of the `S`-polynomial for ``f`` and ``g``. - - This computes: + A strategy for computing a Gröbner basis. - .. MATH:: + INPUT: - LCM(LM(f), LM(g)) / LT(f) \cdot f. + - ``I`` -- the ideal """ - cdef FrozenBitset lmf = leading_supp(f) - cdef FrozenBitset lmg = leading_supp(g) - cdef FrozenBitset D = lmg.difference(lmf) - if side == 0: - ret = build_monomial(E, D) * f + def __init__(self, I): + """ + Initialize ``self``. + """ + self.ideal = I + self.E = I.ring() + if I.side() == "left": + self.side = 0 + elif I.side() == "right": + self.side = 1 + else: + self.side = 2 + + cdef inline FrozenBitset leading_supp(self, CliffordAlgebraElement f): + """ + Return the leading support of the exterior algebra element ``f``. + """ + cdef dict mc = f._monomial_coefficients + return self.int_to_bitset(min(self.bitset_to_int(k) for k in mc)) + + cdef inline partial_S_poly_left(self, CliffordAlgebraElement f, CliffordAlgebraElement g): + """ + 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 = self.leading_supp(f) + cdef FrozenBitset lmg = self.leading_supp(g) + cdef FrozenBitset D = lmg.difference(lmf) + ret = build_monomial(self.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, Parent E, int side): - """ - Perform the preprocessing step. - """ - #print("Start preprocessing:", P) - cdef CliffordAlgebraElement f - 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: + cdef inline partial_S_poly_right(self, CliffordAlgebraElement f, CliffordAlgebraElement g): + """ + Compute one half of the `S`-polynomial for ``f`` and ``g``. + + This computes: + + .. MATH:: + + f \cdot LCM(LM(f), LM(g)) / LT(f). + """ + cdef FrozenBitset lmf = self.leading_supp(f) + cdef FrozenBitset lmg = self.leading_supp(g) + cdef FrozenBitset D = lmg.difference(lmf) + ret = f * build_monomial(self.E, D) + return ret * (~ret[lmf._union(lmg)]) + + cdef inline set preprocessing(self, list P, list G): + """ + Perform the preprocessing step. + """ + #print("Start preprocessing:", P) + cdef CliffordAlgebraElement f, g + + cdef set L = set() + if self.side == 0: + L.update(self.partial_S_poly_left(f0, f1) for f0,f1 in P) + L.update(self.partial_S_poly_left(f1, f0) for f0,f1 in P) + elif self.side == 1: + L.update(self.partial_S_poly_right(f0, f1) for f0,f1 in P) + L.update(self.partial_S_poly_right(f1, f0) for f0,f1 in P) + if self.side == 2: + L.update(self.partial_S_poly_left(f0, f1) for f0,f1 in P) + L.update(self.partial_S_poly_left(f1, f0) for f0,f1 in P) + L.update(self.partial_S_poly_right(f0, f1) for f0,f1 in P) + L.update(self.partial_S_poly_right(f1, f0) for f0,f1 in P) + + cdef set done = set(self.leading_supp(f) for f in L) + cdef set monL = set() + for f in L: + monL.update(f._monomial_coefficients) + monL.difference_update(done) + + while monL: + m = self.int_to_bitset(max(self.bitset_to_int(k) for k in monL)) + done.add(m) + monL.remove(m) + for g in G: + lm = self.leading_supp(g) + if lm <= m: + f = build_monomial(self.E, m.difference(lm)) * g + if f in L: + break + monL.update(set(f._monomial_coefficients) - done) + L.add(f) break - monL.update(set(f._monomial_coefficients) - done) - L.add(f) - break - #print("preprocessing:", L) - return L - -cdef inline list reduction(list P, list G, Parent 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( 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 - """ - cdef Parent E = I.ring() - cdef FrozenBitset p0, p1 - cdef 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( (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) + #print("preprocessing:", L) + return L + + cdef inline list reduction(self, list P, list G): + """ + Perform the reduction of ``P`` mod ``G`` in ``E``. + """ + cdef set L = self.preprocessing(P, G) + cdef Py_ssize_t i + from sage.matrix.constructor import matrix + M = matrix({(i, self.bitset_to_int( 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(self.leading_supp( f) for f in L) + return [self.E.element_class(self.E, {self.int_to_bitset(Integer(j)): c for j,c in M[i].iteritems()}) + for i,p in enumerate(M.pivots()) + if self.int_to_bitset(Integer(p)) not in lead_supports] + + def compute_groebner(self): + """ + Compute the reduced ``side`` Gröbner basis for the ideal ``I``. + """ + cdef FrozenBitset p0, p1 + cdef long deg + cdef Py_ssize_t i, j, k + cdef list G = [f for f in self.ideal.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 = self.leading_supp(f0) + for j in range(i+1, n): + f1 = G[j] + p1 = self.leading_supp(f1) deg = degree( (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 + + while P: + #print("Cur G:", G) + Pp = P.pop(min(P)) # The selection: lowest lcm degree + Gp = self.reduction(Pp, G) + #print("Reduction yielded:", Gp) + G.extend(Gp) + for j in range(n, len(G)): + f1 = G[j] + p1 = self.leading_supp(f1) + for i in range(j): + f0 = G[i] + p0 = self.leading_supp(f0) + deg = degree( (p0._union(p1))) + if deg in P: + P[deg].append((f0, f1)) 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]) + 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 = self.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 = self.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]) + + + cdef inline Integer bitset_to_int(self, 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(self, Integer n): + """ + Convert a nonnegative integer ``n`` to a :class:`FrozenBitset`. + """ + cdef Py_ssize_t i + + if mpz_sgn(n.value) == 0: + return FrozenBitset() + + cdef FrozenBitset ret = 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 + From 368491cc191b499916b48b38ae58e44432b7f4a0 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Thu, 14 Jul 2022 12:03:52 +0900 Subject: [PATCH 10/27] Implementing different term orders for GB. --- src/sage/algebras/clifford_algebra.py | 53 +- .../algebras/clifford_algebra_element.pxd | 12 + .../algebras/clifford_algebra_element.pyx | 595 ++++++++++++++++++ .../algebras/exterior_algebra_groebner.pxd | 19 +- .../algebras/exterior_algebra_groebner.pyx | 235 ++++++- 5 files changed, 886 insertions(+), 28 deletions(-) create mode 100644 src/sage/algebras/clifford_algebra_element.pxd create mode 100644 src/sage/algebras/clifford_algebra_element.pyx diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index bb3387bc88f..e4f211c5512 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -18,6 +18,7 @@ #***************************************************************************** from sage.misc.cachefunc import cached_method +from sage.misc.lazy_attribute import lazy_attribute from sage.structure.unique_representation import UniqueRepresentation from sage.structure.parent import Parent from sage.structure.element import Element @@ -2932,17 +2933,36 @@ class ExteriorAlgebraIdeal(Ideal_nc): """ An ideal of the exterior algebra. """ + def __init__(self, ring, gens, coerce=True, side="twosided"): + """ + Initialize ``self``. + """ + self._groebner_basis = None + self._groebner_strategy = None + self._homogeneous = all(x.is_super_homogeneous() for x in gens) + if self._homogeneous: + side = "twosided" + Ideal_nc.__init__(self, ring, gens, coerce, side) + def reduce(self, f): """ Reduce ``f`` modulo ``self``. """ return f.reduce(self) - @cached_method - def groebner_basis(self): + def groebner_basis(self, term_order="negrevlex"): r""" Return the reduced Gröbner basis of ``self``. + INPUT: + + - ``term_order`` -- the term order used to compute the Gröbner basis; + must be one of the following: + + * ``"negrevlex"`` -- (default) negative reverse lex order + * ``"degrevlex"`` -- degree reverse lex order + * ``"deglex"`` -- degree lex order + EXAMPLES: We compute an example that was checked against Macaulay2:: @@ -2955,11 +2975,36 @@ def groebner_basis(self): ....: b*c*d - a*c*d + a*b*d - a*b*c] sage: I = E.ideal(rels) sage: I.groebner_basis() + (-a*c*d + a*c*e - a*d*e + c*d*e, + -a*b*c + a*b*d - a*c*d + b*c*d, + -a*b*d + a*b*e - a*d*e + b*d*e, + -a*b*c + a*b*e - a*c*e + b*c*e) + + With different term orders:: + + sage: I.groebner_basis("degrevlex") (-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) + + sage: I.groebner_basis("deglex") + (-a*c*d + a*c*e - a*d*e + c*d*e, + -a*b*c + a*b*d - a*c*d + b*c*d, + -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_groebner import GroebnerStrategy - return GroebnerStrategy(self).compute_groebner() + if term_order == "negrevlex": + from sage.algebras.exterior_algebra_groebner import GroebnerStrategyNegRevLex as strategy + elif term_order == "degrevlex": + from sage.algebras.exterior_algebra_groebner import GroebnerStrategyDegRevLex as strategy + elif term_order == "deglex": + from sage.algebras.exterior_algebra_groebner import GroebnerStrategyDegLex as strategy + else: + raise ValueError("invalid term order") + if strategy == self._groebner_strategy: + return self._groebner_basis + self._groebner_strategy = strategy + self._groebner_basis = self._groebner_strategy(self).compute_groebner() + return self._groebner_basis diff --git a/src/sage/algebras/clifford_algebra_element.pxd b/src/sage/algebras/clifford_algebra_element.pxd new file mode 100644 index 00000000000..f7a25582d54 --- /dev/null +++ b/src/sage/algebras/clifford_algebra_element.pxd @@ -0,0 +1,12 @@ +""" +Clifford algebra elements +""" + +from sage.modules.with_basis.indexed_element cimport IndexedFreeModuleElement + +cdef class CliffordAlgebraElement(IndexedFreeModuleElement): + pass + +cdef class ExteriorAlgebraElement(CliffordAlgebraElement): + pass + diff --git a/src/sage/algebras/clifford_algebra_element.pyx b/src/sage/algebras/clifford_algebra_element.pyx new file mode 100644 index 00000000000..93568ebd1d1 --- /dev/null +++ b/src/sage/algebras/clifford_algebra_element.pyx @@ -0,0 +1,595 @@ +""" +Clifford algebra elements + +AUTHORS: + +- Travis Scrimshaw (2013-09-06): Initial version +- Trevor Karn (2022-07-10): Rewrite multiplication using bitsets +""" + +#***************************************************************************** +# Copyright (C) 2022 Trevor Karn +# (C) 2022 Travis Scrimshaw +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# http://www.gnu.org/licenses/ +#***************************************************************************** + +from sage.structure.parent cimport Parent +from sage.data_structures.bitset cimport FrozenBitset, Bitset +from sage.algebras.weyl_algebra import repr_from_monomials +from copy import copy + +cdef class CliffordAlgebraElement(IndexedFreeModuleElement): + """ + An element in a Clifford algebra. + + TESTS:: + + sage: Q = QuadraticForm(ZZ, 3, [1, 2, 3, 4, 5, 6]) + sage: Cl. = CliffordAlgebra(Q) + sage: elt = ((x^3-z)*x + y)^2 + sage: TestSuite(elt).run() + """ + def _repr_(self): + """ + Return a string representation of ``self``. + + TESTS:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl. = CliffordAlgebra(Q) + sage: ((x^3-z)*x + y)^2 + -2*x*y*z - x*z + 5*x - 4*y + 2*z + 2 + sage: Cl.zero() + 0 + """ + return repr_from_monomials(self.list(), self._parent._repr_term) + + def _latex_(self): + r""" + Return a `\LaTeX` representation of ``self``. + + TESTS:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl. = CliffordAlgebra(Q) + sage: latex( ((x^3-z)*x + y)^2 ) + -2 x y z - x z + 5 x - 4 y + 2 z + 2 + sage: Cl. = CliffordAlgebra(Q) + sage: latex( (x1 - x2)*x0 + 5*x0*x1*x2 ) + 5 x_{0} x_{1} x_{2} - x_{0} x_{1} + x_{0} x_{2} - 1 + """ + return repr_from_monomials(self.list(), self._parent._latex_term, True) + + cdef _mul_(self, other): + """ + Return ``self`` multiplied by ``other``. + + INPUT: + + - ``other`` -- element of the same Clifford algebra as ``self`` + + EXAMPLES:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl. = CliffordAlgebra(Q) + sage: (x^3 - z*y)*x*(y*z + x*y*z) + x*y*z + y*z - 24*x + 12*y + 2*z - 24 + sage: y*x + -x*y + 2 + sage: z*x + -x*z + 3 + sage: z*z + 6 + sage: x*0 + 0 + sage: 0*x + 0 + """ + Q = self._parent._quadratic_form + zero = self._parent._base.zero() + cdef dict next_level, cur, d = {} + cdef FrozenBitset ml, mr, t + cdef Py_ssize_t i, j + + for ml,cl in self: + # Distribute the current term ``cl`` * ``ml`` over ``other``. + cur = copy(other._monomial_coefficients) # The current distribution of the term + for i in reversed(ml): + # Distribute the current factor ``e[i]`` (the ``i``-th + # element of the standard basis). + next_level = {} + # At the end of the following for-loop, ``next`` will be + # the dictionary describing the element + # ``e[i]`` * (the element described by the dictionary ``cur``) + # (where ``e[i]`` is the ``i``-th standard basis vector). + for mr,cr in cur.items(): + + # Commute the factor as necessary until we are in order + for j in mr: + if i <= j: + break + # Add the additional term from the commutation + # get a non-frozen bitset to manipulate + t = Bitset(mr) # a mutable copy + t.discard(j) + t = FrozenBitset(t) + next_level[t] = next_level.get(t, zero) + cr * Q[i,j] + # Note: ``Q[i,j] == Q(e[i]+e[j]) - Q(e[i]) - Q(e[j])`` for + # ``i != j``, where ``e[k]`` is the ``k``-th standard + # basis vector. + cr = -cr + if next_level[t] == zero: + del next_level[t] + + # Check to see if we have a squared term or not + mr = Bitset(mr) # temporarily mutable + if i in mr: + mr.discard(i) + cr *= Q[i,i] + # Note: ``Q[i,i] == Q(e[i])`` where ``e[i]`` is the + # ``i``-th standard basis vector. + else: + # mr is implicitly sorted + mr.add(i) + mr = FrozenBitset(mr) # refreeze it + next_level[mr] = next_level.get(mr, zero) + cr + if next_level[mr] == zero: + del next_level[mr] + cur = next_level + + # Add the distributed terms to the total + for index,coeff in cur.items(): + d[index] = d.get(index, zero) + cl * coeff + if d[index] == zero: + del d[index] + + return self.__class__(self.parent(), d) + + def list(self): + """ + Return the list of monomials and their coefficients in ``self`` + (as a list of `2`-tuples, each of which has the form + ``(monomial, coefficient)``). + + EXAMPLES:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl. = CliffordAlgebra(Q) + sage: elt = 5*x + y + sage: elt.list() + [(1, 5), (01, 1)] + """ + return sorted(self._monomial_coefficients.items(), key=lambda m: (-len(m[0]), list(m[0]))) + + def support(self): + """ + Return the support of ``self``. + + This is the list of all monomials which appear with nonzero + coefficient in ``self``. + + EXAMPLES:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl. = CliffordAlgebra(Q) + sage: elt = 5*x + y + sage: elt.support() + [1, 01] + """ + return sorted(self._monomial_coefficients, key=lambda x: (-len(x), list(x))) + + def reflection(self): + r""" + Return the image of the reflection automorphism on ``self``. + + The *reflection automorphism* of a Clifford algebra is defined + as the linear endomorphism of this algebra which maps + + .. MATH:: + + x_1 \wedge x_2 \wedge \cdots \wedge x_m \mapsto + (-1)^m x_1 \wedge x_2 \wedge \cdots \wedge x_m. + + It is an algebra automorphism of the Clifford algebra. + + :meth:`degree_negation` is an alias for :meth:`reflection`. + + EXAMPLES:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl. = CliffordAlgebra(Q) + sage: elt = 5*x + y + x*z + sage: r = elt.reflection(); r + x*z - 5*x - y + sage: r.reflection() == elt + True + + TESTS: + + We check that the reflection is an involution:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl. = CliffordAlgebra(Q) + sage: all(x.reflection().reflection() == x for x in Cl.basis()) + True + """ + return self.__class__(self._parent, {m: (-1)**len(m) * c for m,c in self}) + + degree_negation = reflection + + def transpose(self): + r""" + Return the transpose of ``self``. + + The transpose is an anti-algebra involution of a Clifford algebra + and is defined (using linearity) by + + .. MATH:: + + x_1 \wedge x_2 \wedge \cdots \wedge x_m \mapsto + x_m \wedge \cdots \wedge x_2 \wedge x_1. + + EXAMPLES:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl. = CliffordAlgebra(Q) + sage: elt = 5*x + y + x*z + sage: t = elt.transpose(); t + -x*z + 5*x + y + 3 + sage: t.transpose() == elt + True + sage: Cl.one().transpose() + 1 + + TESTS: + + We check that the transpose is an involution:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl. = CliffordAlgebra(Q) + sage: all(x.transpose().transpose() == x for x in Cl.basis()) + True + + Zero is sent to zero:: + + sage: Cl.zero().transpose() == Cl.zero() + True + """ + P = self._parent + if not self._monomial_coefficients: + return P.zero() + g = P.gens() + return P.sum(c * P.prod(g[i] for i in reversed(m)) for m,c in self) + + def conjugate(self): + r""" + Return the Clifford conjugate of ``self``. + + The Clifford conjugate of an element `x` of a Clifford algebra is + defined as + + .. MATH:: + + \bar{x} := \alpha(x^t) = \alpha(x)^t + + where `\alpha` denotes the :meth:`reflection ` + automorphism and `t` the :meth:`transposition `. + + EXAMPLES:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl. = CliffordAlgebra(Q) + sage: elt = 5*x + y + x*z + sage: c = elt.conjugate(); c + -x*z - 5*x - y + 3 + sage: c.conjugate() == elt + True + + TESTS: + + We check that the conjugate is an involution:: + + sage: Q = QuadraticForm(ZZ, 3, [1,2,3,4,5,6]) + sage: Cl. = CliffordAlgebra(Q) + sage: all(x.conjugate().conjugate() == x for x in Cl.basis()) + True + """ + return self.reflection().transpose() + + clifford_conjugate = conjugate + + +cdef class ExteriorAlgebraElement(CliffordAlgebraElement): + """ + An element of an exterior algebra. + """ + cdef _mul_(self, other): + """ + Return ``self`` multiplied by ``other``. + + INPUT: + + - ``other`` -- element of the same exterior algebra as ``self`` + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: x*y + x*y + sage: y*x + -x*y + sage: z*y*x + -x*y*z + sage: (x*z)*y + -x*y*z + sage: (3*x + y)^2 + 0 + sage: (x - 3*y + z/3)^2 + 0 + sage: (x+y) * (y+z) + x*y + x*z + y*z + + sage: E. = ExteriorAlgebra(QQ) + sage: (x * y) * (w * z) + -x*y*z*w + sage: x * y * w * z + -x*y*z*w + sage: (z * w) * (x * y) + x*y*z*w + """ + cdef Parent P = self._parent + zero = P._base.zero() + cdef dict d = {} + cdef Py_ssize_t n = P.ngens() + cdef ExteriorAlgebraElement rhs = other + + cdef FrozenBitset ml, mr, t + cdef Py_ssize_t num_cross, tot_cross, i, j + + for ml,cl in self._monomial_coefficients.items(): # ml for "monomial on the left" + for mr,cr in rhs._monomial_coefficients.items(): # mr for "monomial on the right" + if ml.intersection(mr): + # if they intersect nontrivially, move along. + continue + + if not mr: + t = ml + else: + t = ml._union(mr) + it = iter(mr) + j = next(it) + + num_cross = 0 # keep track of the number of signs + tot_cross = 0 + for i in ml: + while i > j: + num_cross += 1 + try: + j = next(it) + except StopIteration: + j = n + 1 + tot_cross += num_cross + if tot_cross % 2: + cr = -cr + + d[t] = d.get(t, zero) + cl * cr + if not d[t]: + del d[t] + + 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 + """ + from sage.algebras.clifford_algebra import ExteriorAlgebraIdeal + if isinstance(I, ExteriorAlgebraIdeal): + left = (I.side() == "left") + I = I.groebner_basis() + + f = self + E = self._parent + 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 + ``self`` with respect to ``x`` (that is, the element + `\iota_{x}(\text{self})` of the exterior algebra). + + If `V` is an `R`-module, and if `\alpha` is a fixed element of + `V^*`, then the *interior product* with respect to `\alpha` is + an `R`-linear map + `i_{\alpha} \colon \Lambda(V) \to \Lambda(V)`, determined by + the following requirements: + + - `i_{\alpha}(v) = \alpha(v)` for all `v \in V = \Lambda^1(V)`, + - it is a graded derivation of degree `-1`: all `x` and `y` + in `\Lambda(V)` satisfy + + .. MATH:: + + i_{\alpha}(x \wedge y) = (i_{\alpha} x) \wedge y + + (-1)^{\deg x} x \wedge (i_{\alpha} y). + + It can be shown that this map `i_{\alpha}` is graded of + degree `-1` (that is, sends `\Lambda^k(V)` into + `\Lambda^{k-1}(V)` for every `k`). + + When `V` is a finite free `R`-module, the interior product can + also be defined by + + .. MATH:: + + (i_{\alpha} \omega)(u_1, \ldots, u_k) + = \omega(\alpha, u_1, \ldots, u_k), + + where `\omega \in \Lambda^k(V)` is thought of as an + alternating multilinear mapping from + `V^* \times \cdots \times V^*` to `R`. + + Since Sage is only dealing with exterior powers of modules + of the form `R^d` for some nonnegative integer `d`, the + element `\alpha \in V^*` can be thought of as an element of + `V` (by identifying the standard basis of `V = R^d` with its + dual basis). This is how `\alpha` should be passed to this + method. + + We then extend the interior product to all + `\alpha \in \Lambda (V^*)` by + + .. MATH:: + + i_{\beta \wedge \gamma} = i_{\gamma} \circ i_{\beta}. + + INPUT: + + - ``x`` -- element of (or coercing into) `\Lambda^1(V)` + (for example, an element of `V`); this plays the role of + `\alpha` in the above definition + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: x.interior_product(x) + 1 + sage: (x + x*y).interior_product(2*y) + -2*x + sage: (x*z + x*y*z).interior_product(2*y - x) + -2*x*z - y*z - z + sage: x.interior_product(E.one()) + x + sage: E.one().interior_product(x) + 0 + sage: x.interior_product(E.zero()) + 0 + sage: E.zero().interior_product(x) + 0 + + REFERENCES: + + - :wikipedia:`Exterior_algebra#Interior_product` + """ + P = self._parent + return P.sum([c * cx * P.interior_product_on_basis(m, mx) + for m,c in self for mx,cx in x]) + + antiderivation = interior_product + + def hodge_dual(self): + r""" + Return the Hodge dual of ``self``. + + The Hodge dual of an element `\alpha` of the exterior algebra is + defined as `i_{\alpha} \sigma`, where `\sigma` is the volume + form + (:meth:`~sage.algebras.clifford_algebra.ExteriorAlgebra.volume_form`) + and `i_{\alpha}` denotes the antiderivation function with + respect to `\alpha` (see :meth:`interior_product` for the + definition of this). + + .. NOTE:: + + The Hodge dual of the Hodge dual of a homogeneous element + `p` of `\Lambda(V)` equals `(-1)^{k(n-k)} p`, where + `n = \dim V` and `k = \deg(p) = |p|`. + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: x.hodge_dual() + y*z + sage: (x*z).hodge_dual() + -y + sage: (x*y*z).hodge_dual() + 1 + sage: [a.hodge_dual().hodge_dual() for a in E.basis()] + [1, x, y, z, x*y, x*z, y*z, x*y*z] + sage: (x + x*y).hodge_dual() + y*z + z + sage: (x*z + x*y*z).hodge_dual() + -y + 1 + sage: E = ExteriorAlgebra(QQ, 'wxyz') + sage: [a.hodge_dual().hodge_dual() for a in E.basis()] + [1, -w, -x, -y, -z, w*x, w*y, w*z, x*y, x*z, y*z, + -w*x*y, -w*x*z, -w*y*z, -x*y*z, w*x*y*z] + """ + volume_form = self._parent.volume_form() + return volume_form.interior_product(self) + + def constant_coefficient(self): + """ + Return the constant coefficient of ``self``. + + .. TODO:: + + Define a similar method for general Clifford algebras once + the morphism to exterior algebras is implemented. + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: elt = 5*x + y + x*z + 10 + sage: elt.constant_coefficient() + 10 + sage: x.constant_coefficient() + 0 + """ + return self._monomial_coefficients.get(self._parent.one_basis(), + self._parent._base.zero()) + + def scalar(self, other): + r""" + Return the standard scalar product of ``self`` with ``other``. + + The standard scalar product of `x, y \in \Lambda(V)` is + defined by `\langle x, y \rangle = \langle x^t y \rangle`, where + `\langle a \rangle` denotes the degree-0 term of `a`, and where + `x^t` denotes the transpose + (:meth:`~sage.algebras.clifford_algebra.CliffordAlgebraElement.transpose`) + of `x`. + + .. TODO:: + + Define a similar method for general Clifford algebras once + the morphism to exterior algebras is implemented. + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: elt = 5*x + y + x*z + sage: elt.scalar(z + 2*x) + 0 + sage: elt.transpose() * (z + 2*x) + -2*x*y + 5*x*z + y*z + """ + return (self.transpose() * other).constant_coefficient() + diff --git a/src/sage/algebras/exterior_algebra_groebner.pxd b/src/sage/algebras/exterior_algebra_groebner.pxd index a35f22ae36a..5d9430a1cce 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pxd +++ b/src/sage/algebras/exterior_algebra_groebner.pxd @@ -16,6 +16,9 @@ cdef class GroebnerStrategy: cdef Parent E # the exterior algebra cdef int side cdef MonoidElement ideal + cdef bint homogeneous + + cdef inline bint build_S_poly(self, CliffordAlgebraElement f, CliffordAlgebraElement g) cdef inline FrozenBitset leading_supp(self, CliffordAlgebraElement f) cdef inline partial_S_poly_left(self, CliffordAlgebraElement f, CliffordAlgebraElement g) @@ -24,8 +27,16 @@ cdef class GroebnerStrategy: cdef list reduction(self, list P, list G) # These are the methods that determine the ordering of the monomials. - # Override these for other orderings. - # TODO: Make them abstract methods that must be implemented in subclasses - cdef inline Integer bitset_to_int(self, FrozenBitset X) - cdef inline FrozenBitset int_to_bitset(self, Integer n) + # These must be implemented in subclasses. Declare them as "inline" there. + cdef Integer bitset_to_int(self, FrozenBitset X) + cdef FrozenBitset int_to_bitset(self, Integer n) + +cdef class GroebnerStrategyNegRevLex(GroebnerStrategy): + pass + +cdef class GroebnerStrategyDegRevLex(GroebnerStrategy): + cdef Integer rank + +cdef class GroebnerStrategyDegLex(GroebnerStrategy): + cdef Integer rank diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index df13f4f30b4..29998694656 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -52,7 +52,8 @@ cdef class GroebnerStrategy: """ self.ideal = I self.E = I.ring() - if I.side() == "left": + self.homogeneous = all(x.is_super_homogeneous() for x in I.gens()) + if self.homogeneous or I.side() == "left": self.side = 0 elif I.side() == "right": self.side = 1 @@ -64,7 +65,7 @@ cdef class GroebnerStrategy: Return the leading support of the exterior algebra element ``f``. """ cdef dict mc = f._monomial_coefficients - return self.int_to_bitset(min(self.bitset_to_int(k) for k in mc)) + return self.int_to_bitset(max(self.bitset_to_int(k) for k in mc)) cdef inline partial_S_poly_left(self, CliffordAlgebraElement f, CliffordAlgebraElement g): """ @@ -98,25 +99,46 @@ cdef class GroebnerStrategy: ret = f * build_monomial(self.E, D) return ret * (~ret[lmf._union(lmg)]) + cdef inline bint build_S_poly(self, CliffordAlgebraElement f, CliffordAlgebraElement g): + r""" + Check to see if we should build the `S`-polynomial. + + For homogeneous ideals, we throw out all those pairs `(f, g)` such that + + .. MATH:: + + LCM(LM(f), LM(g)) == LM(f) \cdot LM(g). + """ + if not self.homogeneous: + return True + + return ( self.leading_supp(f).intersection(self.leading_supp(g))).isempty() + cdef inline set preprocessing(self, list P, list G): """ Perform the preprocessing step. """ #print("Start preprocessing:", P) - cdef CliffordAlgebraElement f, g + cdef CliffordAlgebraElement f, g, f0, f1 cdef set L = set() if self.side == 0: - L.update(self.partial_S_poly_left(f0, f1) for f0,f1 in P) - L.update(self.partial_S_poly_left(f1, f0) for f0,f1 in P) + for f0, f1 in P: + if self.build_S_poly(f0, f1): + L.add(self.partial_S_poly_left(f0, f1)) + L.add(self.partial_S_poly_left(f1, f0)) elif self.side == 1: - L.update(self.partial_S_poly_right(f0, f1) for f0,f1 in P) - L.update(self.partial_S_poly_right(f1, f0) for f0,f1 in P) + for f0, f1 in P: + if self.build_S_poly(f0, f1): + L.add(self.partial_S_poly_right(f0, f1) for f0,f1 in P) + L.add(self.partial_S_poly_right(f1, f0) for f0,f1 in P) if self.side == 2: - L.update(self.partial_S_poly_left(f0, f1) for f0,f1 in P) - L.update(self.partial_S_poly_left(f1, f0) for f0,f1 in P) - L.update(self.partial_S_poly_right(f0, f1) for f0,f1 in P) - L.update(self.partial_S_poly_right(f1, f0) for f0,f1 in P) + for f0, f1 in P: + if self.build_S_poly(f0, f1): + L.add(self.partial_S_poly_left(f0, f1)) + L.add(self.partial_S_poly_left(f1, f0)) + L.add(self.partial_S_poly_right(f0, f1)) + L.add(self.partial_S_poly_right(f1, f0)) cdef set done = set(self.leading_supp(f) for f in L) cdef set monL = set() @@ -164,18 +186,12 @@ cdef class GroebnerStrategy: cdef FrozenBitset p0, p1 cdef long deg cdef Py_ssize_t i, j, k - cdef list G = [f for f in self.ideal.gens() if f] # Remove 0s TODO: We should make this unnecessary here + cdef list G = [f for f in self.ideal.gens() if f] # Remove 0s 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 = self.leading_supp(f0) @@ -229,7 +245,7 @@ cdef class GroebnerStrategy: if lm <= s: did_reduction = True mon = self.E.monomial(s - lm) - if side == 0: + if self.side == 0: gp = mon * g f = f - f[s] / gp[s] * gp else: @@ -246,7 +262,88 @@ cdef class GroebnerStrategy: return tuple([f for f in G if f]) + cdef Integer bitset_to_int(self, FrozenBitset X): + raise NotImplementedError + + cdef FrozenBitset int_to_bitset(self, Integer n): + raise NotImplementedError + def sorted_monomials(self, as_dict=False): + """ + Helper function to display the monomials in their term order + used by ``self``. + + EXAMPLES:: + + sage: from sage.algebras.exterior_algebra_groebner import * + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([x, y]) + sage: GroebnerStrategyNegRevLex(I).sorted_monomials() + [1, x, y, x*y, z, x*z, y*z, x*y*z] + sage: GroebnerStrategyDegLex(I).sorted_monomials() + [1, x, y, z, x*y, x*z, y*z, x*y*z] + sage: GroebnerStrategyDegRevLex(I).sorted_monomials() + [1, z, y, x, y*z, x*z, x*y, x*y*z] + + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([a, b]) + sage: GroebnerStrategyNegRevLex(I).sorted_monomials() + [1, + a, + b, a*b, + c, a*c, b*c, a*b*c, + d, a*d, b*d, a*b*d, c*d, a*c*d, b*c*d, a*b*c*d] + sage: GroebnerStrategyDegLex(I).sorted_monomials() + [1, + a, b, c, d, + a*b, a*c, a*d, b*c, b*d, c*d, + a*b*c, a*b*d, a*c*d, b*c*d, + a*b*c*d] + sage: GroebnerStrategyDegRevLex(I).sorted_monomials() + [1, + d, c, b, a, + c*d, b*d, a*d, b*c, a*c, a*b, + b*c*d, a*c*d, a*b*d, a*b*c, + a*b*c*d] + """ + cdef FrozenBitset X + cdef Integer i + cdef list D = [self.bitset_to_int(X) for X in self.E._indices] + D.sort() + if as_dict: + return {i: build_monomial(self.E, self.int_to_bitset(i)) for i in D} + return [build_monomial(self.E, self.int_to_bitset(i)) for i in D] + + def monomial_to_int(self): + """ + Helper function to display the monomials in their term order + used by ``self``. + + EXAMPLES:: + + sage: from sage.algebras.exterior_algebra_groebner import * + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([a, b]) + sage: GroebnerStrategyDegLex(I).sorted_monomials() + {1: 0, + a: 1, b: 2, c: 3, d: 4, + a*b: 5, a*c: 6, a*d: 7, b*c: 8, b*d: 9, c*d: 10, + a*b*c: 11, a*b*d: 12, a*c*d: 13, b*c*d: 14, + a*b*c*d: 15} + sage: GroebnerStrategyDegRevLex(I).monomial_to_int() + {1: 0, + a: 4, b: 3, c: 2, d: 1, + a*b: 10, a*c: 9, a*d: 7, b*c: 8, b*d: 6, c*d: 5, + a*b*c: 14, a*b*d: 13, a*c*d: 12, b*c*d: 11, + a*b*c*d: 15} + """ + B = self.E.basis() + return {B[X]: self.bitset_to_int(X) for X in self.E._indices} + +cdef class GroebnerStrategyNegRevLex(GroebnerStrategy): + """ + Gröbner basis strategy implementing negrevlex ordering. + """ cdef inline Integer bitset_to_int(self, FrozenBitset X): """ Convert ``X`` to an :class:`Integer`. @@ -262,7 +359,7 @@ cdef class GroebnerStrategy: """ Convert a nonnegative integer ``n`` to a :class:`FrozenBitset`. """ - cdef Py_ssize_t i + cdef size_t i if mpz_sgn(n.value) == 0: return FrozenBitset() @@ -274,4 +371,102 @@ cdef class GroebnerStrategy: bitset_set_to(ret._bitset, i, mpz_tstbit(n.value, i)) return ret +cdef class GroebnerStrategyDegRevLex(GroebnerStrategy): + """ + Gröbner basis strategy implementing degree revlex ordering. + """ + def __init__(self, I): + """ + Initialize ``self``. + """ + GroebnerStrategy.__init__(self, I) + self.rank = Integer(self.E.ngens()) + + cdef inline Integer bitset_to_int(self, FrozenBitset X): + """ + Convert ``X`` to an :class:`Integer`. + """ + if X.isempty(): + return Integer(0) + + cdef Integer n = self.rank + cdef long i, deg = degree(X) + cdef long r = 1 + cdef Integer t = Integer(0) + + cdef long elt = bitset_first(X._bitset) + while elt >= 0: + t += Integer(elt).binomial(r) + r += 1 + elt = bitset_next(X._bitset, elt + 1) + return Integer(sum(n.binomial(i) for i in range(deg+1)) - t - 1) + + cdef inline FrozenBitset int_to_bitset(self, Integer n): + """ + Convert a nonnegative integer ``n`` to a :class:`FrozenBitset`. + """ + cdef size_t i + + if mpz_sgn(n.value) == 0: + return FrozenBitset() + + cdef Py_ssize_t deg = 0 + cdef Integer binom = Integer(1) + while n >= binom: + n -= binom + deg += 1 + binom = self.rank.binomial(deg) + + # TODO: Cythonize the from_rank + from sage.combinat.combination import from_rank + return FrozenBitset([self.rank - val - 1 for val in from_rank(n, self.rank, deg)]) + +cdef class GroebnerStrategyDegLex(GroebnerStrategy): + """ + Gröbner basis strategy implementing degree lex ordering. + """ + def __init__(self, I): + """ + Initialize ``self``. + """ + GroebnerStrategy.__init__(self, I) + self.rank = Integer(self.E.ngens()) + + cdef inline Integer bitset_to_int(self, FrozenBitset X): + """ + Convert ``X`` to an :class:`Integer`. + """ + if X.isempty(): + return Integer(0) + + cdef Integer n = self.rank + cdef long i, deg = degree(X) + cdef long r = deg + cdef Integer t = Integer(0) + + cdef long elt = bitset_first(X._bitset) + while elt >= 0: + t += Integer(n - 1 - elt).binomial(r) + r -= 1 + elt = bitset_next(X._bitset, elt + 1) + return Integer(sum(n.binomial(i) for i in range(deg+1)) - t - 1) + + cdef inline FrozenBitset int_to_bitset(self, Integer n): + """ + Convert a nonnegative integer ``n`` to a :class:`FrozenBitset`. + """ + cdef size_t i + + if mpz_sgn(n.value) == 0: + return FrozenBitset() + + cdef Py_ssize_t deg = 0 + cdef Integer binom = Integer(1) + while n >= binom: + n -= binom + deg += 1 + binom = self.rank.binomial(deg) + + from sage.combinat.combination import from_rank + return FrozenBitset(from_rank(n, self.rank, deg)) From f2a25e30e9d71a855d3321d1498825848331e353 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Thu, 14 Jul 2022 12:17:28 +0900 Subject: [PATCH 11/27] Fixing some things and tests. --- src/sage/algebras/clifford_algebra.py | 145 +++--------------- .../algebras/exterior_algebra_groebner.pxd | 2 +- .../algebras/exterior_algebra_groebner.pyx | 19 +-- 3 files changed, 30 insertions(+), 136 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index e4f211c5512..00f15bdbe12 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -1,4 +1,3 @@ -# -*- coding: utf-8 -*- r""" Clifford Algebras @@ -1850,118 +1849,6 @@ def lifted_form(x, y): codomain=self.base_ring(), name="Bilinear Form") - def exterior_bitset_f4(self, I): - r""" - Return a Groebner basis for an ideal `I` of the exterior algebra. - - EXAMPLES:: - - sage: 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: exterior_bitset_f4(I) - (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) - - The example above was computed first using M2: - - E = QQ[a..e, SkewCommutative => true] - I = ideal( 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) - groebnerBasis(I) - - returns: - o3 = | bcd-bce+bde-cde acd-ace+ade-cde abd-abe+ade-bde abc-abe+ace-bce | - """ - - # following https://pi.math.cornell.edu/~djp282/documents/math6140-2017a.pdf - - def get_pairs(pairs): - # this implements the `select` function of Peifer - - # change pairs in here so I don't have to do it after calling get_pairs. - return pairs.pop() # temporarily do Buchbergers - - def f4_sel(P): - # this is the function on page 13 of the original F4 paper. - - return P - - # P is a list of pairs - # TODO: implement the better choice. - - def symbolic_preprocessing(P, G): - # the first/second terms of the S polynomial might - # have to be adjusted a la Stokes' paper - left = set() # the first term of the S polynomial - right = set() # the second term of the S polynomial - L = left.union(right) - done = set(f.monomials()[0] for f in L) - - from itertools import chain - mon_L = set(chain.from_iterable(f.monomials() for f in L)) - - while done != mon_L: - m = sorted(mon_L.difference(done), key = lambda x: (-len(x.support_of_term()), tuple(x.support_of_term())))[0] - done = done.add(m) - for g in G: - if divides(g.monomials()[0], m): - L.add(m * g/g.monomials()[0]) - break - return L - - def f4_reduce(P, G): - # given a current set of pairs P and a - # current basis G, return G' of new basis - - L = symbolic_preprocessing(P, G) - lm_L = set(f.monomials()[0] for f in L) - - d = dict() - for i, f in enumerate(F): - d.update(((i,hash(m)),c) for m,c in f._monomial_coefficients.items()) - M = MatrixSpace(E.base_ring(), len(L), 2^self.ngens(), sparse=True) - poly_matrix = M(d).rref() - - Lprime = set(E._from_dict(dict((FrozenBitset(format(k,'b')[::-1]),v) for k,v in row.dict().items())) for row in poly_matrix) - - Gprime = set() - - for f in Lprime: - if f.monomials()[0] in lm_L: - continue - Gprime.add(f) - - return Gprime - - F = I.gens() - G = set(F) - k = I.ngens() - - from itertools import combinations - pairs = set(combinations(range(k), 2)) # this is Peifer's P - - while pairs: - P = f4_sel(pairs) # this is different from Buchbergers which would be pairs.pop() - Gtemp = f4_reduce(P, G) - pairs.difference_update(P) - - for h in Gtemp: - G.add(h) - k += 1 - pairs.update((i,k) for i in range(k)) - - return G - def _ideal_class_(self, n=0): """ Return the class that is used to implement ideals of ``self``. @@ -2950,7 +2837,7 @@ def reduce(self, f): """ return f.reduce(self) - def groebner_basis(self, term_order="negrevlex"): + def groebner_basis(self, term_order="neglex"): r""" Return the reduced Gröbner basis of ``self``. @@ -2959,13 +2846,13 @@ def groebner_basis(self, term_order="negrevlex"): - ``term_order`` -- the term order used to compute the Gröbner basis; must be one of the following: - * ``"negrevlex"`` -- (default) negative reverse lex order + * ``"neglex"`` -- (default) negative (read right-to-left) lex order * ``"degrevlex"`` -- degree reverse lex order * ``"deglex"`` -- degree lex order EXAMPLES: - We compute an example that was checked against Macaulay2:: + We compute an example:: sage: E. = ExteriorAlgebra(QQ) sage: rels = [c*d*e - b*d*e + b*c*e - b*c*d, @@ -2983,19 +2870,33 @@ def groebner_basis(self, term_order="negrevlex"): With different term orders:: sage: I.groebner_basis("degrevlex") - (-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) + (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) sage: I.groebner_basis("deglex") (-a*c*d + a*c*e - a*d*e + c*d*e, -a*b*c + a*b*d - a*c*d + b*c*d, -a*b*d + a*b*e - a*d*e + b*d*e, -a*b*c + a*b*e - a*c*e + b*c*e) + + The example above was computed first using M2, which agrees with + the ``"degrevlex"`` ordering:: + + E = QQ[a..e, SkewCommutative => true] + I = ideal( 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) + groebnerBasis(I) + + returns: + o3 = | bcd-bce+bde-cde acd-ace+ade-cde abd-abe+ade-bde abc-abe+ace-bce | """ - if term_order == "negrevlex": - from sage.algebras.exterior_algebra_groebner import GroebnerStrategyNegRevLex as strategy + if term_order == "neglex": + from sage.algebras.exterior_algebra_groebner import GroebnerStrategyNegLex as strategy elif term_order == "degrevlex": from sage.algebras.exterior_algebra_groebner import GroebnerStrategyDegRevLex as strategy elif term_order == "deglex": diff --git a/src/sage/algebras/exterior_algebra_groebner.pxd b/src/sage/algebras/exterior_algebra_groebner.pxd index 5d9430a1cce..664d70202fe 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pxd +++ b/src/sage/algebras/exterior_algebra_groebner.pxd @@ -31,7 +31,7 @@ cdef class GroebnerStrategy: cdef Integer bitset_to_int(self, FrozenBitset X) cdef FrozenBitset int_to_bitset(self, Integer n) -cdef class GroebnerStrategyNegRevLex(GroebnerStrategy): +cdef class GroebnerStrategyNegLex(GroebnerStrategy): pass cdef class GroebnerStrategyDegRevLex(GroebnerStrategy): diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index 29998694656..be35add50f7 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -118,7 +118,6 @@ cdef class GroebnerStrategy: """ Perform the preprocessing step. """ - #print("Start preprocessing:", P) cdef CliffordAlgebraElement f, g, f0, f1 cdef set L = set() @@ -159,7 +158,6 @@ cdef class GroebnerStrategy: monL.update(set(f._monomial_coefficients) - done) L.add(f) break - #print("preprocessing:", L) return L cdef inline list reduction(self, list P, list G): @@ -205,10 +203,8 @@ cdef class GroebnerStrategy: P[deg] = [(f0, f1)] while P: - #print("Cur G:", G) Pp = P.pop(min(P)) # The selection: lowest lcm degree Gp = self.reduction(Pp, G) - #print("Reduction yielded:", Gp) G.extend(Gp) for j in range(n, len(G)): f1 = G[j] @@ -223,8 +219,6 @@ cdef class GroebnerStrategy: 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 @@ -254,13 +248,12 @@ cdef class GroebnerStrategy: 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]) + return tuple([~f[self.leading_supp(f)] * f for f in G if f]) cdef Integer bitset_to_int(self, FrozenBitset X): raise NotImplementedError @@ -278,7 +271,7 @@ cdef class GroebnerStrategy: sage: from sage.algebras.exterior_algebra_groebner import * sage: E. = ExteriorAlgebra(QQ) sage: I = E.ideal([x, y]) - sage: GroebnerStrategyNegRevLex(I).sorted_monomials() + sage: GroebnerStrategyNegLex(I).sorted_monomials() [1, x, y, x*y, z, x*z, y*z, x*y*z] sage: GroebnerStrategyDegLex(I).sorted_monomials() [1, x, y, z, x*y, x*z, y*z, x*y*z] @@ -287,7 +280,7 @@ cdef class GroebnerStrategy: sage: E. = ExteriorAlgebra(QQ) sage: I = E.ideal([a, b]) - sage: GroebnerStrategyNegRevLex(I).sorted_monomials() + sage: GroebnerStrategyNegLex(I).sorted_monomials() [1, a, b, a*b, @@ -324,7 +317,7 @@ cdef class GroebnerStrategy: sage: from sage.algebras.exterior_algebra_groebner import * sage: E. = ExteriorAlgebra(QQ) sage: I = E.ideal([a, b]) - sage: GroebnerStrategyDegLex(I).sorted_monomials() + sage: GroebnerStrategyDegLex(I).monomial_to_int() {1: 0, a: 1, b: 2, c: 3, d: 4, a*b: 5, a*c: 6, a*d: 7, b*c: 8, b*d: 9, c*d: 10, @@ -340,9 +333,9 @@ cdef class GroebnerStrategy: B = self.E.basis() return {B[X]: self.bitset_to_int(X) for X in self.E._indices} -cdef class GroebnerStrategyNegRevLex(GroebnerStrategy): +cdef class GroebnerStrategyNegLex(GroebnerStrategy): """ - Gröbner basis strategy implementing negrevlex ordering. + Gröbner basis strategy implementing neglex ordering. """ cdef inline Integer bitset_to_int(self, FrozenBitset X): """ From dba5089f1517642d8654f31d415e42feb6f00f56 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Thu, 14 Jul 2022 13:15:36 +0900 Subject: [PATCH 12/27] Bringing more to the GB strategy classes. --- src/sage/algebras/clifford_algebra.py | 30 +++++++-- .../algebras/clifford_algebra_element.pyx | 11 ++-- .../algebras/exterior_algebra_groebner.pxd | 4 ++ .../algebras/exterior_algebra_groebner.pyx | 64 +++++++++++++------ 4 files changed, 76 insertions(+), 33 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 00f15bdbe12..78d3e212c2e 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -2824,7 +2824,6 @@ def __init__(self, ring, gens, coerce=True, side="twosided"): """ Initialize ``self``. """ - self._groebner_basis = None self._groebner_strategy = None self._homogeneous = all(x.is_super_homogeneous() for x in gens) if self._homogeneous: @@ -2835,7 +2834,24 @@ def reduce(self, f): """ Reduce ``f`` modulo ``self``. """ - return f.reduce(self) + if self._groebner_strategy is None: + self.groebner_basis() + R = self.ring() + return self._groebner_strategy.reduce(R(f)) + + def _contains_(self, f): + r""" + Return ``True`` if ``f`` is in this ideal, + ``False`` otherwise. + + EXAMPLES:: + + .. NOTE:: + + Requires computation of a Groebner basis, which can be a very + expensive operation. + """ + return self.reduce(f).is_zero() def groebner_basis(self, term_order="neglex"): r""" @@ -2903,9 +2919,9 @@ def groebner_basis(self, term_order="neglex"): from sage.algebras.exterior_algebra_groebner import GroebnerStrategyDegLex as strategy else: raise ValueError("invalid term order") - if strategy == self._groebner_strategy: - return self._groebner_basis - self._groebner_strategy = strategy - self._groebner_basis = self._groebner_strategy(self).compute_groebner() - return self._groebner_basis + if strategy == type(self._groebner_strategy): + return self._groebner_strategy.groebner_basis + self._groebner_strategy = strategy(self) + self._groebner_strategy.compute_groebner() + return self._groebner_strategy.groebner_basis diff --git a/src/sage/algebras/clifford_algebra_element.pyx b/src/sage/algebras/clifford_algebra_element.pyx index 93568ebd1d1..0c35b57d983 100644 --- a/src/sage/algebras/clifford_algebra_element.pyx +++ b/src/sage/algebras/clifford_algebra_element.pyx @@ -390,18 +390,19 @@ cdef class ExteriorAlgebraElement(CliffordAlgebraElement): INPUT: - ``I`` -- a list of exterior algebra elements or an ideal - - ``side`` -- the side, ignored if ``I`` is an ideal + - ``left`` -- boolean; if reduce as a left ideal (``True``) + or right ideal (``False``), ignored if ``I`` is an ideal """ from sage.algebras.clifford_algebra import ExteriorAlgebraIdeal if isinstance(I, ExteriorAlgebraIdeal): - left = (I.side() == "left") - I = I.groebner_basis() + return I.reduce(self) f = self E = self._parent - from sage.algebras.exterior_algebra_cython import leading_support + + cdef FrozenBitset lm, s for g in I: - lm = leading_support(g) + lm = g.leading_support() reduction = True while reduction: supp = f.support() diff --git a/src/sage/algebras/exterior_algebra_groebner.pxd b/src/sage/algebras/exterior_algebra_groebner.pxd index 664d70202fe..4e235dae6ef 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pxd +++ b/src/sage/algebras/exterior_algebra_groebner.pxd @@ -17,6 +17,7 @@ cdef class GroebnerStrategy: cdef int side cdef MonoidElement ideal cdef bint homogeneous + cdef public tuple groebner_basis cdef inline bint build_S_poly(self, CliffordAlgebraElement f, CliffordAlgebraElement g) @@ -26,6 +27,9 @@ cdef class GroebnerStrategy: cdef set preprocessing(self, list P, list G) cdef list reduction(self, list P, list G) + cpdef CliffordAlgebraElement reduce(self, CliffordAlgebraElement f) + cdef CliffordAlgebraElement reduce_single(self, CliffordAlgebraElement f, CliffordAlgebraElement g) + # These are the methods that determine the ordering of the monomials. # These must be implemented in subclasses. Declare them as "inline" there. cdef Integer bitset_to_int(self, FrozenBitset X) diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index be35add50f7..cb2dc49c8bf 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -51,6 +51,7 @@ cdef class GroebnerStrategy: Initialize ``self``. """ self.ideal = I + self.groebner_basis = (None,) self.E = I.ring() self.homogeneous = all(x.is_super_homogeneous() for x in I.gens()) if self.homogeneous or I.side() == "left": @@ -221,31 +222,14 @@ cdef class GroebnerStrategy: # 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 tuple 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 = self.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 = self.E.monomial(s - lm) - if self.side == 0: - gp = mon * g - f = f - f[s] / gp[s] * gp - else: - gp = g * mon - f = f - f[s] / gp[s] * gp - break + # TODO: Make this faster by using the previous technique? + f = self.reduce_single(G[i], G[j]) if G[i] != f: G[i] = f if not f: @@ -253,7 +237,45 @@ cdef class GroebnerStrategy: else: pairs.update((k, i) for k in range(n) if k != i) - return tuple([~f[self.leading_supp(f)] * f for f in G if f]) + self.groebner_basis = tuple([~f[self.leading_supp(f)] * f for f in G if f]) + + cpdef CliffordAlgebraElement reduce(self, CliffordAlgebraElement f): + """ + Reduce ``f`` modulo the ideal with Gröbner basis ``G``. + """ + for g in self.groebner_basis: + f = self.reduce_single(f, g) + return f + + cdef CliffordAlgebraElement reduce_single(self, CliffordAlgebraElement f, CliffordAlgebraElement g): + """ + Reduce ``f`` by ``g``. + + .. TODO:: + + Optimize this by doing it in-place and changing the underlying dict of ``f``. + """ + cdef FrozenBitset lm, s + cdef tuple supp + + lm = self.leading_supp(g) + did_reduction = True + while did_reduction: + supp = tuple(f._monomial_coefficients) + did_reduction = False + for s in supp: + if lm <= s: + did_reduction = True + mon = self.E.monomial(s - lm) + if self.side == 0: + gp = mon * g + f -= f[s] / gp[s] * gp + else: + gp = g * mon + f -= f[s] / gp[s] * gp + break + return f + cdef Integer bitset_to_int(self, FrozenBitset X): raise NotImplementedError From 2362a631cf4d34b278f378424d375d2fe3e237c9 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Thu, 14 Jul 2022 13:34:41 +0900 Subject: [PATCH 13/27] Implementing containment of ideals, mostly. --- src/sage/algebras/clifford_algebra.py | 92 ++++++++++++++++++++++++++- 1 file changed, 91 insertions(+), 1 deletion(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 78d3e212c2e..73d461cfbad 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -21,6 +21,8 @@ from sage.structure.unique_representation import UniqueRepresentation from sage.structure.parent import Parent from sage.structure.element import Element +from sage.structure.richcmp import (richcmp_method, op_EQ, op_NE, + op_LT, op_GT, op_LE, op_GE, rich_to_bool) from sage.data_structures.bitset import Bitset, FrozenBitset from sage.algebras.clifford_algebra_element import CliffordAlgebraElement, ExteriorAlgebraElement @@ -2816,6 +2818,7 @@ def chain_complex(self, R=None): return ChainComplex(data, degree=1) +@richcmp_method class ExteriorAlgebraIdeal(Ideal_nc): """ An ideal of the exterior algebra. @@ -2846,12 +2849,99 @@ def _contains_(self, f): EXAMPLES:: + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([x, x*y*z + 2*x*z + 3*y*z]) + sage: I.groebner_basis() + (x, y*z) + sage: x in I + True + sage: y*z in I + True + sage: x + 3*y*z in I + True + sage: x + 3*y in I + False + .. NOTE:: Requires computation of a Groebner basis, which can be a very expensive operation. """ - return self.reduce(f).is_zero() + return not self.reduce(f) + + def __richcmp__(self, other, op): + """ + Compare ``self`` and ``other``. + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([x, x*y*z + 2*x*z + 3*y*z]) + sage: I == I + True + sage: Ip = E.ideal([x, y*z]) + sage: Ip == I + True + sage: Ip <= I + True + sage: Ip < I + False + sage: Ip >= I + True + sage: Ip > I + False + sage: E.ideal([x]) < I + True + sage: E.ideal([x]) <= I + True + sage: I <= E.ideal([x]) + False + """ + if not isinstance(other, ExteriorAlgebraIdeal): + if op == op_EQ: + return False + if op == op_NE: + return True + return NotImplemented + + if self is other: + return rich_to_bool(op, 0) + + # comparison for >= and > : swap the arguments + if op == op_GE: + return other.__richcmp__(self, op_LE) + elif op == op_GT: + return other.__richcmp__(self, op_LT) + + if self.side() == other.side(): + s_gens = set(g for g in self.gens() if g) + o_gens = set(g for g in other.gens() if g) + if set(s_gens) == set(o_gens): + return rich_to_bool(op, 0) + + contained = all(f in other for f in s_gens) + if op == op_LE: + return contained + + contains = all(f in self for f in o_gens) + if op == op_EQ: + return contained and contains + if op == op_NE: + return not (contained and contains) + # remaining case < + return contained and not contains + + + if op in [op_LT, op_LE] and other.side() == "twosided": + if not all(f in other for f in set(self.gens()) if f): + return False + if op == op_LE: + return True + return self.__richcmp__(other, op_NE) + + # Otherwise we will fallback to linear algebra containment + # TODO: Implement this + return NotImplemented def groebner_basis(self, term_order="neglex"): r""" From 6480faaad4f3be343fd55e620e6292b9a177e450 Mon Sep 17 00:00:00 2001 From: "Trevor K. Karn" Date: Tue, 26 Jul 2022 10:28:12 -0500 Subject: [PATCH 14/27] Rebase off of 9.7.beta6 and fix merge conflicts --- src/sage/algebras/clifford_algebra.py | 292 +------------------------- 1 file changed, 1 insertion(+), 291 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 73d461cfbad..900a796881c 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -1628,11 +1628,7 @@ def coproduct_on_basis(self, a): one = self.base_ring().one() L = unshuffle_iterator(tuple(a), one) return self.tensor_square()._from_dict( -<<<<<<< HEAD {tuple(FrozenBitset(e) if e else FrozenBitset('0') for e in t): c for t, c in L if c}, -======= - {tuple(FrozenBitset(e) if e else FrozenBitset() for e in t): c for t,c in L if c}, ->>>>>>> b0f66e328e (Cythonizing the element classes.) coerce=False, remove_zeros=False) @@ -1863,294 +1859,8 @@ def _ideal_class_(self, n=0): """ return ExteriorAlgebraIdeal -<<<<<<< HEAD - class Element(CliffordAlgebraElement): - """ - An element of an exterior algebra. - """ - def _mul_(self, other): - """ - Return ``self`` multiplied by ``other``. - - INPUT: - - - ``other`` -- element of the same exterior algebra as ``self`` - - EXAMPLES:: - - sage: E. = ExteriorAlgebra(QQ) - sage: x*y - x*y - sage: y*x - -x*y - sage: z*y*x - -x*y*z - sage: (x*z)*y - -x*y*z - sage: (3*x + y)^2 - 0 - sage: (x - 3*y + z/3)^2 - 0 - sage: (x+y) * (y+z) - x*y + x*z + y*z - - sage: E. = ExteriorAlgebra(QQ) - sage: (x * y) * (w * z) - -x*y*z*w - sage: x * y * w * z - -x*y*z*w - sage: (z * w) * (x * y) - x*y*z*w - """ - P = self.parent() - zero = P.base_ring().zero() - d = {} - n = P.ngens() - - for ml, cl in self: # ml for "monomial on the left" - for mr, cr in other: # mr for "monomial on the right" - if ml.intersection(mr): - # if they intersect nontrivially, move along. - continue - - if not mr: - t = ml - else: - t = ml.union(mr) - it = iter(mr) - j = next(it) - - num_cross = 0 # keep track of the number of signs - tot_cross = 0 - for i in ml: - num_cross_new = 0 - while i > j: - num_cross_new += 1 - try: - j = next(it) - except StopIteration: - j = n + 1 - tot_cross += num_cross - if tot_cross % 2: - cr = -cr - - d[t] = d.get(t, zero) + cl * cr - if d[t] == zero: - del d[t] - - 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): - left = (I.side() == "left") - I = I.groebner_basis() - - 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 - ``self`` with respect to ``x`` (that is, the element - `\iota_{x}(\text{self})` of the exterior algebra). - - If `V` is an `R`-module, and if `\alpha` is a fixed element of - `V^*`, then the *interior product* with respect to `\alpha` is - an `R`-linear map - `i_{\alpha} \colon \Lambda(V) \to \Lambda(V)`, determined by - the following requirements: - - - `i_{\alpha}(v) = \alpha(v)` for all `v \in V = \Lambda^1(V)`, - - it is a graded derivation of degree `-1`: all `x` and `y` - in `\Lambda(V)` satisfy - - .. MATH:: - - i_{\alpha}(x \wedge y) = (i_{\alpha} x) \wedge y - + (-1)^{\deg x} x \wedge (i_{\alpha} y). - - It can be shown that this map `i_{\alpha}` is graded of - degree `-1` (that is, sends `\Lambda^k(V)` into - `\Lambda^{k-1}(V)` for every `k`). - - When `V` is a finite free `R`-module, the interior product can - also be defined by - - .. MATH:: - - (i_{\alpha} \omega)(u_1, \ldots, u_k) - = \omega(\alpha, u_1, \ldots, u_k), - - where `\omega \in \Lambda^k(V)` is thought of as an - alternating multilinear mapping from - `V^* \times \cdots \times V^*` to `R`. - - Since Sage is only dealing with exterior powers of modules - of the form `R^d` for some nonnegative integer `d`, the - element `\alpha \in V^*` can be thought of as an element of - `V` (by identifying the standard basis of `V = R^d` with its - dual basis). This is how `\alpha` should be passed to this - method. - - We then extend the interior product to all - `\alpha \in \Lambda (V^*)` by - - .. MATH:: - - i_{\beta \wedge \gamma} = i_{\gamma} \circ i_{\beta}. - - INPUT: - - - ``x`` -- element of (or coercing into) `\Lambda^1(V)` - (for example, an element of `V`); this plays the role of - `\alpha` in the above definition - - EXAMPLES:: - - sage: E. = ExteriorAlgebra(QQ) - sage: x.interior_product(x) - 1 - sage: (x + x*y).interior_product(2*y) - -2*x - sage: (x*z + x*y*z).interior_product(2*y - x) - -2*x*z - y*z - z - sage: x.interior_product(E.one()) - x - sage: E.one().interior_product(x) - 0 - sage: x.interior_product(E.zero()) - 0 - sage: E.zero().interior_product(x) - 0 - - REFERENCES: - - - :wikipedia:`Exterior_algebra#Interior_product` - """ - P = self.parent() - return P.sum([c * cx * P.interior_product_on_basis(m, mx) - for m, c in self for mx, cx in x]) - - antiderivation = interior_product - - def hodge_dual(self): - r""" - Return the Hodge dual of ``self``. - - The Hodge dual of an element `\alpha` of the exterior algebra is - defined as `i_{\alpha} \sigma`, where `\sigma` is the volume - form - (:meth:`~sage.algebras.clifford_algebra.ExteriorAlgebra.volume_form`) - and `i_{\alpha}` denotes the antiderivation function with - respect to `\alpha` (see :meth:`interior_product` for the - definition of this). - - .. NOTE:: - - The Hodge dual of the Hodge dual of a homogeneous element - `p` of `\Lambda(V)` equals `(-1)^{k(n-k)} p`, where - `n = \dim V` and `k = \deg(p) = |p|`. - - EXAMPLES:: - - sage: E. = ExteriorAlgebra(QQ) - sage: x.hodge_dual() - y*z - sage: (x*z).hodge_dual() - -y - sage: (x*y*z).hodge_dual() - 1 - sage: [a.hodge_dual().hodge_dual() for a in E.basis()] - [1, x, y, z, x*y, x*z, y*z, x*y*z] - sage: (x + x*y).hodge_dual() - y*z + z - sage: (x*z + x*y*z).hodge_dual() - -y + 1 - sage: E = ExteriorAlgebra(QQ, 'wxyz') - sage: [a.hodge_dual().hodge_dual() for a in E.basis()] - [1, -w, -x, -y, -z, w*x, w*y, w*z, x*y, x*z, y*z, - -w*x*y, -w*x*z, -w*y*z, -x*y*z, w*x*y*z] - """ - volume_form = self.parent().volume_form() - return volume_form.interior_product(self) - - def constant_coefficient(self): - """ - Return the constant coefficient of ``self``. - - .. TODO:: - - Define a similar method for general Clifford algebras once - the morphism to exterior algebras is implemented. - - EXAMPLES:: - - sage: E. = ExteriorAlgebra(QQ) - sage: elt = 5*x + y + x*z + 10 - sage: elt.constant_coefficient() - 10 - sage: x.constant_coefficient() - 0 - """ - return self._monomial_coefficients.get(self.parent().one_basis(), - self.base_ring().zero()) - - def scalar(self, other): - r""" - Return the standard scalar product of ``self`` with ``other``. - - The standard scalar product of `x, y \in \Lambda(V)` is - defined by `\langle x, y \rangle = \langle x^t y \rangle`, where - `\langle a \rangle` denotes the degree-0 term of `a`, and where - `x^t` denotes the transpose - (:meth:`~sage.algebras.clifford_algebra.CliffordAlgebraElement.transpose`) - of `x`. - - .. TODO:: - - Define a similar method for general Clifford algebras once - the morphism to exterior algebras is implemented. - - EXAMPLES:: - - sage: E. = ExteriorAlgebra(QQ) - sage: elt = 5*x + y + x*z - sage: elt.scalar(z + 2*x) - 0 - sage: elt.transpose() * (z + 2*x) - -2*x*y + 5*x*z + y*z - """ - return (self.transpose() * other).constant_coefficient() -======= Element = ExteriorAlgebraElement ->>>>>>> b0f66e328e (Cythonizing the element classes.) + ##################################################################### # Differentials From e95c5466ea763c1cf969cce5efdc069497ba5e97 Mon Sep 17 00:00:00 2001 From: "Trevor K. Karn" Date: Thu, 28 Jul 2022 12:24:21 -0500 Subject: [PATCH 15/27] Add examples --- src/sage/algebras/clifford_algebra.py | 32 +++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 900a796881c..4830c0969cb 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -1851,6 +1851,12 @@ def _ideal_class_(self, n=0): """ Return the class that is used to implement ideals of ``self``. + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: type(E.ideal(x*y - z)) + + TESTS:: sage: E. = ExteriorAlgebra(QQ) @@ -2532,6 +2538,21 @@ def chain_complex(self, R=None): class ExteriorAlgebraIdeal(Ideal_nc): """ An ideal of the exterior algebra. + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal(x*y); I + Twosided Ideal (x*y) of The exterior algebra of rank 3 over Rational Field + + We can also use it to build a quotient:: + + sage: Q = E.quotient(I); Q + Quotient of The exterior algebra of rank 3 over Rational Field by the ideal (x*y) + sage: Q.inject_variables() + Defining xbar, ybar, zbar + sage: xbar * ybar + 0 """ def __init__(self, ring, gens, coerce=True, side="twosided"): """ @@ -2546,6 +2567,17 @@ def __init__(self, ring, gens, coerce=True, side="twosided"): def reduce(self, f): """ Reduce ``f`` modulo ``self``. + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal(x*y); + sage: I.reduce(x*y + x*y*z + z) + z + sage: I.reduce(x*y + x + y) + x + y + sage: I.reduce(x*y + x*y*z) + 0 """ if self._groebner_strategy is None: self.groebner_basis() From b3ef1f706a4acd9687099126b0d7e0a3d6970188 Mon Sep 17 00:00:00 2001 From: "Trevor K. Karn" Date: Thu, 28 Jul 2022 14:43:06 -0500 Subject: [PATCH 16/27] Fix matrix indexing --- src/sage/algebras/exterior_algebra_groebner.pyx | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index cb2dc49c8bf..d76f797bbe8 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -168,19 +168,30 @@ cdef class GroebnerStrategy: cdef set L = self.preprocessing(P, G) cdef Py_ssize_t i from sage.matrix.constructor import matrix - M = matrix({(i, self.bitset_to_int( m)): c + r = 2 ** self.E.ngens() - 1 # r for "rank" or "reverso" + M = matrix({(i, r - self.bitset_to_int( 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(self.leading_supp( f) for f in L) - return [self.E.element_class(self.E, {self.int_to_bitset(Integer(j)): c for j,c in M[i].iteritems()}) + return [self.E.element_class(self.E, {self.int_to_bitset(r - Integer(j)): c for j,c in M[i].iteritems()}) for i,p in enumerate(M.pivots()) - if self.int_to_bitset(Integer(p)) not in lead_supports] + if self.int_to_bitset(r - Integer(p)) not in lead_supports] def compute_groebner(self): """ Compute the reduced ``side`` Gröbner basis for the ideal ``I``. + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([x*y - x, x*y -1]) + sage: I.groebner_basis() + (1,) + sage: J = E.ideal([x*y - x, 2*x*y - 2]) + sage: J.groebner_basis() + (1,) """ cdef FrozenBitset p0, p1 cdef long deg From 692f2f54b95bac1d9694ae571e3f3d9bf8102600 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Tue, 2 Aug 2022 15:40:13 +0900 Subject: [PATCH 17/27] Fixing containment and reduction of exterior algebra ideals. --- src/sage/algebras/clifford_algebra.py | 107 +++++++++++++----- .../algebras/clifford_algebra_element.pyx | 2 +- .../algebras/exterior_algebra_groebner.pxd | 5 +- .../algebras/exterior_algebra_groebner.pyx | 59 +++++----- 4 files changed, 110 insertions(+), 63 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index f9f2e929fe3..dda85132dc7 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -80,15 +80,13 @@ def __init__(self, Qdim): 7 sage: idx._cardinality 128 - sage: idx._maximal_set - 1111111 sage: i = idx.an_element(); i 0 sage: type(i) """ self._nbits = Qdim - self._cardinality = 2**Qdim + self._cardinality = 2 ** Qdim # the if statement here is in case Qdim is 0. category = FiniteEnumeratedSets().Facade() Parent.__init__(self, category=category, facade=True) @@ -2599,9 +2597,24 @@ class ExteriorAlgebraIdeal(Ideal_nc): def __init__(self, ring, gens, coerce=True, side="twosided"): """ Initialize ``self``. + + EXAMPLES: + + We skip the category test because the ideals are not a proper + element class of the monoid of all ideals:: + + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([x*y - x, x*y - 1]) + sage: TestSuite(I).run(skip="_test_category") + + sage: I = E.ideal([x*y - 3, 0, 2*3]) + sage: TestSuite(I).run(skip="_test_category") + + sage: I = E.ideal([]) + sage: TestSuite(I).run(skip="_test_category") """ self._groebner_strategy = None - self._homogeneous = all(x.is_super_homogeneous() for x in gens) + self._homogeneous = all(x.is_super_homogeneous() for x in gens if x) if self._homogeneous: side = "twosided" Ideal_nc.__init__(self, ring, gens, coerce, side) @@ -2620,6 +2633,11 @@ def reduce(self, f): x + y sage: I.reduce(x*y + x*y*z) 0 + + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([a+b*c]) + sage: I.reduce(I.gen(0) * d) + 0 """ if self._groebner_strategy is None: self.groebner_basis() @@ -2680,6 +2698,31 @@ def __richcmp__(self, other, op): True sage: I <= E.ideal([x]) False + + sage: E. = ExteriorAlgebra(QQ) + sage: p = a + b*c + sage: IT = E.ideal([p], side="twosided") + sage: IR = E.ideal([p], side="right") + sage: IL = E.ideal([p], side="left") + sage: IR == IL + False + sage: IR <= IL + False + sage: IR >= IL + False + sage: IL.reduce(p * d) + 2*a*d + sage: IR.reduce(d * p) + -2*a*d + + sage: IR <= IT + True + sage: IL <= IT + True + sage: IT <= IL + False + sage: IT <= IR + False """ if not isinstance(other, ExteriorAlgebraIdeal): if op == op_EQ: @@ -2697,35 +2740,41 @@ def __richcmp__(self, other, op): elif op == op_GT: return other.__richcmp__(self, op_LT) - if self.side() == other.side(): - s_gens = set(g for g in self.gens() if g) - o_gens = set(g for g in other.gens() if g) - if set(s_gens) == set(o_gens): - return rich_to_bool(op, 0) - - contained = all(f in other for f in s_gens) - if op == op_LE: - return contained + s_gens = set(g for g in self.gens() if g) + o_gens = set(g for g in other.gens() if g) - contains = all(f in self for f in o_gens) - if op == op_EQ: - return contained and contains - if op == op_NE: - return not (contained and contains) - # remaining case < - return contained and not contains + if self.side() != other.side(): + if other.side() == "right": + X = set(t * f for t in self.ring().basis() for f in s_gens) + s_gens.update(X) + elif other.side() == "left": + X = set(f * t for t in self.ring().basis() for f in s_gens) + s_gens.update(X) + if set(s_gens) == set(o_gens): + return rich_to_bool(op, 0) - if op in [op_LT, op_LE] and other.side() == "twosided": - if not all(f in other for f in set(self.gens()) if f): - return False - if op == op_LE: - return True - return self.__richcmp__(other, op_NE) + contained = all(f in other for f in s_gens) + if op == op_LE: + return contained + if op == op_NE and not contained: + return True - # Otherwise we will fallback to linear algebra containment - # TODO: Implement this - return NotImplemented + if self.side() != other.side(): + if self.side() == "right": + X = set(t * f for t in self.ring().basis() for f in o_gens) + s_gens.update(X) + elif self.side() == "left": + X = set(f * t for t in self.ring().basis() for f in o_gens) + s_gens.update(X) + + contains = all(f in self for f in o_gens) + if op == op_EQ: + return contained and contains + if op == op_NE: + return not (contained and contains) + # remaining case < + return contained and not contains def groebner_basis(self, term_order="neglex"): r""" diff --git a/src/sage/algebras/clifford_algebra_element.pyx b/src/sage/algebras/clifford_algebra_element.pyx index 0c35b57d983..ba2d89b880d 100644 --- a/src/sage/algebras/clifford_algebra_element.pyx +++ b/src/sage/algebras/clifford_algebra_element.pyx @@ -8,7 +8,7 @@ AUTHORS: """ #***************************************************************************** -# Copyright (C) 2022 Trevor Karn +# Copyright (C) 2022 Trevor K. Karn # (C) 2022 Travis Scrimshaw # # This program is free software: you can redistribute it and/or modify diff --git a/src/sage/algebras/exterior_algebra_groebner.pxd b/src/sage/algebras/exterior_algebra_groebner.pxd index 4e235dae6ef..37ac62d08cc 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pxd +++ b/src/sage/algebras/exterior_algebra_groebner.pxd @@ -17,6 +17,7 @@ cdef class GroebnerStrategy: cdef int side cdef MonoidElement ideal cdef bint homogeneous + cdef Integer rank cdef public tuple groebner_basis cdef inline bint build_S_poly(self, CliffordAlgebraElement f, CliffordAlgebraElement g) @@ -39,8 +40,8 @@ cdef class GroebnerStrategyNegLex(GroebnerStrategy): pass cdef class GroebnerStrategyDegRevLex(GroebnerStrategy): - cdef Integer rank + pass cdef class GroebnerStrategyDegLex(GroebnerStrategy): - cdef Integer rank + pass diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index d76f797bbe8..771af149931 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -6,11 +6,11 @@ of exterior algebra. AUTHORS: -- Trevor Karn, Travis Scrimshaw (July 2022): Initial implementation +- Trevor K. Karn, Travis Scrimshaw (July 2022): Initial implementation """ #***************************************************************************** -# Copyright (C) 2022 Trevor Karn +# Copyright (C) 2022 Trevor K. Karn # (C) 2022 Travis Scrimshaw # # This program is free software: you can redistribute it and/or modify @@ -53,7 +53,8 @@ cdef class GroebnerStrategy: self.ideal = I self.groebner_basis = (None,) self.E = I.ring() - self.homogeneous = all(x.is_super_homogeneous() for x in I.gens()) + self.homogeneous = I._homogeneous + self.rank = Integer(self.E.ngens()) if self.homogeneous or I.side() == "left": self.side = 0 elif I.side() == "right": @@ -120,25 +121,24 @@ cdef class GroebnerStrategy: Perform the preprocessing step. """ cdef CliffordAlgebraElement f, g, f0, f1 + cdef set additions cdef set L = set() - if self.side == 0: + if self.side == 1: for f0, f1 in P: if self.build_S_poly(f0, f1): - L.add(self.partial_S_poly_left(f0, f1)) - L.add(self.partial_S_poly_left(f1, f0)) - elif self.side == 1: - for f0, f1 in P: - if self.build_S_poly(f0, f1): - L.add(self.partial_S_poly_right(f0, f1) for f0,f1 in P) - L.add(self.partial_S_poly_right(f1, f0) for f0,f1 in P) - if self.side == 2: + L.add(self.partial_S_poly_right(f0, f1)) + L.add(self.partial_S_poly_right(f1, f0)) + else: # We compute a left Gröbner basis for two-sided ideals for f0, f1 in P: if self.build_S_poly(f0, f1): L.add(self.partial_S_poly_left(f0, f1)) L.add(self.partial_S_poly_left(f1, f0)) - L.add(self.partial_S_poly_right(f0, f1)) - L.add(self.partial_S_poly_right(f1, f0)) + + if self.side == 2 and not self.homogeneous: + # Add in all S-poly times positive degree monomials + additions = set(f * t for t in self.E.basis() for f in L) + L.update(f for f in additions if f) cdef set done = set(self.leading_supp(f) for f in L) cdef set monL = set() @@ -168,7 +168,7 @@ cdef class GroebnerStrategy: cdef set L = self.preprocessing(P, G) cdef Py_ssize_t i from sage.matrix.constructor import matrix - r = 2 ** self.E.ngens() - 1 # r for "rank" or "reverso" + cdef Integer r = Integer(2) ** self.rank - Integer(1) # r for "rank" or "reverso" M = matrix({(i, r - self.bitset_to_int( m)): c for i,f in enumerate(L) for m,c in ( f)._monomial_coefficients.items()}, @@ -186,18 +186,29 @@ cdef class GroebnerStrategy: EXAMPLES:: sage: E. = ExteriorAlgebra(QQ) - sage: I = E.ideal([x*y - x, x*y -1]) + sage: I = E.ideal([x*y - x, x*y - 1], side="left") sage: I.groebner_basis() (1,) - sage: J = E.ideal([x*y - x, 2*x*y - 2]) + sage: J = E.ideal([x*y - x, 2*x*y - 2], side="left") sage: J.groebner_basis() (1,) + + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([a+b*c]) + sage: I.groebner_basis() + (b*c + a, a*c, a*b, a*d) """ cdef FrozenBitset p0, p1 cdef long deg cdef Py_ssize_t i, j, k + cdef set additions cdef list G = [f for f in self.ideal.gens() if f] # Remove 0s + if self.side == 2 and not self.homogeneous: + # Add in all S-poly times positive degree monomials + additions = set(f * t for t in self.E.basis() for f in G) + G.extend(f for f in additions if f) + cdef Py_ssize_t n = len(G) cdef dict P = {} cdef list Gp @@ -401,13 +412,6 @@ cdef class GroebnerStrategyDegRevLex(GroebnerStrategy): """ Gröbner basis strategy implementing degree revlex ordering. """ - def __init__(self, I): - """ - Initialize ``self``. - """ - GroebnerStrategy.__init__(self, I) - self.rank = Integer(self.E.ngens()) - cdef inline Integer bitset_to_int(self, FrozenBitset X): """ Convert ``X`` to an :class:`Integer`. @@ -451,13 +455,6 @@ cdef class GroebnerStrategyDegLex(GroebnerStrategy): """ Gröbner basis strategy implementing degree lex ordering. """ - def __init__(self, I): - """ - Initialize ``self``. - """ - GroebnerStrategy.__init__(self, I) - self.rank = Integer(self.E.ngens()) - cdef inline Integer bitset_to_int(self, FrozenBitset X): """ Convert ``X`` to an :class:`Integer`. From 6c92f180c49697fac6d2d85013cbc4da58b405e1 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Tue, 2 Aug 2022 17:06:29 +0900 Subject: [PATCH 18/27] Implementing multiplication of exterior algebra ideals. --- src/sage/algebras/clifford_algebra.py | 72 +++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index dda85132dc7..88046a341a7 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -2776,6 +2776,78 @@ def __richcmp__(self, other, op): # remaining case < return contained and not contains + def __mul__(self, other): + """ + Return the product of ``self`` with ``other``. + + .. WARNING:: + + If ``self`` is a right ideal and ``other`` is a left ideal, + this returns a submodule rather than an ideal. + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + + sage: I = E.ideal([a + 1], side="left") + sage: J = I * I; J + Left Ideal (2*a + 1, a, b, c, d, a*b, a*c, a*d, 2*a*b*c + b*c, 2*a*b*d + b*d, + 2*a*c*d + c*d, a*b*c, a*b*d, a*c*d, b*c*d, a*b*c*d) + of The exterior algebra of rank 4 over Rational Field + sage: J.groebner_basis() + (1,) + sage: I.gen(0)^2 + 2*a + 1 + + sage: J = E.ideal([b+c]) + sage: I * J + Twosided Ideal (a*b + a*c + b + c) of The exterior algebra of rank 4 over Rational Field + sage: J * I + Left Ideal (-a*b - a*c + b + c) of The exterior algebra of rank 4 over Rational Field + + sage: K = J * I + sage: K + Left Ideal (-a*b - a*c + b + c) of The exterior algebra of rank 4 over Rational Field + sage: E.ideal([J.gen(0) * d * I.gen(0)], side="left") <= K + True + + sage: J = E.ideal([b + c*d], side="right") + sage: I * J + Twosided Ideal (a*c*d + a*b + c*d + b) of The exterior algebra of rank 4 over Rational Field + sage: X = J * I; X + Free module generated by {0, 1, 2, 3, 4, 5, 6, 7} over Rational Field + sage: [X.lift(b) for b in X.basis()] + [c*d + b, -a*c*d + a*b, b*c, b*d, a*b*c, a*b*d, b*c*d, a*b*c*d] + sage: p = X.lift(X.basis()[0]) + sage: p + c*d + b + sage: a * p # not a left ideal + a*c*d + a*b + + sage: I = E.ideal([a + 1], side="right") + sage: E.ideal([1]) * I + Twosided Ideal (a + 1) of The exterior algebra of rank 4 over Rational Field + sage: I * E.ideal([1]) + Right Ideal (a + 1) of The exterior algebra of rank 4 over Rational Field + """ + if not isinstance(other, ExteriorAlgebraIdeal) or self.ring() != other.ring(): + return super().__mul__(other) + + if self._homogeneous or other._homogeneous or (self.side() == "left" and other.side() == "right"): + gens = (x * y for x in self.gens() for y in other.gens()) + else: + gens = (x * t * y for t in self.ring().basis() for x in self.gens() for y in other.gens()) + gens = [z for z in gens if z] + + if self.side() == "right" and other.side() == "left": + return self.ring().submodule(gens) + + if self.side() == "left" or self.side() == "twosided": + if other.side() == "right" or other.side() == "twosided": + return self.ring().ideal(gens, side="twosided") + return self.ring().ideal(gens, side="left") + return self.ring().ideal(gens, side="right") + def groebner_basis(self, term_order="neglex"): r""" Return the reduced Gröbner basis of ``self``. From 698e051b3e7a8775a2146839110464ecb322d808 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Tue, 2 Aug 2022 17:22:21 +0900 Subject: [PATCH 19/27] Some last details for pyflakes and full coverage. --- src/sage/algebras/clifford_algebra.py | 13 ++++++- .../algebras/clifford_algebra_element.pyx | 13 +++++++ .../algebras/exterior_algebra_groebner.pyx | 34 +++++++++++++++++-- 3 files changed, 56 insertions(+), 4 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 88046a341a7..0969856f4c6 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -19,7 +19,6 @@ #***************************************************************************** from sage.misc.cachefunc import cached_method -from sage.misc.lazy_attribute import lazy_attribute from sage.structure.unique_representation import UniqueRepresentation from sage.structure.parent import Parent from sage.structure.element import Element @@ -30,6 +29,7 @@ from sage.algebras.clifford_algebra_element import CliffordAlgebraElement, ExteriorAlgebraElement from sage.categories.algebras_with_basis import AlgebrasWithBasis from sage.categories.hopf_algebras_with_basis import HopfAlgebrasWithBasis +from sage.categories.fields import Fields from sage.categories.finite_enumerated_sets import FiniteEnumeratedSets from sage.modules.with_basis.morphism import ModuleMorphismByLinearity from sage.categories.poor_man_map import PoorManMap @@ -2905,7 +2905,18 @@ def groebner_basis(self, term_order="neglex"): returns: o3 = | bcd-bce+bde-cde acd-ace+ade-cde abd-abe+ade-bde abc-abe+ace-bce | + + TESTS:: + + sage: E. = ExteriorAlgebra(ZZ) + sage: I = E.ideal([a+1, b*c+d]) + sage: I.groebner_basis() + Traceback (most recent call last): + ... + NotImplementedError: only implemented over fields """ + if self.ring().base_ring() not in Fields(): + raise NotImplementedError("only implemented over fields") if term_order == "neglex": from sage.algebras.exterior_algebra_groebner import GroebnerStrategyNegLex as strategy elif term_order == "degrevlex": diff --git a/src/sage/algebras/clifford_algebra_element.pyx b/src/sage/algebras/clifford_algebra_element.pyx index ba2d89b880d..ac74f313f76 100644 --- a/src/sage/algebras/clifford_algebra_element.pyx +++ b/src/sage/algebras/clifford_algebra_element.pyx @@ -392,6 +392,19 @@ cdef class ExteriorAlgebraElement(CliffordAlgebraElement): - ``I`` -- a list of exterior algebra elements or an ideal - ``left`` -- boolean; if reduce as a left ideal (``True``) or right ideal (``False``), ignored if ``I`` is an ideal + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: f = (a + b*c) * d + sage: f.reduce([a + b*c], True) + 2*a*d + sage: f.reduce([a + b*c], False) + 0 + + sage: I = E.ideal([a + b*c]) + sage: f.reduce(I) + 0 """ from sage.algebras.clifford_algebra import ExteriorAlgebraIdeal if isinstance(I, ExteriorAlgebraIdeal): diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index 771af149931..3cd7cf3911e 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -49,6 +49,14 @@ cdef class GroebnerStrategy: def __init__(self, I): """ Initialize ``self``. + + EXAMPLES:: + + sage: from sage.algebras.exterior_algebra_groebner import GroebnerStrategy + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([a + 1], side="left") + sage: G = GroebnerStrategy(I) + sage: TestSuite(G).run(skip="_test_pickling") """ self.ideal = I self.groebner_basis = (None,) @@ -187,15 +195,15 @@ cdef class GroebnerStrategy: sage: E. = ExteriorAlgebra(QQ) sage: I = E.ideal([x*y - x, x*y - 1], side="left") - sage: I.groebner_basis() + sage: I.groebner_basis() # indirect doctest (1,) sage: J = E.ideal([x*y - x, 2*x*y - 2], side="left") - sage: J.groebner_basis() + sage: J.groebner_basis() # indirect doctest (1,) sage: E. = ExteriorAlgebra(QQ) sage: I = E.ideal([a+b*c]) - sage: I.groebner_basis() + sage: I.groebner_basis() # indirect doctest (b*c + a, a*c, a*b, a*d) """ cdef FrozenBitset p0, p1 @@ -264,6 +272,26 @@ cdef class GroebnerStrategy: cpdef CliffordAlgebraElement reduce(self, CliffordAlgebraElement f): """ Reduce ``f`` modulo the ideal with Gröbner basis ``G``. + + EXAMPLES:: + + sage: 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.reduce(a*b*e) + a*b*e + sage: I.reduce(b*d*e) + a*b*d - a*b*e + a*d*e + sage: I.reduce(c*d*e) + a*c*d - a*c*e + a*d*e + sage: I.reduce(a*b*c*d*e) + 0 + sage: I.reduce(a*b*c*d) + 0 """ for g in self.groebner_basis: f = self.reduce_single(f, g) From 12d76e7aff60783f4c9099ff1b7319d540926409 Mon Sep 17 00:00:00 2001 From: "Trevor K. Karn" Date: Thu, 4 Aug 2022 19:52:56 -0500 Subject: [PATCH 20/27] Revert change in morphisms --- src/sage/modules/with_basis/morphism.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sage/modules/with_basis/morphism.py b/src/sage/modules/with_basis/morphism.py index 2efcc4f25e9..b8fe98111b2 100644 --- a/src/sage/modules/with_basis/morphism.py +++ b/src/sage/modules/with_basis/morphism.py @@ -404,10 +404,10 @@ def __call__(self, *args): if self._is_module_with_basis_over_same_base_ring: return self.codomain().linear_combination( (self._on_basis(*(before+(index,)+after)), coeff ) - for (index, coeff) in mc.items() if self._on_basis(index)) + for (index, coeff) in mc.items()) else: return sum((coeff * self._on_basis(*(before+(index,)+after)) - for (index, coeff) in mc.items() if self._on_basis(index)), self._zero) + for (index, coeff) in mc.items()), self._zero) # As per the specs of Map, we should in fact implement _call_. # However we currently need to abuse Map.__call__ (which strict From ea9610d76ff08beca1b2916bdadf621cd61581be Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Sat, 6 Aug 2022 15:08:22 +0900 Subject: [PATCH 21/27] Caching the leading support FrozenBitset and Integer. --- .../algebras/exterior_algebra_groebner.pxd | 16 +- .../algebras/exterior_algebra_groebner.pyx | 161 ++++++++++++------ 2 files changed, 122 insertions(+), 55 deletions(-) diff --git a/src/sage/algebras/exterior_algebra_groebner.pxd b/src/sage/algebras/exterior_algebra_groebner.pxd index 37ac62d08cc..7708930faa8 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pxd +++ b/src/sage/algebras/exterior_algebra_groebner.pxd @@ -11,6 +11,11 @@ from sage.structure.element cimport MonoidElement cdef long degree(FrozenBitset X) cdef CliffordAlgebraElement build_monomial(Parent E, FrozenBitset supp) +cdef class GBElement: + cdef CliffordAlgebraElement elt + cdef FrozenBitset ls # leading support as a bitset + cdef Integer lsi # leading support as an Integer + # Grobner basis functions cdef class GroebnerStrategy: cdef Parent E # the exterior algebra @@ -20,11 +25,14 @@ cdef class GroebnerStrategy: cdef Integer rank cdef public tuple groebner_basis - cdef inline bint build_S_poly(self, CliffordAlgebraElement f, CliffordAlgebraElement g) + cdef inline GBElement build_elt(self, CliffordAlgebraElement f) + cdef inline GBElement prod_GB_term(self, GBElement f, FrozenBitset t) + cdef inline GBElement prod_term_GB(self, FrozenBitset t, GBElement f) + cdef inline bint build_S_poly(self, GBElement f, GBElement g) - cdef inline FrozenBitset leading_supp(self, CliffordAlgebraElement f) - cdef inline partial_S_poly_left(self, CliffordAlgebraElement f, CliffordAlgebraElement g) - cdef inline partial_S_poly_right(self, CliffordAlgebraElement f, CliffordAlgebraElement g) + cdef inline FrozenBitset leading_support(self, CliffordAlgebraElement f) + cdef inline partial_S_poly_left(self, GBElement f, GBElement g) + cdef inline partial_S_poly_right(self, GBElement f, GBElement g) cdef set preprocessing(self, list P, list G) cdef list reduction(self, list P, list G) diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index 3cd7cf3911e..5a1858b08e0 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -1,4 +1,5 @@ -""" +# cython: linetrace=True +r""" Exterior algebras Gröbner bases This contains the backend implementations in Cython for the Gröbner bases @@ -24,6 +25,7 @@ from sage.libs.gmp.mpz cimport mpz_sizeinbase, mpz_setbit, mpz_tstbit, mpz_cmp_s from sage.data_structures.bitset_base cimport (bitset_t, bitset_init, bitset_first, bitset_next, bitset_set_to, bitset_len) from sage.structure.parent cimport Parent +from sage.structure.richcmp cimport richcmp cdef inline long degree(FrozenBitset X): """ @@ -38,6 +40,18 @@ cdef inline CliffordAlgebraElement build_monomial(Parent E, FrozenBitset supp): """ return E.element_class(E, {supp: E._base.one()}) +cdef class GBElement: + def __init__(self, CliffordAlgebraElement x, FrozenBitset ls, Integer n): + self.elt = x + self.ls = ls + self.lsi = n + + def __hash__(self): + return int(self.lsi) + + def __richcmp__(self, other, int op): + return richcmp(self.elt, ( other).elt, op) + cdef class GroebnerStrategy: """ A strategy for computing a Gröbner basis. @@ -70,14 +84,14 @@ cdef class GroebnerStrategy: else: self.side = 2 - cdef inline FrozenBitset leading_supp(self, CliffordAlgebraElement f): + cdef inline FrozenBitset leading_support(self, CliffordAlgebraElement f): """ Return the leading support of the exterior algebra element ``f``. """ cdef dict mc = f._monomial_coefficients return self.int_to_bitset(max(self.bitset_to_int(k) for k in mc)) - cdef inline partial_S_poly_left(self, CliffordAlgebraElement f, CliffordAlgebraElement g): + cdef inline partial_S_poly_left(self, GBElement f, GBElement g): """ Compute one half of the `S`-polynomial for ``f`` and ``g``. @@ -87,13 +101,14 @@ cdef class GroebnerStrategy: LCM(LM(f), LM(g)) / LT(f) \cdot f. """ - cdef FrozenBitset lmf = self.leading_supp(f) - cdef FrozenBitset lmg = self.leading_supp(g) - cdef FrozenBitset D = lmg.difference(lmf) - ret = build_monomial(self.E, D) * f - return (~ret[lmf._union(lmg)]) * ret + cdef FrozenBitset D = g.ls.difference(f.ls) + cdef GBElement ret = self.prod_term_GB(D, f) + inv = ~ret.elt[ret.ls] + for k in ret.elt._monomial_coefficients: + ret.elt._monomial_coefficients[k] *= inv + return ret - cdef inline partial_S_poly_right(self, CliffordAlgebraElement f, CliffordAlgebraElement g): + cdef inline partial_S_poly_right(self, GBElement f, GBElement g): """ Compute one half of the `S`-polynomial for ``f`` and ``g``. @@ -103,13 +118,48 @@ cdef class GroebnerStrategy: f \cdot LCM(LM(f), LM(g)) / LT(f). """ - cdef FrozenBitset lmf = self.leading_supp(f) - cdef FrozenBitset lmg = self.leading_supp(g) - cdef FrozenBitset D = lmg.difference(lmf) - ret = f * build_monomial(self.E, D) - return ret * (~ret[lmf._union(lmg)]) + cdef FrozenBitset D = g.ls.difference(f.ls) + cdef GBElement ret = self.prod_GB_term(f, D) + inv = ~ret.elt[ret.ls] + for k in ret.elt._monomial_coefficients: + ret.elt._monomial_coefficients[k] *= inv + return ret + + cdef inline GBElement build_elt(self, CliffordAlgebraElement f): + """ + Convert ``f`` into a ``GBElement``. + """ + cdef dict mc = f._monomial_coefficients + #if not mc: + # return GBElement(f, FrozenBitset(), -1) + cdef Integer r = max(self.bitset_to_int(k) for k in mc) + return GBElement(f, self.int_to_bitset(r), r) + + cdef inline GBElement prod_GB_term(self, GBElement f, FrozenBitset t): + """ + Return the GBElement corresponding to ``f * t``. + + .. WARNING:: + + This assumes the leading support is ``f.ls._union(t)``. + """ + ret = f.elt * build_monomial(self.E, t) + cdef FrozenBitset ls = f.ls._union(t) + return GBElement(ret, ls, self.bitset_to_int(ls)) + + cdef inline GBElement prod_term_GB(self, FrozenBitset t, GBElement f): + """ + Return the GBElement corresponding to ``t * f``. + + .. WARNING:: + + This assumes the leading support is ``f.ls._union(t)``. + """ + ret = build_monomial(self.E, t) * f.elt + cdef FrozenBitset ls = f.ls._union(t) + return GBElement(ret, ls, self.bitset_to_int(ls)) - cdef inline bint build_S_poly(self, CliffordAlgebraElement f, CliffordAlgebraElement g): + cdef inline bint build_S_poly(self, GBElement f, GBElement g): r""" Check to see if we should build the `S`-polynomial. @@ -122,13 +172,13 @@ cdef class GroebnerStrategy: if not self.homogeneous: return True - return ( self.leading_supp(f).intersection(self.leading_supp(g))).isempty() + return ( f.ls.intersection(g.ls)).isempty() cdef inline set preprocessing(self, list P, list G): """ Perform the preprocessing step. """ - cdef CliffordAlgebraElement f, g, f0, f1 + cdef GBElement f, g, f0, f1 cdef set additions cdef set L = set() @@ -145,13 +195,13 @@ cdef class GroebnerStrategy: if self.side == 2 and not self.homogeneous: # Add in all S-poly times positive degree monomials - additions = set(f * t for t in self.E.basis() for f in L) - L.update(f for f in additions if f) + additions = set(( f).elt * t for t in self.E.basis() for f in L) + L.update(self.build_elt(f) for f in additions if f) - cdef set done = set(self.leading_supp(f) for f in L) + cdef set done = set(( f).ls for f in L) cdef set monL = set() for f in L: - monL.update(f._monomial_coefficients) + monL.update(f.elt._monomial_coefficients) monL.difference_update(done) while monL: @@ -159,12 +209,12 @@ cdef class GroebnerStrategy: done.add(m) monL.remove(m) for g in G: - lm = self.leading_supp(g) + lm = g.ls if lm <= m: - f = build_monomial(self.E, m.difference(lm)) * g + f = self.prod_term_GB( m.difference(lm), g) if f in L: break - monL.update(set(f._monomial_coefficients) - done) + monL.update(set(f.elt._monomial_coefficients) - done) L.add(f) break return L @@ -179,13 +229,15 @@ cdef class GroebnerStrategy: cdef Integer r = Integer(2) ** self.rank - Integer(1) # r for "rank" or "reverso" M = matrix({(i, r - self.bitset_to_int( m)): c for i,f in enumerate(L) - for m,c in ( f)._monomial_coefficients.items()}, + for m,c in ( f).elt._monomial_coefficients.items()}, sparse=True) M.echelonize() # Do this in place - lead_supports = set(self.leading_supp( f) for f in L) - return [self.E.element_class(self.E, {self.int_to_bitset(r - Integer(j)): c for j,c in M[i].iteritems()}) - for i,p in enumerate(M.pivots()) - if self.int_to_bitset(r - Integer(p)) not in lead_supports] + lead_supports = set(( f).lsi for f in L) + return [GBElement(self.E.element_class(self.E, {self.int_to_bitset(r - Integer(j)): c for j,c in M[i].iteritems()}), + self.int_to_bitset(Integer(r - p)), + Integer(r - p)) + for i, p in enumerate(M.pivots()) + if r - Integer(p) not in lead_supports] def compute_groebner(self): """ @@ -210,23 +262,24 @@ cdef class GroebnerStrategy: cdef long deg cdef Py_ssize_t i, j, k cdef set additions - cdef list G = [f for f in self.ideal.gens() if f] # Remove 0s + cdef GBElement f0, f1 + cdef list G = [self.build_elt(f) for f in self.ideal.gens() if f] # Remove 0s + cdef list Gp if self.side == 2 and not self.homogeneous: # Add in all S-poly times positive degree monomials - additions = set(f * t for t in self.E.basis() for f in G) - G.extend(f for f in additions if f) + additions = set(( f0).elt * t for t in self.E.basis() for f0 in G) + G.extend(self.build_elt(f) for f in additions if f) cdef Py_ssize_t n = len(G) cdef dict P = {} - cdef list Gp for i in range(n): - f0 = G[i] - p0 = self.leading_supp(f0) + f0 = G[i] + p0 = f0.ls for j in range(i+1, n): - f1 = G[j] - p1 = self.leading_supp(f1) + f1 = G[j] + p1 = f1.ls deg = degree( (p0._union(p1))) if deg in P: P[deg].append((f0, f1)) @@ -239,10 +292,10 @@ cdef class GroebnerStrategy: G.extend(Gp) for j in range(n, len(G)): f1 = G[j] - p1 = self.leading_supp(f1) + p1 = f1.ls for i in range(j): f0 = G[i] - p0 = self.leading_supp(f0) + p0 = f0.ls deg = degree( (p0._union(p1))) if deg in P: P[deg].append((f0, f1)) @@ -250,24 +303,31 @@ cdef class GroebnerStrategy: P[deg] = [(f0, f1)] n = len(G) + G.sort(key=lambda x: ( x).lsi) + + for i in range(len(G)): + G[i] = ( G[i]).elt + # 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 set pairs = set((i, j) for i in range(n) for j in range(i+1,n)) cdef tuple supp cdef bint did_reduction cdef FrozenBitset lm, s while pairs: i,j = pairs.pop() + f = G[i] + g = G[j] # We perform the classical reduction algorithm here on each pair # TODO: Make this faster by using the previous technique? - f = self.reduce_single(G[i], G[j]) - if G[i] != f: - G[i] = f - if not f: + fp = self.reduce_single(f, g) + if f != fp: + G[i] = fp + if not fp: pairs.difference_update((k, i) for k in range(n)) else: pairs.update((k, i) for k in range(n) if k != i) - self.groebner_basis = tuple([~f[self.leading_supp(f)] * f for f in G if f]) + self.groebner_basis = tuple([~f[self.leading_support(f)] * f for f in G if f]) cpdef CliffordAlgebraElement reduce(self, CliffordAlgebraElement f): """ @@ -294,7 +354,7 @@ cdef class GroebnerStrategy: 0 """ for g in self.groebner_basis: - f = self.reduce_single(f, g) + f = self.reduce_single(f, g) return f cdef CliffordAlgebraElement reduce_single(self, CliffordAlgebraElement f, CliffordAlgebraElement g): @@ -305,20 +365,19 @@ cdef class GroebnerStrategy: Optimize this by doing it in-place and changing the underlying dict of ``f``. """ - cdef FrozenBitset lm, s + cdef FrozenBitset lm = self.leading_support(g), s + cdef bint did_reduction = True cdef tuple supp - lm = self.leading_supp(g) - did_reduction = True while did_reduction: supp = tuple(f._monomial_coefficients) did_reduction = False for s in supp: if lm <= s: did_reduction = True - mon = self.E.monomial(s - lm) + mon = build_monomial(self.E, s - lm) if self.side == 0: - gp = mon * g + gp = mon * g.elt f -= f[s] / gp[s] * gp else: gp = g * mon From d21d2e82c724caf14234c964cb108ae0291041dc Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Sat, 6 Aug 2022 15:11:06 +0900 Subject: [PATCH 22/27] Forgot to remove the Cython profiling directive. --- src/sage/algebras/exterior_algebra_groebner.pyx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index 5a1858b08e0..3a2e73de41c 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -1,4 +1,3 @@ -# cython: linetrace=True r""" Exterior algebras Gröbner bases From 57345553ed8801593ab22ea4842afbaa5579b4c8 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Sun, 7 Aug 2022 18:53:34 +0900 Subject: [PATCH 23/27] Fixing a number of small bugs and adding option for non-reduced GB. --- src/sage/algebras/clifford_algebra.py | 72 ++++++++++--- .../algebras/exterior_algebra_groebner.pxd | 1 + .../algebras/exterior_algebra_groebner.pyx | 101 +++++++++++++----- 3 files changed, 133 insertions(+), 41 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 0969856f4c6..55c3957ca3c 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -2614,6 +2614,7 @@ def __init__(self, ring, gens, coerce=True, side="twosided"): sage: TestSuite(I).run(skip="_test_category") """ self._groebner_strategy = None + self._reduced = False self._homogeneous = all(x.is_super_homogeneous() for x in gens if x) if self._homogeneous: side = "twosided" @@ -2652,7 +2653,7 @@ def _contains_(self, f): EXAMPLES:: sage: E. = ExteriorAlgebra(QQ) - sage: I = E.ideal([x, x*y*z + 2*x*z + 3*y*z]) + sage: I = E.ideal([x, x*y*z + 2*x*z + 3*y*z], side="left") sage: I.groebner_basis() (x, y*z) sage: x in I @@ -2663,6 +2664,10 @@ def _contains_(self, f): True sage: x + 3*y in I False + sage: x*y in I + True + sage: x + x*y + y*z + x*z in I + True .. NOTE:: @@ -2848,9 +2853,9 @@ def __mul__(self, other): return self.ring().ideal(gens, side="left") return self.ring().ideal(gens, side="right") - def groebner_basis(self, term_order="neglex"): + def groebner_basis(self, term_order=None, reduced=True): r""" - Return the reduced Gröbner basis of ``self``. + Return the (reduced) Gröbner basis of ``self``. INPUT: @@ -2861,6 +2866,9 @@ def groebner_basis(self, term_order="neglex"): * ``"degrevlex"`` -- degree reverse lex order * ``"deglex"`` -- degree lex order + - ``reduced`` -- (default: ``True``) whether or not to return the + reduced Gröbner basis + EXAMPLES: We compute an example:: @@ -2873,10 +2881,10 @@ def groebner_basis(self, term_order="neglex"): ....: b*c*d - a*c*d + a*b*d - a*b*c] sage: I = E.ideal(rels) sage: I.groebner_basis() - (-a*c*d + a*c*e - a*d*e + c*d*e, - -a*b*c + a*b*d - a*c*d + b*c*d, + (-a*b*c + a*b*d - a*c*d + b*c*d, + -a*b*c + a*b*e - a*c*e + b*c*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) + -a*c*d + a*c*e - a*d*e + c*d*e) With different term orders:: @@ -2887,10 +2895,10 @@ def groebner_basis(self, term_order="neglex"): a*b*c - a*b*e + a*c*e - b*c*e) sage: I.groebner_basis("deglex") - (-a*c*d + a*c*e - a*d*e + c*d*e, - -a*b*c + a*b*d - a*c*d + b*c*d, + (-a*b*c + a*b*d - a*c*d + b*c*d, + -a*b*c + a*b*e - a*c*e + b*c*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) + -a*c*d + a*c*e - a*d*e + c*d*e) The example above was computed first using M2, which agrees with the ``"degrevlex"`` ordering:: @@ -2906,6 +2914,25 @@ def groebner_basis(self, term_order="neglex"): returns: o3 = | bcd-bce+bde-cde acd-ace+ade-cde abd-abe+ade-bde abc-abe+ace-bce | + By default, the Gröbner basis is reduced, but we can get non-reduced + Gröber bases (which are not unique):: + + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([x+y*z]) + sage: I.groebner_basis(reduced=False) + (x*y, x*z, y*z + x, y*z + x, x*y*z) + sage: I.groebner_basis(reduced=True) + (x*y, x*z, y*z + x) + + However, if we have already computed a reduced Gröbner basis (with + a given term order), then we return that:: + + sage: I = E.ideal([x+y*z]) # A fresh ideal + sage: I.groebner_basis() + (x*y, x*z, y*z + x) + sage: I.groebner_basis(reduced=False) + (x*y, x*z, y*z + x) + TESTS:: sage: E. = ExteriorAlgebra(ZZ) @@ -2917,17 +2944,28 @@ def groebner_basis(self, term_order="neglex"): """ if self.ring().base_ring() not in Fields(): raise NotImplementedError("only implemented over fields") - if term_order == "neglex": - from sage.algebras.exterior_algebra_groebner import GroebnerStrategyNegLex as strategy - elif term_order == "degrevlex": - from sage.algebras.exterior_algebra_groebner import GroebnerStrategyDegRevLex as strategy - elif term_order == "deglex": - from sage.algebras.exterior_algebra_groebner import GroebnerStrategyDegLex as strategy + if term_order is None: + if self._groebner_strategy is not None: + strategy = type(self._groebner_strategy) + else: + from sage.algebras.exterior_algebra_groebner import GroebnerStrategyNegLex as strategy else: - raise ValueError("invalid term order") + if term_order == "neglex": + from sage.algebras.exterior_algebra_groebner import GroebnerStrategyNegLex as strategy + elif term_order == "degrevlex": + from sage.algebras.exterior_algebra_groebner import GroebnerStrategyDegRevLex as strategy + elif term_order == "deglex": + from sage.algebras.exterior_algebra_groebner import GroebnerStrategyDegLex as strategy + else: + raise ValueError("invalid term order") if strategy == type(self._groebner_strategy): + if self._reduced or not reduced: + return self._groebner_strategy.groebner_basis + self._reduced = reduced + self._groebner_strategy.reduce_computed_gb() return self._groebner_strategy.groebner_basis self._groebner_strategy = strategy(self) - self._groebner_strategy.compute_groebner() + self._groebner_strategy.compute_groebner(reduced=reduced) + self._reduced = reduced return self._groebner_strategy.groebner_basis diff --git a/src/sage/algebras/exterior_algebra_groebner.pxd b/src/sage/algebras/exterior_algebra_groebner.pxd index 7708930faa8..0c001363130 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pxd +++ b/src/sage/algebras/exterior_algebra_groebner.pxd @@ -38,6 +38,7 @@ cdef class GroebnerStrategy: cpdef CliffordAlgebraElement reduce(self, CliffordAlgebraElement f) cdef CliffordAlgebraElement reduce_single(self, CliffordAlgebraElement f, CliffordAlgebraElement g) + cdef int reduced_gb(self, list G) except -1 # These are the methods that determine the ordering of the monomials. # These must be implemented in subclasses. Declare them as "inline" there. diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index 3a2e73de41c..5a0ce6ccb6c 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -20,6 +20,7 @@ AUTHORS: # http://www.gnu.org/licenses/ #***************************************************************************** +from cysignals.signals cimport sig_check 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, bitset_len) @@ -238,9 +239,9 @@ cdef class GroebnerStrategy: for i, p in enumerate(M.pivots()) if r - Integer(p) not in lead_supports] - def compute_groebner(self): - """ - Compute the reduced ``side`` Gröbner basis for the ideal ``I``. + def compute_groebner(self, reduced=True): + r""" + Compute the (``reduced``) left/right Gröbner basis for the ideal ``I``. EXAMPLES:: @@ -253,9 +254,12 @@ cdef class GroebnerStrategy: (1,) sage: E. = ExteriorAlgebra(QQ) - sage: I = E.ideal([a+b*c]) + sage: I = E.ideal([a+b*c], side="left") + sage: I.groebner_basis() # indirect doctest + (b*c + a,) + sage: I = E.ideal([a+b*c], side="twosided") sage: I.groebner_basis() # indirect doctest - (b*c + a, a*c, a*b, a*d) + (a*b, a*c, b*c + a, a*d) """ cdef FrozenBitset p0, p1 cdef long deg @@ -264,6 +268,7 @@ cdef class GroebnerStrategy: cdef GBElement f0, f1 cdef list G = [self.build_elt(f) for f in self.ideal.gens() if f] # Remove 0s cdef list Gp + cdef CliffordAlgebraElement f if self.side == 2 and not self.homogeneous: # Add in all S-poly times positive degree monomials @@ -286,6 +291,7 @@ cdef class GroebnerStrategy: P[deg] = [(f0, f1)] while P: + sig_check() Pp = P.pop(min(P)) # The selection: lowest lcm degree Gp = self.reduction(Pp, G) G.extend(Gp) @@ -304,29 +310,71 @@ cdef class GroebnerStrategy: G.sort(key=lambda x: ( x).lsi) - for i in range(len(G)): - G[i] = ( G[i]).elt + if not reduced: + self.groebner_basis = tuple([( f0).elt for f0 in G if ( f0).elt]) + return + self.reduced_gb(G) + + cdef int reduced_gb(self, list G) except -1: + """ + Convert the Gröbner basis ``G`` into a reduced Gröbner basis. + """ + cdef Py_ssize_t i, j, k + cdef GBElement f0, f1 # 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(i+1,n)) cdef tuple supp cdef bint did_reduction cdef FrozenBitset lm, s + cdef Integer r + cdef Py_ssize_t num_zeros = 0 + cdef Py_ssize_t n = len(G) + cdef set pairs = set((i, j) for i in range(n) for j in range(n) if i != j) + while pairs: + sig_check() i,j = pairs.pop() - f = G[i] - g = G[j] + f0 = G[i] + f1 = G[j] # We perform the classical reduction algorithm here on each pair # TODO: Make this faster by using the previous technique? - fp = self.reduce_single(f, g) - if f != fp: - G[i] = fp - if not fp: - pairs.difference_update((k, i) for k in range(n)) - else: + f = self.reduce_single(f0.elt, f1.elt) + if f0.elt != f: + if f: + G[i] = self.build_elt(f) pairs.update((k, i) for k in range(n) if k != i) + else: + G[i] = GBElement(f, FrozenBitset(), Integer(2)**self.rank + 1) + num_zeros += 1 + pairs.difference_update((k, i) for k in range(n) if k != i) + pairs.difference_update((i, k) for k in range(n) if k != i) + + G.sort(key=lambda x: ( x).lsi) + for i in range(len(G)-num_zeros): + f0 = G[i] + if f0.elt: + inv = ~f0.elt[f0.ls] + for key in f0.elt._monomial_coefficients: + f0.elt._monomial_coefficients[key] *= inv + self.groebner_basis = tuple([f0.elt for f0 in G[:len(G)-num_zeros]]) + return 0 - self.groebner_basis = tuple([~f[self.leading_support(f)] * f for f in G if f]) + def reduce_computed_gb(self): + """ + Convert the computed Gröbner basis to a reduced Gröbner basis. + + sage: E. = ExteriorAlgebra(QQ) + sage: I = E.ideal([x+y*z]) + sage: I.groebner_basis(reduced=False) + (x*y, x*z, y*z + x, y*z + x, x*y*z) + sage: I._groebner_strategy.reduce_computed_gb() + sage: I._groebner_strategy.groebner_basis + (x*y, x*z, y*z + x) + """ + if self.groebner_basis == [(None,)]: + raise ValueError("Gröbner basis not yet computed") + cdef list G = [self.build_elt(f) for f in self.groebner_basis] + self.reduced_gb(G) cpdef CliffordAlgebraElement reduce(self, CliffordAlgebraElement f): """ @@ -341,23 +389,28 @@ cdef class GroebnerStrategy: ....: 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.reduce(a*b*e) + sage: GB = I.groebner_basis() + sage: I._groebner_strategy.reduce(a*b*e) a*b*e - sage: I.reduce(b*d*e) + sage: I._groebner_strategy.reduce(b*d*e) a*b*d - a*b*e + a*d*e - sage: I.reduce(c*d*e) + sage: I._groebner_strategy.reduce(c*d*e) a*c*d - a*c*e + a*d*e - sage: I.reduce(a*b*c*d*e) + sage: I._groebner_strategy.reduce(a*b*c*d*e) + 0 + sage: I._groebner_strategy.reduce(a*b*c*d) 0 - sage: I.reduce(a*b*c*d) + sage: I._groebner_strategy.reduce(E.zero()) 0 """ + if not f: + return f for g in self.groebner_basis: f = self.reduce_single(f, g) return f cdef CliffordAlgebraElement reduce_single(self, CliffordAlgebraElement f, CliffordAlgebraElement g): - """ + r""" Reduce ``f`` by ``g``. .. TODO:: @@ -376,7 +429,7 @@ cdef class GroebnerStrategy: did_reduction = True mon = build_monomial(self.E, s - lm) if self.side == 0: - gp = mon * g.elt + gp = mon * g f -= f[s] / gp[s] * gp else: gp = g * mon From 140ac6c9b2fced6e9ff87336040d3a906a248a61 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Mon, 8 Aug 2022 14:45:26 +0900 Subject: [PATCH 24/27] Migration of RAAG cohomology to Cython. --- .../algebras/clifford_algebra_element.pxd | 3 + .../algebras/clifford_algebra_element.pyx | 60 +++++++++++++++++++ src/sage/groups/raag.py | 50 +--------------- 3 files changed, 65 insertions(+), 48 deletions(-) diff --git a/src/sage/algebras/clifford_algebra_element.pxd b/src/sage/algebras/clifford_algebra_element.pxd index f7a25582d54..a274f65205a 100644 --- a/src/sage/algebras/clifford_algebra_element.pxd +++ b/src/sage/algebras/clifford_algebra_element.pxd @@ -10,3 +10,6 @@ cdef class CliffordAlgebraElement(IndexedFreeModuleElement): cdef class ExteriorAlgebraElement(CliffordAlgebraElement): pass +cdef class CohomologyRAAGElement(CliffordAlgebraElement): + pass + diff --git a/src/sage/algebras/clifford_algebra_element.pyx b/src/sage/algebras/clifford_algebra_element.pyx index ac74f313f76..dc24c92ae03 100644 --- a/src/sage/algebras/clifford_algebra_element.pyx +++ b/src/sage/algebras/clifford_algebra_element.pyx @@ -607,3 +607,63 @@ cdef class ExteriorAlgebraElement(CliffordAlgebraElement): """ return (self.transpose() * other).constant_coefficient() +cdef class CohomologyRAAGElement(CliffordAlgebraElement): + """ + An element in the cohomology of a right-angled Artin group. + + .. SEEALSO:: + + :class:`~sage.groups.raag.CohomologyRAAG` + """ + cdef _mul_(self, other): + """ + Return ``self`` multiplied by ``other``. + + EXAMPLES:: + + sage: C4 = graphs.CycleGraph(4) + sage: A = groups.misc.RightAngledArtin(C4) + sage: H = A.cohomology() + sage: b = sum(H.basis()) + sage: b * b + 2*e0*e2 + 2*e1*e3 + 2*e0 + 2*e1 + 2*e2 + 2*e3 + 1 + """ + zero = self._parent._base.zero() + cdef frozenset I = frozenset(self._parent._indices) + cdef dict d = {} + cdef list t + cdef tuple tp + cdef tuple ml, mr + cdef Py_ssize_t pos, i, j + cdef bint negate + + for ml, cl in self._monomial_coefficients.items(): + for mr, cr in other._monomial_coefficients.items(): + # Create the next term + tp = tuple(sorted(mr + ml)) + if any(tp[i] == tp[i+1] for i in range(len(tp)-1)): # e_i ^ e_i = 0 + continue + if tp not in I: # not an independent set, so this term is also 0 + continue + + t = list(mr) + negate = False + for i in reversed(ml): + pos = 0 + for j in t: + assert i != j + if i < j: + break + pos += 1 + negate = not negate + t.insert(pos, i) + + if negate: + cr = -cr + + d[tp] = d.get(tp, zero) + cl * cr + if d[tp] == zero: + del d[tp] + + return self.__class__(self._parent, d) + diff --git a/src/sage/groups/raag.py b/src/sage/groups/raag.py index 7e0dffe4419..ecb36fc8b18 100644 --- a/src/sage/groups/raag.py +++ b/src/sage/groups/raag.py @@ -39,7 +39,7 @@ from sage.combinat.free_module import CombinatorialFreeModule from sage.categories.fields import Fields from sage.categories.algebras_with_basis import AlgebrasWithBasis -from sage.algebras.clifford_algebra import CliffordAlgebraElement +from sage.algebras.clifford_algebra_element import CohomologyRAAGElement from sage.typeset.ascii_art import ascii_art from sage.typeset.unicode_art import unicode_art @@ -853,54 +853,8 @@ def degree_on_basis(self, I): """ return len(I) - class Element(CliffordAlgebraElement): + class Element(CohomologyRAAGElement): """ An element in the cohomology ring of a right-angled Artin group. """ - def _mul_(self, other): - """ - Return ``self`` multiplied by ``other``. - - EXAMPLES:: - - sage: C4 = graphs.CycleGraph(4) - sage: A = groups.misc.RightAngledArtin(C4) - sage: H = A.cohomology() - sage: b = sum(H.basis()) - sage: b * b - 2*e0*e2 + 2*e1*e3 + 2*e0 + 2*e1 + 2*e2 + 2*e3 + 1 - """ - zero = self.parent().base_ring().zero() - I = self.parent()._indices - d = {} - - for ml,cl in self: - for mr,cr in other: - # Create the next term - t = list(mr) - for i in reversed(ml): - pos = 0 - for j in t: - if i == j: - pos = None - break - if i < j: - break - pos += 1 - cr = -cr - if pos is None: - t = None - break - t.insert(pos, i) - - if t is None: # The next term is 0, move along - continue - - t = tuple(t) - if t not in I: # not an independent set, so this term is also 0 - continue - d[t] = d.get(t, zero) + cl * cr - if d[t] == zero: - del d[t] - return self.__class__(self.parent(), d) From b2279734577506e0e4983e07dd7488c7106573dc Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Tue, 9 Aug 2022 09:19:42 +0900 Subject: [PATCH 25/27] Speeding up multiplication even further. --- .../algebras/clifford_algebra_element.pyx | 35 ++++++++++++++++--- 1 file changed, 30 insertions(+), 5 deletions(-) diff --git a/src/sage/algebras/clifford_algebra_element.pyx b/src/sage/algebras/clifford_algebra_element.pyx index dc24c92ae03..f2a6e09ad86 100644 --- a/src/sage/algebras/clifford_algebra_element.pyx +++ b/src/sage/algebras/clifford_algebra_element.pyx @@ -341,25 +341,50 @@ cdef class ExteriorAlgebraElement(CliffordAlgebraElement): -x*y*z*w sage: (z * w) * (x * y) x*y*z*w + + sage: E. = ExteriorAlgebra(QQ) + sage: r = sum(E.basis()) + sage: r*r + 4*a*b*c*d + 4*a*b*c + 4*a*b*d + 4*a*c*d + 4*b*c*d + + 2*a*b + 2*a*c + 2*a*d + 2*b*c + 2*b*d + 2*c*d + + 2*a + 2*b + 2*c + 2*d + 1 """ cdef Parent P = self._parent zero = P._base.zero() - cdef dict d = {} + cdef dict d cdef Py_ssize_t n = P.ngens() cdef ExteriorAlgebraElement rhs = other + cdef list to_remove cdef FrozenBitset ml, mr, t cdef Py_ssize_t num_cross, tot_cross, i, j + ml = FrozenBitset() + + if ml in self._monomial_coefficients: + const_coeff = self._monomial_coefficients[ml] + d = dict(rhs._monomial_coefficients) # Make a shallow copy + to_remove = [] + if const_coeff != P._base.one(): + for k in d: + d[k] *= const_coeff + if not d[k]: # there might be zero divisors + to_remove.append(k) + for k in to_remove: + del d[k] + else: + d = {} + for ml,cl in self._monomial_coefficients.items(): # ml for "monomial on the left" + if not ml: # We already handled the trivial element + continue for mr,cr in rhs._monomial_coefficients.items(): # mr for "monomial on the right" - if ml.intersection(mr): - # if they intersect nontrivially, move along. - continue - if not mr: t = ml else: + if not ml.isdisjoint(mr): + # if they intersect nontrivially, move along. + continue t = ml._union(mr) it = iter(mr) j = next(it) From e7e84b73a87aac57620190aa776448f422f0c6cd Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Tue, 9 Aug 2022 11:05:26 +0900 Subject: [PATCH 26/27] Special casing multiplication by a term. --- .../algebras/clifford_algebra_element.pxd | 4 +- .../algebras/clifford_algebra_element.pyx | 311 +++++++++++++++++- .../algebras/exterior_algebra_groebner.pyx | 29 +- 3 files changed, 322 insertions(+), 22 deletions(-) diff --git a/src/sage/algebras/clifford_algebra_element.pxd b/src/sage/algebras/clifford_algebra_element.pxd index a274f65205a..64ac7253351 100644 --- a/src/sage/algebras/clifford_algebra_element.pxd +++ b/src/sage/algebras/clifford_algebra_element.pxd @@ -3,9 +3,11 @@ Clifford algebra elements """ from sage.modules.with_basis.indexed_element cimport IndexedFreeModuleElement +from sage.data_structures.bitset cimport FrozenBitset cdef class CliffordAlgebraElement(IndexedFreeModuleElement): - pass + cdef CliffordAlgebraElement _mul_self_term(self, FrozenBitset supp, coeff) + cdef CliffordAlgebraElement _mul_term_self(self, FrozenBitset supp, coeff) cdef class ExteriorAlgebraElement(CliffordAlgebraElement): pass diff --git a/src/sage/algebras/clifford_algebra_element.pyx b/src/sage/algebras/clifford_algebra_element.pyx index f2a6e09ad86..e065b6b70b9 100644 --- a/src/sage/algebras/clifford_algebra_element.pyx +++ b/src/sage/algebras/clifford_algebra_element.pyx @@ -19,8 +19,9 @@ AUTHORS: #***************************************************************************** from sage.structure.parent cimport Parent -from sage.data_structures.bitset cimport FrozenBitset, Bitset +from sage.data_structures.bitset cimport Bitset from sage.algebras.weyl_algebra import repr_from_monomials +from sage.data_structures.blas_dict cimport scal from copy import copy cdef class CliffordAlgebraElement(IndexedFreeModuleElement): @@ -95,8 +96,25 @@ cdef class CliffordAlgebraElement(IndexedFreeModuleElement): cdef dict next_level, cur, d = {} cdef FrozenBitset ml, mr, t cdef Py_ssize_t i, j + cdef CliffordAlgebraElement rhs = other - for ml,cl in self: + # Special case when multiplying by 0 + if not self._monomial_coefficients: + return self + if not rhs._monomial_coefficients: + return rhs + + # Special case when multiplying by an element of the base ring + if len(self._monomial_coefficients) == 1: + ml, cl = next(iter(self._monomial_coefficients.items())) + if ml.isempty(): + return rhs._mul_term_self(ml, cl) + if len(rhs._monomial_coefficients) == 1: + mr, cr = next(iter(self._monomial_coefficients.items())) + if mr.isempty(): + return self._mul_self_term(mr, cr) + + for ml, cl in self: # Distribute the current term ``cl`` * ``ml`` over ``other``. cur = copy(other._monomial_coefficients) # The current distribution of the term for i in reversed(ml): @@ -150,6 +168,76 @@ cdef class CliffordAlgebraElement(IndexedFreeModuleElement): return self.__class__(self.parent(), d) + cdef CliffordAlgebraElement _mul_self_term(self, FrozenBitset supp, coeff): + r""" + Multiply ``self * term`` with the ``term`` having support ``supp`` + and coefficient ``coeff``. + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: r = sum(E.basis()) + sage: x * y # indirect doctest + x*y + sage: y * x # indirect doctest + -x*y + sage: r * x # indirect doctest + x*y*z - x*y - x*z + x + sage: r * -x # indirect doctest + -x*y*z + x*y + x*z - x + sage: r * (2*x) # indirect doctest + 2*x*y*z - 2*x*y - 2*x*z + 2*x + sage: r * y # indirect doctest + -x*y*z + x*y - y*z + y + sage: r * z # indirect doctest + x*y*z + x*z + y*z + z + sage: r * (x*y) # indirect doctest + x*y*z + x*y + sage: r * (-x*y) # indirect doctest + -x*y*z - x*y + sage: r * (x*y*z) # indirect doctest + x*y*z + sage: r * 1 == r # indirect doctest + True + sage: r * -1 == -r # indirect doctest + True + sage: r * 2 # indirect doctest + 2*x*y*z + 2*x*y + 2*x*z + 2*y*z + 2*x + 2*y + 2*z + 2 + """ + cdef dict d + cdef list to_remove + cdef Py_ssize_t num_cross, tot_cross, i, j + cdef FrozenBitset ml + + if supp.isempty(): # Multiplication by a base ring element + if coeff == self._parent._base.one(): + return self + if coeff == -self._parent._base.one(): + return self._neg_() + + return type(self)(self._parent, + scal(coeff, self._monomial_coefficients, + factor_on_left=False)) + + return type(self)(self._parent, {supp: coeff}) * self + + cdef CliffordAlgebraElement _mul_term_self(self, FrozenBitset supp, coeff): + r""" + Multiply ``term * self`` with the ``term`` having support ``supp`` + and coefficient ``coeff``. + """ + if supp.isempty(): # Multiplication by a base ring element + if coeff == self._parent._base.one(): + return self + if coeff == -self._parent._base.one(): + return self._neg_() + + return type(self)(self._parent, + scal(coeff, self._monomial_coefficients, + factor_on_left=True)) + + return type(self)(self._parent, {supp: coeff}) * self + def list(self): """ Return the list of monomials and their coefficients in ``self`` @@ -352,15 +440,30 @@ cdef class ExteriorAlgebraElement(CliffordAlgebraElement): cdef Parent P = self._parent zero = P._base.zero() cdef dict d - cdef Py_ssize_t n = P.ngens() cdef ExteriorAlgebraElement rhs = other cdef list to_remove cdef FrozenBitset ml, mr, t - cdef Py_ssize_t num_cross, tot_cross, i, j + cdef Py_ssize_t n, num_cross, tot_cross, i, j - ml = FrozenBitset() + # Special case: one of them is zero + if not self._monomial_coefficients: + return self + if not rhs._monomial_coefficients: + return rhs + + # Special case: other is a single term + if len(rhs._monomial_coefficients) == 1: + mr, cr = next(iter(rhs._monomial_coefficients.items())) + return self._mul_self_term(mr, cr) + # Special case: self is a single term + if len(self._monomial_coefficients) == 1: + ml, cl = next(iter(self._monomial_coefficients.items())) + return rhs._mul_term_self(ml, cl) + + # Do some special processing for the constant monomial in ml + ml = FrozenBitset() if ml in self._monomial_coefficients: const_coeff = self._monomial_coefficients[ml] d = dict(rhs._monomial_coefficients) # Make a shallow copy @@ -375,11 +478,12 @@ cdef class ExteriorAlgebraElement(CliffordAlgebraElement): else: d = {} - for ml,cl in self._monomial_coefficients.items(): # ml for "monomial on the left" - if not ml: # We already handled the trivial element + n = P.ngens() + for ml, cl in self._monomial_coefficients.items(): # ml for "monomial on the left" + if ml.isempty(): # We already handled the trivial element continue for mr,cr in rhs._monomial_coefficients.items(): # mr for "monomial on the right" - if not mr: + if mr.isempty(): t = ml else: if not ml.isdisjoint(mr): @@ -402,12 +506,199 @@ cdef class ExteriorAlgebraElement(CliffordAlgebraElement): if tot_cross % 2: cr = -cr - d[t] = d.get(t, zero) + cl * cr - if not d[t]: + val = d.get(t, zero) + cl * cr + if not val: del d[t] + else: + d[t] = val return self.__class__(P, d) + cdef CliffordAlgebraElement _mul_self_term(self, FrozenBitset supp, coeff): + r""" + Multiply ``self * term`` with the ``term`` having support ``supp`` + and coefficient ``coeff``. + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: r = sum(E.basis()) + sage: x * y # indirect doctest + x*y + sage: y * x # indirect doctest + -x*y + sage: r * x # indirect doctest + x*y*z - x*y - x*z + x + sage: r * -x # indirect doctest + -x*y*z + x*y + x*z - x + sage: r * (2*x) # indirect doctest + 2*x*y*z - 2*x*y - 2*x*z + 2*x + sage: r * y # indirect doctest + -x*y*z + x*y - y*z + y + sage: r * z # indirect doctest + x*y*z + x*z + y*z + z + sage: r * (x*y) # indirect doctest + x*y*z + x*y + sage: r * (-x*y) # indirect doctest + -x*y*z - x*y + sage: r * (x*y*z) # indirect doctest + x*y*z + sage: r * 1 == r # indirect doctest + True + sage: r * -1 == -r # indirect doctest + True + sage: r * 2 # indirect doctest + 2*x*y*z + 2*x*y + 2*x*z + 2*y*z + 2*x + 2*y + 2*z + 2 + """ + cdef dict d + cdef list to_remove + cdef Py_ssize_t num_cross, tot_cross, i, j + cdef FrozenBitset ml + + if supp.isempty(): # Multiplication by a base ring element + if coeff == self._parent._base.one(): + return self + if coeff == -self._parent._base.one(): + return self._neg_() + + return type(self)(self._parent, + scal(coeff, self._monomial_coefficients, + factor_on_left=False)) + + n = self._parent.ngens() + d = {} + for ml, cl in self._monomial_coefficients.items(): # ml for "monomial on the left" + if not ml.isdisjoint(supp): + # if they intersect nontrivially, move along. + continue + t = ml._union(supp) + it = iter(supp) + j = next(it) + + num_cross = 0 # keep track of the number of signs + tot_cross = 0 + for i in ml: + while i > j: + num_cross += 1 + try: + j = next(it) + except StopIteration: + j = n + 1 + tot_cross += num_cross + if tot_cross % 2: + d[t] = -cl + else: + d[t] = cl + + if coeff == -self._parent._base.one(): + for k in d: + d[k] = -d[k] + elif coeff != self._parent._base.one(): + to_remove = [] + for k in d: + d[k] *= coeff + if not d[k]: # there might be zero divisors + to_remove.append(k) + for k in to_remove: + del d[k] + return type(self)(self._parent, d) + + cdef CliffordAlgebraElement _mul_term_self(self, FrozenBitset supp, coeff): + r""" + Multiply ``term * self`` with the ``term`` having support ``supp`` + and coefficient ``coeff``. + + EXAMPLES:: + + sage: E. = ExteriorAlgebra(QQ) + sage: r = sum(E.basis()) + sage: x * r # indirect doctest + x*y*z + x*y + x*z + x + sage: (-x) * r # indirect doctest + -x*y*z - x*y - x*z - x + sage: (2*x) * r # indirect doctest + 2*x*y*z + 2*x*y + 2*x*z + 2*x + sage: y * r # indirect doctest + -x*y*z - x*y + y*z + y + sage: z * r # indirect doctest + x*y*z - x*z - y*z + z + sage: (x*y) * r # indirect doctest + x*y*z + x*y + sage: (-x*y) * r # indirect doctest + -x*y*z - x*y + sage: (x*y*z) * r # indirect doctest + x*y*z + sage: 1 * r == r # indirect doctest + True + sage: -1 * r == -r # indirect doctest + True + sage: 2 * r # indirect doctest + 2*x*y*z + 2*x*y + 2*x*z + 2*y*z + 2*x + 2*y + 2*z + 2 + """ + cdef dict d + cdef list to_remove + cdef Py_ssize_t n, num_cross, tot_cross, i, j + cdef FrozenBitset mr, t + + if supp.isempty(): # Multiplication by a base ring element + if coeff == self._parent._base.one(): + return self + if coeff == -self._parent._base.one(): + return self._neg_() + + return type(self)(self._parent, + scal(coeff, self._monomial_coefficients, + factor_on_left=True)) + + n = self._parent.ngens() + d = {} + mr = FrozenBitset() + # We need to special case the constant coefficient + const_coeff = None + if mr in self._monomial_coefficients: + const_coeff = self._monomial_coefficients.pop(mr) + d[supp] = const_coeff + + for mr, cr in self._monomial_coefficients.items(): # mr for "monomial on the right" + if not supp.isdisjoint(mr): + # if they intersect nontrivially, move along. + continue + t = supp._union(mr) + it = iter(mr) + j = next(it) # We assume mr is non-empty here + + num_cross = 0 # keep track of the number of signs + tot_cross = 0 + for i in supp: + while i > j: + num_cross += 1 + try: + j = next(it) + except StopIteration: + j = n + 1 + tot_cross += num_cross + if tot_cross % 2: + d[t] = -cr + else: + d[t] = cr + + if coeff == -self._parent._base.one(): + for k in d: + d[k] = -d[k] + elif coeff != self._parent._base.one(): + to_remove = [] + for k in d: + d[k] = coeff * d[k] # This will work for non-commutative base rings + if not d[k]: # there might be zero divisors + to_remove.append(k) + for k in to_remove: + del d[k] + + # Add back the constant coefficient since we removed it for the special case + if const_coeff is not None: + self._monomial_coefficients[FrozenBitset()] = const_coeff + return type(self)(self._parent, d) + def reduce(self, I, left=True): r""" Reduce ``self`` with respect to the elements in ``I``. diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index 5a0ce6ccb6c..ad514216abb 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -26,6 +26,7 @@ from sage.data_structures.bitset_base cimport (bitset_t, bitset_init, bitset_fir bitset_next, bitset_set_to, bitset_len) from sage.structure.parent cimport Parent from sage.structure.richcmp cimport richcmp +from sage.data_structures.blas_dict cimport iaxpy, remove_zeros cdef inline long degree(FrozenBitset X): """ @@ -143,9 +144,9 @@ cdef class GroebnerStrategy: This assumes the leading support is ``f.ls._union(t)``. """ - ret = f.elt * build_monomial(self.E, t) + ret = f.elt._mul_self_term(self.E, self.E._base.one()) cdef FrozenBitset ls = f.ls._union(t) - return GBElement(ret, ls, self.bitset_to_int(ls)) + return GBElement( ret, ls, self.bitset_to_int(ls)) cdef inline GBElement prod_term_GB(self, FrozenBitset t, GBElement f): """ @@ -155,9 +156,9 @@ cdef class GroebnerStrategy: This assumes the leading support is ``f.ls._union(t)``. """ - ret = build_monomial(self.E, t) * f.elt + ret = f.elt._mul_term_self(t, self.E._base.one()) cdef FrozenBitset ls = f.ls._union(t) - return GBElement(ret, ls, self.bitset_to_int(ls)) + return GBElement( ret, ls, self.bitset_to_int(ls)) cdef inline bint build_S_poly(self, GBElement f, GBElement g): r""" @@ -195,7 +196,8 @@ cdef class GroebnerStrategy: if self.side == 2 and not self.homogeneous: # Add in all S-poly times positive degree monomials - additions = set(( f).elt * t for t in self.E.basis() for f in L) + one = self.E._base.one() + additions = set(( f).elt._mul_self_term(t, one) for t in self.E._indices for f in L) L.update(self.build_elt(f) for f in additions if f) cdef set done = set(( f).ls for f in L) @@ -272,7 +274,8 @@ cdef class GroebnerStrategy: if self.side == 2 and not self.homogeneous: # Add in all S-poly times positive degree monomials - additions = set(( f0).elt * t for t in self.E.basis() for f0 in G) + one = self.E._base.one() + additions = set(( f0).elt._mul_self_term(t, one) for t in self.E._indices for f0 in G) G.extend(self.build_elt(f) for f in additions if f) cdef Py_ssize_t n = len(G) @@ -420,20 +423,24 @@ cdef class GroebnerStrategy: cdef FrozenBitset lm = self.leading_support(g), s cdef bint did_reduction = True cdef tuple supp + cdef CliffordAlgebraElement gp + one = self.E._base.one() while did_reduction: supp = tuple(f._monomial_coefficients) did_reduction = False for s in supp: if lm <= s: did_reduction = True - mon = build_monomial(self.E, s - lm) if self.side == 0: - gp = mon * g - f -= f[s] / gp[s] * gp + gp = g._mul_term_self(s - lm, one) else: - gp = g * mon - f -= f[s] / gp[s] * gp + gp = g._mul_self_term(s - lm, one) + coeff = f[s] / gp._monomial_coefficients[s] + for k in gp._monomial_coefficients: + gp._monomial_coefficients[k] *= coeff + remove_zeros(gp._monomial_coefficients) + f -= gp break return f From 36b6b201b9be9350b284653cd21dcac55d461b71 Mon Sep 17 00:00:00 2001 From: Travis Scrimshaw Date: Tue, 9 Aug 2022 12:18:51 +0900 Subject: [PATCH 27/27] Doing reduction in place; being careful about duplicates. --- src/sage/algebras/clifford_algebra.py | 2 +- .../algebras/exterior_algebra_groebner.pxd | 2 +- .../algebras/exterior_algebra_groebner.pyx | 143 ++++++++++++++---- 3 files changed, 112 insertions(+), 35 deletions(-) diff --git a/src/sage/algebras/clifford_algebra.py b/src/sage/algebras/clifford_algebra.py index 85b6e902f0d..0def7ae087b 100644 --- a/src/sage/algebras/clifford_algebra.py +++ b/src/sage/algebras/clifford_algebra.py @@ -2946,7 +2946,7 @@ def groebner_basis(self, term_order=None, reduced=True): sage: E. = ExteriorAlgebra(QQ) sage: I = E.ideal([x+y*z]) sage: I.groebner_basis(reduced=False) - (x*y, x*z, y*z + x, y*z + x, x*y*z) + (x*y, x*z, y*z + x, x*y*z) sage: I.groebner_basis(reduced=True) (x*y, x*z, y*z + x) diff --git a/src/sage/algebras/exterior_algebra_groebner.pxd b/src/sage/algebras/exterior_algebra_groebner.pxd index 0c001363130..d00fb2e8560 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pxd +++ b/src/sage/algebras/exterior_algebra_groebner.pxd @@ -37,7 +37,7 @@ cdef class GroebnerStrategy: cdef list reduction(self, list P, list G) cpdef CliffordAlgebraElement reduce(self, CliffordAlgebraElement f) - cdef CliffordAlgebraElement reduce_single(self, CliffordAlgebraElement f, CliffordAlgebraElement g) + cdef bint reduce_single(self, CliffordAlgebraElement f, CliffordAlgebraElement g) except -1 cdef int reduced_gb(self, list G) except -1 # These are the methods that determine the ordering of the monomials. diff --git a/src/sage/algebras/exterior_algebra_groebner.pyx b/src/sage/algebras/exterior_algebra_groebner.pyx index ad514216abb..98e338d3468 100644 --- a/src/sage/algebras/exterior_algebra_groebner.pyx +++ b/src/sage/algebras/exterior_algebra_groebner.pyx @@ -25,8 +25,9 @@ from sage.libs.gmp.mpz cimport mpz_sizeinbase, mpz_setbit, mpz_tstbit, mpz_cmp_s from sage.data_structures.bitset_base cimport (bitset_t, bitset_init, bitset_first, bitset_next, bitset_set_to, bitset_len) from sage.structure.parent cimport Parent -from sage.structure.richcmp cimport richcmp -from sage.data_structures.blas_dict cimport iaxpy, remove_zeros +from sage.structure.richcmp cimport richcmp, rich_to_bool +from sage.data_structures.blas_dict cimport iaxpy +from copy import copy cdef inline long degree(FrozenBitset X): """ @@ -42,16 +43,59 @@ cdef inline CliffordAlgebraElement build_monomial(Parent E, FrozenBitset supp): return E.element_class(E, {supp: E._base.one()}) cdef class GBElement: + """ + Helper class for storing an element with its leading support both as + a :class:`FrozenBitset` and an :class:`Integer`. + """ def __init__(self, CliffordAlgebraElement x, FrozenBitset ls, Integer n): + """ + Initialize ``self``. + + EXAMPLES:: + + sage: from sage.algebras.exterior_algebra_groebner import GBElement + sage: E. = ExteriorAlgebra(QQ) + sage: X = GBElement(a, a.leading_support(), 1) + sage: TestSuite(X).run(skip="_test_pickling") + """ self.elt = x self.ls = ls self.lsi = n def __hash__(self): + """ + Return the hash of ``self``. + + EXAMPLES:: + + sage: from sage.algebras.exterior_algebra_groebner import GBElement + sage: E. = ExteriorAlgebra(QQ) + sage: X = GBElement(a, a.leading_support(), 1) + sage: hash(X) == 1 + True + """ return int(self.lsi) def __richcmp__(self, other, int op): - return richcmp(self.elt, ( other).elt, op) + """ + Rich compare ``self`` with ``other`` by ``op``. + + EXAMPLES:: + + sage: from sage.algebras.exterior_algebra_groebner import GBElement + sage: E. = ExteriorAlgebra(QQ) + sage: X = GBElement(a, a.leading_support(), 1) + sage: Y = GBElement(a*b, (a*b).leading_support(), 3) + sage: X == X + True + sage: X == Y + False + sage: X != Y + True + """ + if self is other: + return rich_to_bool(op, 0) + return richcmp(self.elt, ( other).elt, op) cdef class GroebnerStrategy: """ @@ -268,15 +312,38 @@ cdef class GroebnerStrategy: cdef Py_ssize_t i, j, k cdef set additions cdef GBElement f0, f1 - cdef list G = [self.build_elt(f) for f in self.ideal.gens() if f] # Remove 0s - cdef list Gp + cdef list G = [], Gp + cdef dict constructed = {} cdef CliffordAlgebraElement f + for f in self.ideal.gens(): + if not f: # Remove 0s + continue + f0 = self.build_elt(f) + if f0.lsi in constructed: + if f0 in constructed[f0.lsi]: # Already there + continue + constructed[f0.lsi].add(f0) + else: + constructed[f0.lsi] = set([f0]) + G.append(f0) + if self.side == 2 and not self.homogeneous: # Add in all S-poly times positive degree monomials one = self.E._base.one() - additions = set(( f0).elt._mul_self_term(t, one) for t in self.E._indices for f0 in G) - G.extend(self.build_elt(f) for f in additions if f) + for t in self.E._indices: + for f0 in G: + f = f0.elt._mul_self_term(t, one) + if not f: + continue + f1 = self.build_elt(f) + if f1.lsi in constructed: + if f1 in constructed[f1.lsi]: # Already there + continue + constructed[f1.lsi].add(f1) + else: + constructed[f1.lsi] = set([f1]) + G.append(f1) cdef Py_ssize_t n = len(G) cdef dict P = {} @@ -297,7 +364,16 @@ cdef class GroebnerStrategy: sig_check() Pp = P.pop(min(P)) # The selection: lowest lcm degree Gp = self.reduction(Pp, G) - G.extend(Gp) + # Add the elements Gp to G when a new element is found + for f0 in Gp: + if f0.lsi in constructed: + if f0 in constructed[f0.lsi]: # Already there + continue + constructed[f0.lsi].add(f0) + else: + constructed[f0.lsi] = set([f0]) + G.append(f0) + # Find the degress of the new pairs for j in range(n, len(G)): f1 = G[j] p1 = f1.ls @@ -339,15 +415,15 @@ cdef class GroebnerStrategy: i,j = pairs.pop() f0 = G[i] f1 = G[j] + assert f0.elt._monomial_coefficients is not f1.elt._monomial_coefficients, (i,j) # We perform the classical reduction algorithm here on each pair # TODO: Make this faster by using the previous technique? - f = self.reduce_single(f0.elt, f1.elt) - if f0.elt != f: - if f: - G[i] = self.build_elt(f) + if self.reduce_single(f0.elt, f1.elt): + if f0.elt: + G[i] = self.build_elt(f0.elt) pairs.update((k, i) for k in range(n) if k != i) else: - G[i] = GBElement(f, FrozenBitset(), Integer(2)**self.rank + 1) + G[i] = GBElement(f0.elt, FrozenBitset(), Integer(2)**self.rank + 1) num_zeros += 1 pairs.difference_update((k, i) for k in range(n) if k != i) pairs.difference_update((i, k) for k in range(n) if k != i) @@ -369,7 +445,7 @@ cdef class GroebnerStrategy: sage: E. = ExteriorAlgebra(QQ) sage: I = E.ideal([x+y*z]) sage: I.groebner_basis(reduced=False) - (x*y, x*z, y*z + x, y*z + x, x*y*z) + (x*y, x*z, y*z + x, x*y*z) sage: I._groebner_strategy.reduce_computed_gb() sage: I._groebner_strategy.groebner_basis (x*y, x*z, y*z + x) @@ -408,41 +484,42 @@ cdef class GroebnerStrategy: """ if not f: return f + # Make a copy to mutate + f = type(f)(f._parent, copy(f._monomial_coefficients)) for g in self.groebner_basis: - f = self.reduce_single(f, g) + self.reduce_single(f, g) return f - cdef CliffordAlgebraElement reduce_single(self, CliffordAlgebraElement f, CliffordAlgebraElement g): + cdef bint reduce_single(self, CliffordAlgebraElement f, CliffordAlgebraElement g) except -1: r""" Reduce ``f`` by ``g``. - .. TODO:: + .. WARNING:: - Optimize this by doing it in-place and changing the underlying dict of ``f``. + This modifies the element ``f``. """ - cdef FrozenBitset lm = self.leading_support(g), s - cdef bint did_reduction = True + cdef FrozenBitset lm = self.leading_support(g), s, t + cdef bint did_reduction = True, was_reduced=False cdef tuple supp cdef CliffordAlgebraElement gp one = self.E._base.one() while did_reduction: - supp = tuple(f._monomial_coefficients) did_reduction = False - for s in supp: - if lm <= s: + for s in f._monomial_coefficients: + if lm.issubset(s): + t = s did_reduction = True - if self.side == 0: - gp = g._mul_term_self(s - lm, one) - else: - gp = g._mul_self_term(s - lm, one) - coeff = f[s] / gp._monomial_coefficients[s] - for k in gp._monomial_coefficients: - gp._monomial_coefficients[k] *= coeff - remove_zeros(gp._monomial_coefficients) - f -= gp + was_reduced = True break - return f + if did_reduction: + if self.side == 0: + gp = g._mul_term_self(t - lm, one) + else: + gp = g._mul_self_term(t - lm, one) + coeff = f[t] / gp._monomial_coefficients[t] + iaxpy(-coeff, gp._monomial_coefficients, f._monomial_coefficients) + return was_reduced cdef Integer bitset_to_int(self, FrozenBitset X):