sagemathgh-39632: Implement matrix Schubert variety ideals
    
<!-- ^ 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 sagemath#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 sagemath#12345". -->

We provide an implementation that returns the defining ideal of the
matrix Schubert variety corresponding to the permutation `w` using the
Fulton generators for the ideal.

### 📝 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 and checked the documentation
preview.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on. For example,
-->
<!-- - sagemath#12345: short description why this is a dependency -->
<!-- - sagemath#34567: ... -->
    
URL: sagemath#39632
Reported by: Travis Scrimshaw
Reviewer(s): Frédéric Chapoton, Travis Scrimshaw

src/doc/en/reference/references/index.rst

@@ -2812,7 +2812,10 @@ REFERENCES:
 .. [Ful1989] \W. Fulton. Algebraic curves: an introduction to algebraic geometry. Addison-Wesley,
              Redwood City CA (1989).
 
-.. [Ful1993] Wiliam Fulton, *Introduction to Toric Varieties*,
+.. [Ful1992] William Fulton. *Flags, Schubert polynomials, degeneracy loci, and
+             determinantal formulas*. Duke Math. J. **65** (1992), no. 3, 381-420.
+
+.. [Ful1993] William Fulton, *Introduction to Toric Varieties*,
             Princeton University Press, 1993.
 
 .. [Ful1997] William Fulton,

src/sage/combinat/diagram.py

@@ -32,6 +32,7 @@
 from sage.combinat.skew_partition import SkewPartition
 from sage.combinat.skew_tableau import SkewTableaux
 from sage.combinat.tableau import Tableau
+from sage.misc.cachefunc import cached_method
 from sage.misc.inherit_comparison import InheritComparisonClasscallMetaclass
 from sage.misc.lazy_import import lazy_import
 from sage.structure.element import Matrix
@@ -534,6 +535,28 @@ def specht_module_dimension(self, base_ring=None):
         from sage.combinat.specht_module import specht_module_rank
         return specht_module_rank(self, base_ring)
 
+    @cached_method
+    def essential_set(self):
+        r"""
+        Return the essential set of ``self`` as defined by Fulton.
+
+        Let `D` be a diagram. Then the *essential set* of `D` are the
+        cells `(i, j) \in D` such that `(i+1, j) \notin D` and
+        `(i, j+1) \notin D`; that is, the maximally southwest elements
+        in each connected component of `D`.
+
+        EXAMPLES::
+
+            sage: w = Permutation([2, 1, 5, 4, 3])
+            sage: D = w.rothe_diagram()
+            sage: D.essential_set()
+            ((0, 0), (2, 3), (3, 2))
+        """
+        ret = [c for c in self._cells if (c[0]+1, c[1]) not in self._cells
+               and (c[0], c[1]+1) not in self._cells]
+        ret.sort()
+        return tuple(ret)
+
 
 class Diagrams(UniqueRepresentation, Parent):
     r"""

src/sage/combinat/permutation.py

@@ -3158,6 +3158,24 @@ def reduced_word_lexmin(self) -> list[int]:
 
         return rw
 
+    def number_of_reduced_words(self):
+        r"""
+        Return the number of reduced words of ``self`` without explicitly
+        computing them all.
+
+        EXAMPLES::
+
+            sage: p = Permutation([6,4,2,5,1,8,3,7])
+            sage: len(p.reduced_words()) == p.number_of_reduced_words()                 # needs sage.combinat
+            True
+        """
+        Tx = self.rothe_diagram().peelable_tableaux()
+        return sum(map(_tableau_contribution, Tx))
+
+    ##################
+    # Rothe diagrams #
+    ##################
+
     def rothe_diagram(self):
         r"""
         Return the Rothe diagram of ``self``.
@@ -3176,20 +3194,100 @@ def rothe_diagram(self):
         from sage.combinat.diagram import RotheDiagram
         return RotheDiagram(self)
 
-    def number_of_reduced_words(self):
+    def rank_matrix(self):
         r"""
-        Return the number of reduced words of ``self`` without explicitly
-        computing them all.
+        Return the rank matrix of ``self``.
+
+        Let `P = [P_{ij}]_{i, j=1}^n` be the permutation matrix of `w \in S_n`.
+        The rank matrix is the `n \times n` matrix with entries
+
+        .. MATH::
+
+            r_w(i, j) = \sum_{a=1}^i \sum_{b=1}^j P_{ij}.
 
         EXAMPLES::
 
-            sage: p = Permutation([6,4,2,5,1,8,3,7])
-            sage: len(p.reduced_words()) == p.number_of_reduced_words()                 # needs sage.combinat
-            True
+            sage: w = Permutation([2, 1, 5, 4, 3])
+            sage: w.to_matrix()
+            [0 1 0 0 0]
+            [1 0 0 0 0]
+            [0 0 0 0 1]
+            [0 0 0 1 0]
+            [0 0 1 0 0]
+            sage: w.rank_matrix()
+            [0 1 1 1 1]
+            [1 2 2 2 2]
+            [1 2 2 2 3]
+            [1 2 2 3 4]
+            [1 2 3 4 5]
         """
-        Tx = self.rothe_diagram().peelable_tableaux()
+        ret = self.to_matrix()
+        n = ret.nrows()
+        for j in range(1, n):
+            ret[0, j] += ret[0, j-1]
+        for i in range(1, n):
+            ret[i, 0] += ret[i-1, 0]
+            for j in range(1, n):
+                # Compute by inclusion-exclusion
+                ret[i, j] += ret[i-1, j] + ret[i, j-1] - ret[i-1, j-1]
+        return ret
 
-        return sum(map(_tableau_contribution, Tx))
+    def schubert_determinant_ideal(self):
+        r"""
+        Return the Schubert determinant ideal of ``self``.
+
+        From Lemma 3.10 [Ful1992]_, the matrix Schubert variety of
+        a permutation `w \in S_n` is defined by the ideal
+
+        .. MATH::
+
+            I_w = \sum_{(i, j)} I_{r_w(i, j)+1}(i, j),
+
+        where the sum is over the :meth:`essential set
+        <sage.combinat.diagram.Diagram.essential_set>`, `[r_w(i, j)]_{i, j}`
+        is the :meth:`rank_matrix` of `w`, and `I_k(i, j)` is the ideal
+        generated by the `k \times k` minors of the `i \times j`
+        matrix `[z_{ab}]_{a, b}`.
+
+        These ideals are known to be prime of codimension equal to `\ell(w)`
+        (the length of `w`) such that `R / I_w` is Cohen-Macaulay, where `R`
+        is the polynomial ring `\QQ[z_{ab} | 1 \leq a, b \leq n]` from
+        Prop. 3.10 of [Ful1992]_.
+
+        EXAMPLES::
+
+            sage: w = Permutation([3, 1, 4, 2])
+            sage: Iw = w.schubert_determinant_ideal()
+            sage: Iw.gens()
+            [z00, z01, -z01*z10 + z00*z11, -z01*z20 + z00*z21, -z11*z20 + z10*z21]
+            sage: Iw.dimension()
+            13
+            sage: Iw.dimension() + w.length()
+            16
+            sage: Iw.is_prime()
+            True
+
+            sage: w = Permutation([2, 1, 5, 4, 3])
+            sage: Iw = w.schubert_determinant_ideal()
+            sage: Iw.dimension()
+            21
+            sage: w.length() + Iw.dimension()
+            25
+            sage: Iw.is_prime()
+            True
+        """
+        from sage.rings.rational_field import QQ
+        from sage.matrix.constructor import matrix
+        n = len(self)
+        PR = PolynomialRing(QQ, n, var_array='z')
+        z = PR.gens()
+        Z = matrix(PR, [[z[r*n+c] for c in range(n)] for r in range(n)])
+        rk = self.rank_matrix()
+        gens = []
+        for i, j in self.rothe_diagram().essential_set():
+            # we apply the transpose to the rank matrix to match conventions
+            gens.extend(Z.submatrix(0, 0, i+1, j+1).minors(rk[j, i] + 1))
+        return PR.ideal(gens)
 
     ################
     # Fixed Points #