Trac #34518: Expand Schubert polynomials in the Key basis

Expand Schubert polynomials in the Key Basis using the peelable tableaux
algorithm of Reiner and Shimozono.

URL: https://trac.sagemath.org/34518
Reported by: tkarn
Ticket author(s): Trevor K. Karn
Reviewer(s): Travis Scrimshaw

Author: Release Manager

src/sage/combinat/diagram.py

@@ -1013,15 +1013,28 @@ def peelable_tableaux(self):
 
         This implementation uses the algorithm suggested in Remark 25
         of [RS1995]_.
+
+        TESTS:
+
+        Corner case::
+
+            sage: from sage.combinat.diagram import NorthwestDiagram
+            sage: D = NorthwestDiagram([])
+            sage: D.peelable_tableaux()
+            {[]}
         """
         # TODO: There is a condition on the first column (if the rows in Dhat
         # are a subset of the rows in the first column) which simplifies the
         # description without performing JDT, so we should implement that
 
+        # empty diagram case
+        if not self:
+            return set([Tableau([])])
+
         # if there is a single column in the diagram then there is only
         # one posslbe peelable tableau.
         if self._n_nonempty_cols == 1:
-            return {Tableau([[i+1] for i, j in self.cells()])}
+            return set([Tableau([[i+1] for i, j in self.cells()])])
 
         first_col = min(j for i, j in self._cells)
 

src/sage/combinat/key_polynomial.py

@@ -479,7 +479,7 @@ def _coerce_map_from_(self, R):
         EXAMPLES::
 
             sage: k = KeyPolynomials(QQ)
-            sage: m1 = k([3,2,4,0]); m1
+            sage: m1 = k([3, 2, 4, 0]); m1
             k[3, 2, 4]
             sage: m2 = k(Composition([3, 2, 4])); m2
             k[3, 2, 4]
@@ -490,10 +490,19 @@ def _coerce_map_from_(self, R):
             sage: z = R.gen()
             sage: z[0] * k([4, 3, 3, 2])
             k[5, 3, 3, 2]
+
+            sage: X = SchubertPolynomialRing(QQ)
+            sage: k(X([4, 3, 2, 1]))
+            k[3, 2, 1]
         """
         P = self._polynomial_ring
         if R is P:
             return self.from_polynomial
+
+        from sage.combinat.schubert_polynomial import SchubertPolynomialRing_xbasis
+        if isinstance(R, SchubertPolynomialRing_xbasis):
+            return self.from_schubert_polynomial
+
         phi = P.coerce_map_from(R)
         if phi is not None:
             return self.coerce_map_from(P) * phi
@@ -645,6 +654,76 @@ def from_polynomial(self, f):
 
         return out
 
+    def from_schubert_polynomial(self, x):
+        r"""
+        Expand a Schubert polynomial in the key basis.
+
+        EXAMPLES::
+
+            sage: k = KeyPolynomials(ZZ)
+            sage: X = SchubertPolynomialRing(ZZ)
+            sage: f = X([2,1,5,4,3])
+            sage: k.from_schubert_polynomial(f)
+            k[1, 0, 2, 1] + k[2, 0, 2] + k[3, 0, 0, 1]
+            sage: k.from_schubert_polynomial(2)
+            2*k[]
+            sage: k(f)
+            k[1, 0, 2, 1] + k[2, 0, 2] + k[3, 0, 0, 1]
+
+            sage: k = KeyPolynomials(GF(7), 4)
+            sage: k.from_schubert_polynomial(f)
+            k[1, 0, 2, 1] + k[2, 0, 2, 0] + k[3, 0, 0, 1]
+
+        TESTS::
+
+            sage: k = KeyPolynomials(ZZ)
+            sage: k.from_schubert_polynomial(k([3,2]))
+            Traceback (most recent call last):
+            ...
+            ValueError: not a Schubert polynomial
+
+            sage: k = KeyPolynomials(ZZ)
+            sage: X = SchubertPolynomialRing(ZZ)
+            sage: it = iter(Compositions())
+            sage: for _ in range(50):
+            ....:     C = next(it)
+            ....:     assert k.from_schubert_polynomial(X(k[C])) == k[C], C
+
+            sage: k = KeyPolynomials(ZZ, 4)
+            sage: X = SchubertPolynomialRing(ZZ)
+            sage: it = iter(k.basis().keys())
+            sage: for _ in range(50):
+            ....:     C = next(it)
+            ....:     assert k.from_schubert_polynomial(X(k[C])) == k[C], C
+        """
+        if x in self.base_ring():
+            return self(x)
+
+        from sage.combinat.schubert_polynomial import SchubertPolynomial_class
+        if not isinstance(x, SchubertPolynomial_class):
+            raise ValueError('not a Schubert polynomial')
+
+        from sage.combinat.diagram import RotheDiagram
+        out = self.zero()
+        if self._k is not None:
+            def build_elt(wt):
+                wt = list(wt)
+                wt += [0] * (self._k - len(wt))
+                return self[wt]
+        else:
+            def build_elt(wt):
+                return self[wt]
+
+        for m, c in x.monomial_coefficients().items():
+            D = RotheDiagram(m)
+            a = self.zero()
+            for d in D.peelable_tableaux():
+                a += build_elt(d.left_key_tableau().weight())
+            out += c * a
+
+        return out
+
+
 def divided_difference(f, i):
     r"""
     Apply the ``i``-th divided difference operator to the polynomial ``f``.

src/sage/combinat/schubert_polynomial.py

@@ -73,15 +73,18 @@
 #
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
-from sage.combinat.free_module import CombinatorialFreeModule
 from sage.categories.all import GradedAlgebrasWithBasis
+from sage.combinat.free_module import CombinatorialFreeModule
+from sage.combinat.key_polynomial import KeyPolynomial
+from sage.combinat.permutation import Permutations, Permutation
+from sage.misc.cachefunc import cached_method
 from sage.rings.integer import Integer
 from sage.rings.integer_ring import ZZ
+from sage.rings.polynomial.infinite_polynomial_element import InfinitePolynomial_sparse
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.rings.polynomial.multi_polynomial import is_MPolynomial
-from sage.combinat.permutation import Permutations, Permutation
+
 import sage.libs.symmetrica.all as symmetrica
-from sage.misc.cachefunc import cached_method
 
 
 def SchubertPolynomialRing(R):
@@ -414,6 +417,17 @@ def _element_constructor_(self, x):
             sage: X(x1^2*x2)
             X[3, 2, 1]
 
+            sage: S.<x> = InfinitePolynomialRing(QQ)
+            sage: X(x[0]^2*x[1])
+            X[3, 2, 1]
+            sage: X(x[0]*x[1]^2*x[2]^2*x[3] + x[0]^2*x[1]^2*x[2]*x[3] + x[0]^2*x[1]*x[2]^2*x[3])
+            X[2, 4, 5, 3, 1]
+
+            sage: from sage.combinat.key_polynomial import KeyPolynomialBasis
+            sage: k = KeyPolynomialBasis(QQ)
+            sage: X(k([3,2,1]))
+            X[4, 3, 2, 1]
+
         TESTS:
 
         We check that :trac:`12924` is fixed::
@@ -429,6 +443,15 @@ def _element_constructor_(self, x):
 
             sage: X([])
             X[1]
+
+        Check the round trip from key polynomials::
+
+            sage: k = KeyPolynomials(ZZ)
+            sage: X = SchubertPolynomialRing(ZZ)
+            sage: it = iter(Permutations())
+            sage: for _ in range(50):
+            ....:     P = next(it)
+            ....:     assert X(k(X(P))) == X(P), P
         """
         if isinstance(x, list):
             # checking the input to avoid symmetrica crashing Sage, see trac 12924
@@ -441,6 +464,14 @@ def _element_constructor_(self, x):
             return self._from_dict({perm: self.base_ring().one()})
         elif is_MPolynomial(x):
             return symmetrica.t_POLYNOM_SCHUBERT(x)
+        elif isinstance(x, InfinitePolynomial_sparse):
+            R = x.polynomial().parent()
+            # massage the term order to be what symmetrica expects
+            S = PolynomialRing(R.base_ring(),
+                               names=list(map(repr, reversed(R.gens()))))
+            return symmetrica.t_POLYNOM_SCHUBERT(S(x.polynomial()))
+        elif isinstance(x, KeyPolynomial):
+            return self(x.expand())
         else:
             raise TypeError
 
@@ -468,3 +499,4 @@ def product_on_basis(self, left, right):
             X[4, 2, 1, 3]
         """
         return symmetrica.mult_schubert_schubert(left, right)
+