gh-37697: Implement the solvable and nil radicals of a finite dimensi…

…onal Lie algebra

    
<!-- ^ Please provide a concise and informative title. -->
<!-- ^ Don't put issue numbers in the title, do this in the PR
description below. -->
<!-- ^ For example, instead of "Fixes #12345" use "Introduce new method
to calculate 1 + 2". -->
<!-- v Describe your changes below in detail. -->
<!-- v Why is this change required? What problem does it solve? -->
<!-- v If this PR resolves an open issue, please link to it here. For
example, "Fixes #12345". -->

We implement the algorithms in de Graaf to compute the solvable and
nilradicals of a finite dimensional Lie algebra. To do so, we also make
further improvements to how subalgebras and ideals of Lie algebras are
constructed and handled.

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

- [x] The title is concise and informative.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on. For example,
-->
<!-- - #12345: short description why this is a dependency -->
<!-- - #34567: ... -->

- #37693 Uses this in the computation
- #37682 Needs the improvements here
    
URL: #37697
Reported by: Travis Scrimshaw
Reviewer(s): Matthias Köppe, Travis Scrimshaw

src/sage/algebras/lie_algebras/lie_algebra_element.pxd

@@ -20,7 +20,8 @@ cdef class LieSubalgebraElementWrapper(LieAlgebraElementWrapper):
 cdef class StructureCoefficientsElement(LieAlgebraMatrixWrapper):
     cpdef bracket(self, right)
     cpdef _bracket_(self, right)
-    cpdef to_vector(self, bint sparse=*)
+    cpdef _vector_(self, bint sparse=*, order=*)
+    cpdef to_vector(self, bint sparse=*, order=*)
     cpdef dict monomial_coefficients(self, bint copy=*)
     # cpdef lift(self)
 

src/sage/algebras/lie_algebras/lie_algebra_element.pyx

@@ -543,7 +543,7 @@ cdef class LieSubalgebraElementWrapper(LieAlgebraElementWrapper):
         x_lift = (<LieSubalgebraElementWrapper> x).value
         return type(self)(self._parent, self.value._bracket_(x_lift))
 
-    def to_vector(self, order=None, sparse=False):
+    def _vector_(self, sparse=False, order=None):
         r"""
         Return the vector in ``g.module()`` corresponding to the
         element ``self`` of ``g`` (where ``g`` is the parent of ``self``).
@@ -573,6 +573,8 @@ cdef class LieSubalgebraElementWrapper(LieAlgebraElementWrapper):
         """
         return self._parent.module()(self.value.to_vector(sparse=sparse))
 
+    to_vector = _vector_
+
     cpdef dict monomial_coefficients(self, bint copy=True):
         r"""
         Return a dictionary whose keys are indices of basis elements
@@ -833,21 +835,34 @@ cdef class StructureCoefficientsElement(LieAlgebraMatrixWrapper):
             if v != zero:
                 yield (I[i], v)
 
-    cpdef to_vector(self, bint sparse=False):
+    cpdef _vector_(self, bint sparse=False, order=None):
         """
         Return ``self`` as a vector.
 
         EXAMPLES::
 
             sage: L.<x,y,z> = LieAlgebra(QQ, {('x','y'): {'z':1}})
-            sage: a = x + 3*y - z/2
-            sage: a.to_vector()
-            (1, 3, -1/2)
+            sage: a = x + 3*y - z/5
+            sage: vector(a)
+            (1, 3, -1/5)
         """
         if sparse:
             return self.value.sparse_vector()
         return self.value
 
+    cpdef to_vector(self, bint sparse=False, order=None):
+        """
+        Return ``self`` as a vector.
+
+        EXAMPLES::
+
+            sage: L.<x,y,z> = LieAlgebra(QQ, {('x','y'): {'z':1}})
+            sage: a = x + 3*y - z/2
+            sage: a.to_vector()
+            (1, 3, -1/2)
+        """
+        return self._vector_(sparse=sparse)
+
     def lift(self):
         """
         Return the lift of ``self`` to the universal enveloping algebra.

src/sage/algebras/lie_algebras/quotient.py

@@ -218,11 +218,14 @@ def __classcall_private__(cls, I, ambient=None, names=None,
         index_set = [i for i in sorted_indices if i not in I_supp]
 
         if names is None:
-            if ambient._names is not None:
-                # ambient has assigned variable names, so use those
+            try:
                 amb_names = dict(zip(sorted_indices, ambient.variable_names()))
                 names = [amb_names[i] for i in index_set]
-        elif isinstance(names, str):
+            except (ValueError, KeyError):
+                # ambient has not assigned variable names
+                # or the names are for the generators rather than the basis
+                names = 'e'
+        if isinstance(names, str):
             if len(index_set) == 1:
                 names = [names]
             else:
@@ -271,6 +274,7 @@ def __init__(self, I, L, names, index_set, category=None):
         self._ambient = L
         self._I = I
         self._sm = sm
+        self._triv_ideal = bool(I.dimension() == 0)
 
         LieAlgebraWithStructureCoefficients.__init__(
             self, L.base_ring(), s_coeff, names, index_set, category=category)
@@ -423,14 +427,23 @@ def from_vector(self, v, order=None, coerce=False):
             sage: el.parent() == Q
             True
 
-        An element from a vector of the ambient module
+        An element from a vector of the ambient module::
 
             sage: el = Q.from_vector([1, 2, 3]); el
             -2*X - Y
             sage: el.parent() == Q
             True
+
+        Check for the trivial ideal::
+
+            sage: L.<x,y,z> = LieAlgebra(GF(3), {('x','z'): {'x':1, 'y':1}, ('y','z'): {'y':1}})
+            sage: I = L.ideal([])
+            sage: Q = L.quotient(I)
+            sage: v = Q.an_element().to_vector()
+            sage: Q.from_vector(v)
+            x + y + z
         """
-        if len(v) == self.ambient().dimension():
+        if not self._triv_ideal and len(v) == self.ambient().dimension():
             return self.retract(self.ambient().from_vector(v))
 
         return super().from_vector(v)

src/sage/algebras/lie_algebras/subalgebra.py

@@ -24,6 +24,7 @@
 from sage.sets.family import Family
 from sage.sets.finite_enumerated_set import FiniteEnumeratedSet
 from sage.structure.parent import Parent
+from sage.structure.element import parent
 from sage.structure.unique_representation import UniqueRepresentation
 
 
@@ -388,7 +389,7 @@ def _an_element_(self):
             sage: S._an_element_()
             X
         """
-        return self.element_class(self, self.lie_algebra_generators()[0])
+        return self.lie_algebra_generators()[0]
 
     def _element_constructor_(self, x):
         """
@@ -425,14 +426,13 @@ def _element_constructor_(self, x):
             sage: S([S(x), S(y)]) == S(L[x, y])
             True
         """
-        try:
-            P = x.parent()
-            if P is self:
-                return x
-            if P == self._ambient:
-                return self.retract(x)
-        except AttributeError:
-            pass
+        P = parent(x)
+        if P is self:
+            return x
+        if P == self._ambient:
+            return self.retract(x)
+        if isinstance(P, type(self)) and P._ambient == self._ambient:
+            return self.retract(x.value)
 
         if x in self.module():
             return self.from_vector(x)

src/sage/categories/finite_dimensional_algebras_with_basis.py

@@ -143,6 +143,17 @@ def radical_basis(self):
                 (B[1] + 6*B[xbar^6], B[xbar] + 6*B[xbar^6], B[xbar^2] + 6*B[xbar^6],
                  B[xbar^3] + 6*B[xbar^6], B[xbar^4] + 6*B[xbar^6], B[xbar^5] + 6*B[xbar^6])
 
+            We compute the radical basis in a subalgebra using
+            the inherited product::
+
+                sage: scoeffs = {('a','e'): {'a':1}, ('b','e'): {'a':1, 'b':1},
+                ....:            ('c','d'): {'a':1}, ('c','e'): {'c':1}}
+                sage: L.<a,b,c,d,e> = LieAlgebra(QQ, scoeffs)
+                sage: MS = MatrixSpace(QQ, 5)
+                sage: A = MS.subalgebra([bg.adjoint_matrix() for bg in L.lie_algebra_generators()])
+                sage: A.radical_basis()
+                (B[1], B[2], B[3], B[4], B[5])
+
             TESTS::
 
                 sage: A = KleinFourGroup().algebra(GF(2))                               # needs sage.groups sage.modules
@@ -162,17 +173,19 @@ def radical_basis(self):
             from sage.matrix.constructor import matrix
             from sage.modules.free_module_element import vector
 
-            product_on_basis = self.product_on_basis
-
             if p == 0:
-                keys = list(self.basis().keys())
-                cache = [{(i,j): c
-                    for i in keys
-                    for j,c in product_on_basis(y,i)}
-                    for y in keys]
-                mat = [ [ sum(x.get((j, i), 0) * c for (i,j),c in y.items())
-                    for x in cache]
-                    for y in cache]
+                B = self.basis()
+                product_on_basis = self.product_on_basis
+                if product_on_basis is NotImplemented:
+                    def product_on_basis(i, j):
+                        return B[i] * B[j]
+
+                keys = B.keys()
+                cache = [{(i, j): c for i in keys for j, c in product_on_basis(y, i)}
+                         for y in keys]
+                mat = [[sum(x.get((j, i), 0) * c for (i,j), c in y.items())
+                        for x in cache]
+                       for y in cache]
 
                 mat = matrix(self.base_ring(), mat)
                 rad_basis = mat.kernel().basis()
@@ -184,24 +197,32 @@ def radical_basis(self):
                 # I imagine that ``pth_root`` would be fastest, but it is not
                 # always available....
                 if hasattr(self.base_ring().one(), 'nth_root'):
-                    root_fcn = lambda s, x : x.nth_root(s)
+                    def root_fcn(s, x):
+                        return x.nth_root(s)
+
                 else:
-                    root_fcn = lambda s, x : x**(1/s)
+                    def root_fcn(s, x):
+                        return x ** (1 / s)
 
-                s, n = 1, self.dimension()
+                s = 1
+                n = self.dimension()
                 B = [b.on_left_matrix() for b in self.basis()]
                 I = B[0].parent().one()
                 while s <= n:
-                    BB = B + [I]
-                    G = matrix([ [(-1)**s * (b*bb).characteristic_polynomial()[n-s]
-                                    for bb in BB] for b in B])
-                    C = G.left_kernel().basis()
+                    # we use that p_{AB}(x) = p_{BA}(x) here
+                    data = [[None]*(len(B)+1) for _ in B]
+                    for i, b in enumerate(B):
+                        for j, bb in enumerate(B[i:], start=i):
+                            val = (-1)**s * (b*bb).charpoly()[n-s]
+                            data[i][j] = data[j][i] = val
+                        data[i][-1] = (-1)**s * b.charpoly()[n-s]
+                    C = matrix(data).left_kernel().basis()
                     if 1 < s < F.order():
                         C = [vector(F, [root_fcn(s, ci) for ci in c]) for c in C]
-                    B = [ sum(ci*b for (ci,b) in zip(c,B)) for c in C ]
+                    B = [sum(ci * b for (ci, b) in zip(c, B)) for c in C]
                     s = p * s
                 e = vector(self.one())
-                rad_basis = [b*e for b in B]
+                rad_basis = [b * e for b in B]
 
             return tuple([self.from_vector(vec) for vec in rad_basis])
 
@@ -260,7 +281,7 @@ def radical(self):
                 sage: # needs sage.graphs sage.modules
                 sage: TestSuite(radical).run()
             """
-            category = AssociativeAlgebras(self.base_ring()).WithBasis().FiniteDimensional().Subobjects()
+            category = AssociativeAlgebras(self.category().base_ring()).WithBasis().FiniteDimensional().Subobjects()
             radical = self.submodule(self.radical_basis(),
                                      category=category,
                                      already_echelonized=True)
@@ -397,7 +418,89 @@ def center(self):
             center.rename("Center of {}".format(self))
             return center
 
-        def principal_ideal(self, a, side='left'):
+        def subalgebra(self, gens, category=None, *args, **opts):
+            r"""
+            Return the subalgebra of ``self`` generated by ``gens``.
+
+            EXAMPLES::
+
+                sage: scoeffs = {('a','e'): {'a':1}, ('b','e'): {'a':1, 'b':1},
+                ....:            ('c','d'): {'a':1}, ('c','e'): {'c':1}}
+                sage: L.<a,b,c,d,e> = LieAlgebra(QQ, scoeffs)
+                sage: MS = MatrixSpace(QQ, 5)
+                sage: A = MS.subalgebra([bg.adjoint_matrix() for bg in L.lie_algebra_generators()])
+                sage: A.dimension()
+                7
+
+                sage: L.<x,y,z> = LieAlgebra(GF(3), {('x','z'): {'x':1, 'y':1}, ('y','z'): {'y':1}})
+                sage: MS = MatrixSpace(L.base_ring(), L.dimension())
+                sage: gens = [b.adjoint_matrix() for b in L.basis()]
+                sage: A = MS.subalgebra(gens)
+                sage: A.dimension()
+                5
+            """
+            # add the unit to make sure it is unital
+            basis = []
+            new_elts = [self(g) for g in gens] + [self.one()]
+            while new_elts:
+                basis = self.echelon_form(basis + new_elts)
+                trailsupp = {b.trailing_support(): b for b in basis}
+                sortsupp = sorted(trailsupp)
+                new_elts = []
+                # We (re)implement the reduction here
+                for b in basis:
+                    for bp in basis:
+                        elt = b * bp
+                        for s in sortsupp:
+                            c = elt[s]
+                            if c:
+                                elt -= c / trailsupp[s].trailing_coefficient() * trailsupp[s]
+                        if elt:
+                            new_elts.append(elt)
+            C = FiniteDimensionalAlgebrasWithBasis(self.category().base_ring())
+            category = C.Subobjects().or_subcategory(category)
+            return self.submodule(basis, check=False, already_echelonized=True,
+                                  category=category)
+
+        def ideal_submodule(self, gens, side='left', category=None, *args, **opts):
+            r"""
+            Return the ``side`` ideal of ``self`` generated by ``gens``
+            as a submodule.
+
+            .. TODO::
+
+                This is not generally compatible with the implementation of
+                the ideals. This method should be folded into the ``ideal``
+                method after the corresponding classes are refactored to
+                be compatible.
+
+            EXAMPLES::
+
+                sage: scoeffs = {('a','e'): {'a':1}, ('b','e'): {'a':1, 'b':1},
+                ....:            ('c','d'): {'a':1}, ('c','e'): {'c':1}}
+                sage: L.<a,b,c,d,e> = LieAlgebra(QQ, scoeffs)
+                sage: MS = MatrixSpace(QQ, 5)
+                sage: I = MS.ideal_submodule([bg.adjoint_matrix() for bg in L.lie_algebra_generators()])
+                sage: I.dimension()
+                25
+            """
+            C = AssociativeAlgebras(self.category().base_ring()).WithBasis().FiniteDimensional()
+            category = C.Subobjects().or_subcategory(category)
+            if gens in self:
+                gens = [self(gens)]
+            if side == 'left':
+                return self.submodule([b * self(g) for b in self.basis() for g in gens],
+                                      category=category, *args, **opts)
+            if side == 'right':
+                return self.submodule([self(g) * b for b in self.basis() for g in gens],
+                                      category=category, *args, **opts)
+            if side == 'twosided':
+                return self.submodule([b * self(g) * bp for b in self.basis()
+                                       for bp in self.basis() for g in gens],
+                                      category=category, *args, **opts)
+            raise ValueError("side must be either 'left', 'right', or 'twosided'")
+
+        def principal_ideal(self, a, side='left', *args, **opts):
             r"""
             Construct the ``side`` principal ideal generated by ``a``.
 
@@ -445,7 +548,7 @@ def principal_ideal(self, a, side='left'):
                 - :meth:`peirce_summand`
             """
             return self.submodule([(a * b if side == 'right' else b * a)
-                                   for b in self.basis()])
+                                   for b in self.basis()], *args, **opts)
 
         @cached_method
         def orthogonal_idempotents_central_mod_radical(self):