Allow returning all Casimir elements of a fixed order.

Author: tscrim

src/sage/algebras/lie_algebras/poincare_birkhoff_witt.py

@@ -496,7 +496,7 @@ def degree_on_basis(self, m):
         """
         return m.length()
 
-    def casimir_element(self, order=2):
+    def casimir_element(self, order=2, *args, **kwds):
         r"""
         Return the Casimir element of ``self``.
 
@@ -534,7 +534,7 @@ def casimir_element(self, order=2):
         from sage.rings.infinity import Infinity
         if self._g.dimension() == Infinity:
             raise ValueError("the Lie algebra must be finite dimensional")
-        return self._g.casimir_element(order=order, UEA=self)
+        return self._g.casimir_element(order=order, UEA=self, *args, **kwds)
 
     class Element(CombinatorialFreeModule.Element):
         def _act_on_(self, x, self_on_left):

src/sage/categories/finite_dimensional_lie_algebras_with_basis.py

@@ -1745,10 +1745,10 @@ def universal_commutative_algebra(self):
             R = P[0].parent()
             return R.quotient(P)
 
-        def casimir_element(self, order=2, UEA=None, force_generic=False):
+        def casimir_element(self, order=2, UEA=None, force_generic=False, basis=False):
             r"""
-            Return the Casimir element in the universal enveloping algebra
-            of ``self``.
+            Return a Casimir element of order ``order`` in the universal
+            enveloping algebra of ``self``.
 
             A *Casimir element* of order `k` is a distinguished basis element
             for the center of `U(\mathfrak{g})` of homogeneous degree `k`
@@ -1760,11 +1760,13 @@ def casimir_element(self, order=2, UEA=None, force_generic=False):
             INPUT:
 
             - ``order`` -- (default: ``2``) the order of the Casimir element
-            - ``UEA`` -- (optional) the universal enveloping algebra to
-              return the result in
+            - ``UEA`` -- (optional) the universal enveloping algebra
+              implementation to return the result in
             - ``force_generic`` -- (default: ``False``) if ``True`` for the
               quadratic order, then this uses the default algorithm; otherwise
               this is ignored
+            - ``basis`` -- (default: ``False``) if ``True``, this returns a
+              basis of all Casimir elements of order ``order`` as a list
 
             ALGORITHM:
 
@@ -1832,6 +1834,13 @@ def casimir_element(self, order=2, UEA=None, force_generic=False):
                 sage: L.casimir_element()
                 0
 
+                sage: # needs sage.combinat sage.modules
+                sage: g = LieAlgebra(QQ, cartan_type=['D',2])
+                sage: U = g.pbw_basis()
+                sage: U.casimir_element(2, basis=True)
+                [2*PBW[alpha[2]]*PBW[-alpha[2]] + 1/2*PBW[alphacheck[2]]^2 - PBW[alphacheck[2]],
+                 2*PBW[alpha[1]]*PBW[-alpha[1]] + 1/2*PBW[alphacheck[1]]^2 - PBW[alphacheck[1]]]
+
             TESTS::
 
                 sage: # needs sage.combinat sage.modules
@@ -1856,7 +1865,7 @@ def casimir_element(self, order=2, UEA=None, force_generic=False):
 
             B = self.basis()
 
-            if order == 2 and not force_generic:
+            if order == 2 and not force_generic and not basis:
                 # Special case for the quadratic using the Killing form
                 try:
                     K = self.killing_form_matrix().inverse()
@@ -1896,11 +1905,10 @@ def casimir_element(self, order=2, UEA=None, force_generic=False):
             if ker.dimension() == 0:
                 return self.zero()
 
-            tens = ker.basis()[0]
             del eqns  # no need to hold onto the matrix
 
-            def to_prod(index):
-                coeff = tens[index]
+            def to_prod(vec, index):
+                coeff = vec[index]
                 p = [0] * order
                 base = dim ** (order-1)
                 for i in range(order):
@@ -1910,7 +1918,15 @@ def to_prod(index):
                 p.reverse()
                 return coeff * UEA.prod(UEA(B[keys[i]]) for i in p)
 
-            return UEA.sum(to_prod(index) for index in tens.support())
+            tens = ker.basis()
+
+            if not basis:
+                vec = tens[0]
+                return UEA.sum(to_prod(vec, index) for index in vec.support())
+
+            return [UEA.sum(to_prod(vec, index) for index in vec.support())
+                    for vec in tens]
+
 
     class ElementMethods:
         def adjoint_matrix(self, sparse=False): # In #11111 (more or less) by using matrix of a morphism