Fix divided difference

Author: trevorkarn

src/sage/combinat/key_polynomial.py

@@ -26,110 +26,23 @@
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
 
+from sage.categories.graded_algebras_with_basis import GradedAlgebrasWithBasis
 from sage.combinat.free_module import CombinatorialFreeModule
-from sage.combinat.compositions import Compositions
+from sage.combinat.composition import Composition, Compositions
 from sage.misc.inherit_comparison import InheritComparisonClasscallMetaclass
-
-from sage.rings.rational_field import QQ
 from sage.rings.polynomial.infinite_polynomial_ring import InfinitePolynomialRing
+from sage.rings.rational_field import QQ
 
-class KeyPolynomial(CombinatorialFreeModule.Element, 
-                    InheritComparisonClasscallMetaclass):
-    r"""
-    Key polynomials indexed by (weak) compositions.
-    """
-
-    @staticmethod
-    def __classcall_private__(cls, alpha):
-        r"""
-        Normalize the input
-        """
-        R = QQ
-        return KeyPolynomialRing(R)(alpha)
-
-    def pi_operator(self, w):
-        r"""
-        Pi_i or produt
-        """
-        return
-
-    def polynomial_ring(self):
-        return self.parent().polynomial_ring()
-
-    def polynomial_expand(self):
-        R = self.polynomial()
-        out = R.zero()
-
-        for comp, coeff in self.monomial_coefficients():
-            p = [i for i,j in sorted(list(enumerate(comp._list, start=1)), 
-                         key=lambda x: (x[1], x[0]), reverse=True)]
-            P = Permutation(p)  # the minimum sorting permutation
-            m = R.prod(x[i]^e for i,e in enumerate(reversed(P.action(comp))))
-            out += self.pi_operator(w)()
-
-    def atom_expand(self):
-        return
+class KeyPolynomial(CombinatorialFreeModule.Element):
 
-    @cachedmethod
-    def from_polynomial(self):
+    def expand(self):
         r"""
-
+        Return ``self`` written in the monomial basis.
         """
+        out = self.parent().zero()
 
-    def divided_difference(self, i):
-        r"""
-        Return the ``i``-th divided difference operator, applied to ``self``.
-        """
-        if not self:  # if self is 0
-            return self
-        Perms = Permutations()
-        if i in ZZ:
-            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)
-        elif i in Perms:
-            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:
-            raise TypeError("i must either be an integer or permutation")
+        for a, c in self.monomial_coefficients():
+            out += c * self.parent()._pi_i(a)
 
 
 class KeyPolynomialRing(CombinatorialFreeModule):
@@ -149,14 +62,149 @@ def __init__(self, R):
         self._name = "Ring of key polynomials"
         self._repr_option_bracket = False
         CombinatorialFreeModule.__init__(self, R, Compositions(),
-                                         #category=GradedAlgebrasWithBasis(R),
+                                         category=GradedAlgebrasWithBasis(R),
                                          prefix='k')
 
-    def _element_constructor_(self, parent, alpha):
+    def _element_constructor_(self, alpha):
+        #todo: add a case where alpha is an InfinitePolynomialRing element
         return self.element_class(self, alpha)
 
     def polynomial_ring(self):
         r"""
         Return the polynomial ring associated to ``self``.
         """
-        return InfinitePolynomialRing(self.base_ring(), 'x')
+        return InfinitePolynomialRing(self.base_ring(), 'z')
+
+    def from_polynomial(self, f):
+        r"""
+        Expand a polynomial in terms of the key basis.
+        """
+        if not f in self.polynomial_ring():  #todo: change this to accept coerce
+            raise ValueError(f"f must be an element of {self.polynomial_ring()}")
+
+        out = self.zero()
+
+        while f.monomials():
+            m = f.monomials()[0]  #i'm assuming these are sorted in lex order but that is an educated guess.
+            c = f.monomial_coefficient(m)
+            new_term = c * self(Composition(m))
+            f -= new_term.expand()
+            out += new_term
+
+        return out
+
+    def pi(self, w, f):
+        r"""
+
+        """
+        if not hasattr(w, "__iter__"):
+            w = [w]
+        for i in reversed(w):
+            f = self._pi_i(i, f)
+        return f
+
+    def _pi_i(self, i, f):
+        r"""
+        
+        """
+        R = self.polynomial_ring()
+        z, = R.gens()
+        return self._divided_difference(i, z[i] * f)
+
+    def _divided_difference(self, i, f):
+        r"""
+
+        EXAMPLES::
+
+            sage: from sage.combinat.key_polynomial import KeyPolynomialRing
+            sage: K = KeyPolynomialRing(QQ)
+            sage: R = K.polynomial_ring(); R
+            Infinite polynomial ring in x over Rational Field
+            sage: z, = R.gens()
+            sage: f = z[1]*z[2]^3 + z[1]*z[2]*z[3]
+            sage: K._divided_difference(2, f)
+            z[1]*z[2]^2 - z[1]*z[3]^2
+
+        """
+        out = self.polynomial_ring().zero()
+        # linearly extend the divided difference on monomials
+        for m in f.monomials():
+            out += f.monomial_coefficient(m) * self._divided_difference_on_monomial(i, m)
+        return out
+
+    def _divided_difference_on_monomial(self, i, m):
+        r"""
+        Apply the `i`th divided difference operator to `m`.
+
+        NOTE:: `i` can be 0. 
+
+        EXAMPLES::
+
+            sage: from sage.combinat.key_polynomial import KeyPolynomialRing
+            sage: K = KeyPolynomialRing(QQ)
+            sage: R = K.polynomial_ring(); R
+            Infinite polynomial ring in x over Rational Field
+            sage: z, = R.gens()
+            sage: K._divided_difference_on_monomial(1, z[1]*z[2]^3)
+            z_2^2*z_1 + z_2*z_1^2
+            sage: K._divided_difference_on_monomial(2, z[1]*z[2]*z[3])
+            0
+            sage: K._divided_difference_on_monomial(3, z[1]*z[2]*z[3])
+            z[2]*z[1]
+            sage: K._divided_difference_on_monomial(3, z[1]*z[2]*z[4])
+            -z[2]*z[1]
+
+        TODO::
+            add examples where there are multiple factors of xi-x(i+1)
+
+        """
+        # The polynomial ring and generators
+        R = self.polynomial_ring()
+        x, = R.gens()
+
+        # The exponent vector for the monomial
+        exp = list(reversed(m.exponents()[0]))
+
+        if i >= len(exp) - 1:
+            if i == len(exp) - 1:
+                exp = exp + [0]
+            else:
+                # if the transposition acts on two varibles which aren't
+                # present, then the numerator is f - f == 0.
+                return R.zero()
+
+        # Apply the transposition s_i to the exponent vector
+        si_exp = exp[:i] + [exp[i+1], exp[i]] + exp[i+2:]
+
+        # Create the corresponding list of variables in the monomial
+        terms = list(x[i]**j for i,j in enumerate(si_exp) if j)
+
+        # Create the monomial from the list of variables
+        si_m = R.prod(terms)
+
+        # Create the numerator of the divided difference operator
+        f = si_m - m
+
+        # Division using the / operator is not implemented for 
+        # InfinitePolynomialRing, so we want to remove the factor of
+        # x[i+1] - x[i]. If x[i+1] - x[i] is not a factor, it is because
+        # the numerator is zero.
+        try:
+            factors = f.factor()
+            factors_dict = dict(factors)
+        except ArithmeticError:  # if f is zero already
+            return R.zero()
+
+        factors_dict[x[i + 1] - x[i]] = factors_dict[x[i + 1] - x[i]] - 1
+        return R.prod(k**v for k,v in factors_dict.items()) * factors.unit()
+
+def _sorting_permutation(alpha):
+    r"""
+    Get the sorting permutation for a composition ``alpha``.
+
+    The sorting permutation is the permutation of minimal length
+    which turns ``alpha`` into a partition.
+    """
+    key = lambda x: (x[1], x[0])
+    p = [i for i,j in sorted(list(enumerate(a._list, start=1)), key=key, reverse=True)]
+    return Permutation(p)