Trac #23423: divided_difference in SchubertPolynomialRing should not …

…throw errors on unused variables

{{{
sage: X = SchubertPolynomialRing(ZZ)
sage: a = X([1,2,3,4,5])
sage: a.divided_difference(4)
}}}
This should yield `0`, not an error. Elements of X are polynomials in
infinitely many variables; when they don't use a variable, it doesn't
mean that this variable cannot be divided-differenced over.

Followup to #23403.

URL: https://trac.sagemath.org/23423
Reported by: darij
Ticket author(s): Darij Grinberg
Reviewer(s): Travis Scrimshaw

src/sage/combinat/schubert_polynomial.py

@@ -20,16 +20,17 @@
 from sage.combinat.free_module import CombinatorialFreeModule
 from sage.combinat.combinatorial_algebra import CombinatorialAlgebra
 from sage.categories.all import GradedAlgebrasWithBasis
-from sage.rings.all import Integer, PolynomialRing
+from sage.rings.all import Integer, PolynomialRing, ZZ
 from sage.rings.polynomial.multi_polynomial import is_MPolynomial
-from . import permutation
+from sage.combinat.permutation import Permutations, Permutation
 import sage.libs.symmetrica.all as symmetrica
 
 from sage.combinat.permutation import Permutations
 
 def SchubertPolynomialRing(R):
     """
-    Returns the Schubert polynomial ring over R on the X basis.
+    Returns the Schubert polynomial ring over ``R`` on the X basis
+    (i.e., the basis of the Schubert polynomials).
 
     EXAMPLES::
 
@@ -51,7 +52,8 @@ def SchubertPolynomialRing(R):
 
 def is_SchubertPolynomial(x):
     """
-    Returns True if x is a Schubert polynomial and False otherwise.
+    Returns ``True`` if ``x`` is a Schubert polynomial, and
+    ``False`` otherwise.
 
     EXAMPLES::
 
@@ -69,7 +71,6 @@ def is_SchubertPolynomial(x):
     """
     return isinstance(x, SchubertPolynomial_class)
 
-
 class SchubertPolynomial_class(CombinatorialFreeModule.Element):
     def expand(self):
         """
@@ -108,12 +109,25 @@ def expand(self):
             p = R(p)
         return p
 
-    def divided_difference(self, i):
+    def divided_difference(self, i, algorithm="sage"):
         r"""
         Return the ``i``-th divided difference operator, applied to
         ``self``.
         Here, ``i`` can be either a permutation or a positive integer.
 
+        INPUT:
+        
+        - ``i`` -- permutation or positive integer
+
+        - ``algorithm`` -- (default: ``'sage'``) either ``'sage'``
+          or ``'symmetrica'``; this determines which software is
+          called for the computation
+
+        OUTPUT:
+
+        The result of applying the ``i``-th divided difference
+        operator to ``self``.
+        
         If `i` is a positive integer, then the `i`-th divided
         difference operator `\delta_i` is the linear operator sending
         each polynomial `f = f(x_1, x_2, \ldots, x_n)` (in
@@ -150,17 +164,10 @@ def divided_difference(self, i):
             X[2, 3, 1]
             sage: a.divided_difference([3,2,1])
             X[1]
-
-        If the divided difference operator is not defined for the
-        current number of variables, an error is raised::
-
             sage: a.divided_difference(5)
-            Traceback (most recent call last):
-            ...
-            ValueError: cannot apply \delta_{5} to a (= X[3, 2, 1])
+            0
 
-        For convenience, we define the divided difference of `0`
-        to be `0`::
+        Any divided difference of `0` is `0`::
 
             sage: X.zero().divided_difference(2)
             0
@@ -180,6 +187,19 @@ def divided_difference(self, i):
             X[3, 2, 4, 1]
             sage: a.divided_difference(1).divided_difference(2)
             X[3, 2, 4, 1]
+            sage: a.divided_difference([4,1,3,2])
+            X[1, 4, 2, 3]
+            sage: b = X([4, 1, 3, 2])
+            sage: b.divided_difference(1).divided_difference(2)
+            X[1, 3, 4, 2]
+            sage: b.divided_difference(1).divided_difference(2).divided_difference(3)
+            X[1, 3, 2]
+            sage: b.divided_difference(1).divided_difference(2).divided_difference(3).divided_difference(2)
+            X[1]
+            sage: b.divided_difference(1).divided_difference(2).divided_difference(3).divided_difference(3)
+            0
+            sage: b.divided_difference(1).divided_difference(2).divided_difference(1)
+            0
 
         TESTS:
 
@@ -198,16 +218,64 @@ def divided_difference(self, i):
         """
         if not self: # if self is 0
             return self
-        if isinstance(i, Integer):
-            return symmetrica.divdiff_schubert(i, self)
-        elif i in permutation.Permutations():
-            return symmetrica.divdiff_perm_schubert(i, self)
+        Perms = Permutations()
+        if i in ZZ:
+            if algorithm == "sage":
+                if i <= 0:
+                    raise ValueError(r"cannot apply \delta_{%s} to a (= %s)" % (i, self))
+                # The operator `\delta_i` sends the Schubert
+                # polynomial `X_\pi` (where `\pi` is a finitely supported
+                # permutation of `\{1, 2, 3, \ldots\}`) to:
+                # - the Schubert polynomial X_\sigma`, where `\sigma` is
+                #   obtained from `\pi` by switching the values at `i` and `i+1`,
+                #   if `i` is a descent of `\pi` (that is, `\pi(i) > \pi(i+1)`);
+                # - `0` otherwise.
+                # Notice that distinct `\pi`s lead to distinct `\sigma`s,
+                # so we can use `_from_dict` here.
+                res_dict = {}
+                for pi, coeff in self:
+                    pi = pi[:]
+                    n = len(pi)
+                    if n <= i:
+                        continue
+                    if pi[i-1] < pi[i]:
+                        continue
+                    pi[i-1], pi[i] = pi[i], pi[i-1]
+                    pi = Perms(pi).remove_extra_fixed_points()
+                    res_dict[pi] = coeff
+                return self.parent()._from_dict(res_dict)
+            else: # if algorithm == "symmetrica":
+                return symmetrica.divdiff_schubert(i, self)
+        elif i in Perms:
+            if algorithm == "sage":
+                i = Permutation(i)
+                redw = i.reduced_word()
+                res_dict = {}
+                for pi, coeff in self:
+                    next_pi = False
+                    pi = pi[:]
+                    n = len(pi)
+                    for j in redw:
+                        if n <= j:
+                            next_pi = True
+                            break
+                        if pi[j-1] < pi[j]:
+                            next_pi = True
+                            break
+                        pi[j-1], pi[j] = pi[j], pi[j-1]
+                    if next_pi:
+                        continue
+                    pi = Perms(pi).remove_extra_fixed_points()
+                    res_dict[pi] = coeff
+                return self.parent()._from_dict(res_dict)
+            else: # if algorithm == "symmetrica":
+                return symmetrica.divdiff_perm_schubert(i, self)
         else:
             raise TypeError("i must either be an integer or permutation")
 
     def scalar_product(self, x):
         """
-        Returns the standard scalar product of self and x.
+        Returns the standard scalar product of ``self`` and ``x``.
 
         EXAMPLES::
 
@@ -234,8 +302,8 @@ def scalar_product(self, x):
 
     def multiply_variable(self, i):
         """
-        Returns the Schubert polynomial obtained by multiplying self by the
-        variable `x_i`.
+        Returns the Schubert polynomial obtained by multiplying ``self``
+        by the variable `x_i`.
 
         EXAMPLES::
 
@@ -271,8 +339,8 @@ def __init__(self, R):
         """
         self._name = "Schubert polynomial ring with X basis"
         self._repr_option_bracket = False
-        self._one = permutation.Permutation([1])
-        CombinatorialAlgebra.__init__(self, R, cc = permutation.Permutations(), category = GradedAlgebrasWithBasis(R))
+        self._one = Permutations()([1])
+        CombinatorialAlgebra.__init__(self, R, cc = Permutations(), category = GradedAlgebrasWithBasis(R))
         self.print_options(prefix='X')
 
     def _element_constructor_(self, x):
@@ -305,13 +373,13 @@ def _element_constructor_(self, x):
             #checking the input to avoid symmetrica crashing Sage, see trac 12924
             if not x in Permutations():
                 raise ValueError("The input %s is not a valid permutation"%(x))
-            perm = permutation.Permutation(x).remove_extra_fixed_points()
-            return self._from_dict({ perm: self.base_ring()(1) })
-        elif isinstance(x, permutation.Permutation):
+            perm = Permutation(x).remove_extra_fixed_points()
+            return self._from_dict({ perm: self.base_ring().one() })
+        elif isinstance(x, Permutation):
             if not list(x) in Permutations():
                 raise ValueError("The input %s is not a valid permutation"%(x))
             perm = x.remove_extra_fixed_points()
-            return self._from_dict({ perm: self.base_ring()(1) })
+            return self._from_dict({ perm: self.base_ring().one() })
         elif is_MPolynomial(x):
             return symmetrica.t_POLYNOM_SCHUBERT(x)
         else: