implement unitary DFT of symmetric group over finite fields

unitary DFT of symmetric group for number fields and finite fields

- perform the unitary DFT of the symmetric group over finite fields and number fields
- compute unitary representations of the symmetric group over finite fields
- use real orthogonal representations as unitary representations for symmetric group

Co-Authored-By: Travis Scrimshaw <clfrngrown@aol.com>
Co-Authored-By: Dima Pasechnik <dimpase@gmail.com>

src/sage/combinat/symmetric_group_algebra.py

@@ -247,7 +247,7 @@ def __init__(self, R, W, category):
 
             sage: S = SymmetricGroup(4)
             sage: SGA = S.algebra(QQ)
-            sage: TestSuite(SGA).run(skip="_test_cellular")
+            sage: TestSuite(SGA).run(skip='_test_cellular')
             sage: SGA._test_cellular() # long time
 
         Checking that coercion works between equivalent indexing sets::
@@ -1043,10 +1043,8 @@ def central_orthogonal_idempotents(self):
 
             - :meth:`central_orthogonal_idempotent`
         """
-        out = []
-        for key in sorted(self._blocks_dictionary, reverse=True):
-            out.append(self.central_orthogonal_idempotent(key))
-        return out
+        return [self.central_orthogonal_idempotent(key)
+                for key in sorted(self._blocks_dictionary, reverse=True)]
 
     def central_orthogonal_idempotent(self, la, block=True):
         r"""
@@ -1986,7 +1984,7 @@ def seminormal_basis(self, mult='l2r'):
 
         INPUT:
 
-        - ``mult`` -- string (default: ``'l2r'``). If set to ``'r2l'``,
+        - ``mult`` -- string (default: ``'l2r'``); if set to ``'r2l'``,
           this causes the method to return the list of the
           antipodes (:meth:`antipode`) of all `\epsilon(T, S)`
           instead of the `\epsilon(T, S)` themselves.
@@ -2016,9 +2014,9 @@ def seminormal_basis(self, mult='l2r'):
         basis = []
         for part in Partitions_n(self.n):
             stp = StandardTableaux_shape(part)
-            for t1 in stp:
-                for t2 in stp:
-                    basis.append(self.epsilon_ik(t1, t2, mult=mult))
+            basis.extend(self.epsilon_ik(t1, t2, mult=mult)
+                         for t1 in stp
+                         for t2 in stp)
         return basis
 
     def dft(self, form=None, mult='l2r'):
@@ -2032,7 +2030,7 @@ def dft(self, form=None, mult='l2r'):
 
         INPUT:
 
-        - ``mult`` -- string (default: `l2r`). If set to `r2l`,
+        - ``mult`` -- string (default: `l2r`); if set to `r2l`,
           this causes the method to use the antipodes
           (:meth:`antipode`) of the seminormal basis instead of
           the seminormal basis.
@@ -2068,15 +2066,81 @@ def dft(self, form=None, mult='l2r'):
             return self._dft_seminormal(mult=mult)
         if form == "modular":
             return self._dft_modular()
+        if form == "unitary":
+            return self._dft_unitary()
         raise ValueError("invalid form (= %s)" % form)
 
+    def _dft_unitary(self):
+        """
+        Return the unitary form of the discrete Fourier transform for ``self``.
+
+        EXAMPLES::
+
+            sage: QS3 = SymmetricGroupAlgebra(QQ, 3)
+            sage: QS3._dft_unitary()
+            [-1/6*sqrt3*sqrt2 -1/6*sqrt3*sqrt2 -1/6*sqrt3*sqrt2 -1/6*sqrt3*sqrt2 -1/6*sqrt3*sqrt2 -1/6*sqrt3*sqrt2]
+            [       1/3*sqrt3        1/6*sqrt3       -1/3*sqrt3       -1/6*sqrt3       -1/6*sqrt3        1/6*sqrt3]
+            [               0              1/2                0              1/2             -1/2             -1/2]
+            [               0              1/2                0             -1/2              1/2             -1/2]
+            [       1/3*sqrt3       -1/6*sqrt3        1/3*sqrt3       -1/6*sqrt3       -1/6*sqrt3       -1/6*sqrt3]
+            [-1/6*sqrt3*sqrt2  1/6*sqrt3*sqrt2  1/6*sqrt3*sqrt2 -1/6*sqrt3*sqrt2 -1/6*sqrt3*sqrt2  1/6*sqrt3*sqrt2]
+
+        TESTS::
+
+            sage: QS3 = SymmetricGroupAlgebra(QQ, 3)
+            sage: U = QS3._dft_unitary()
+            sage: U*U.H == identity_matrix(QS3.group().cardinality())
+            True
+            sage: GF5S3 = SymmetricGroupAlgebra(GF(5**2), 3)
+            sage: U = GF5S3._dft_unitary()
+            sage: U*U.H == identity_matrix(GF5S3.group().cardinality())
+            True
+        """
+        from sage.matrix.special import diagonal_matrix
+        F = self.base_ring()
+        G = self.group()
+
+        if F.characteristic() == 0:
+            from sage.misc.functional import sqrt
+            from sage.rings.number_field.number_field import NumberField
+            dft_matrix = self.dft()
+            diag = (dft_matrix * dft_matrix.H).diagonal()
+            primes_needed = {factor for d in diag for factor, _ in d.squarefree_part().factor()}
+            names = [f"sqrt{factor}" for factor in primes_needed]
+            x = PolynomialRing(QQ, 'x').gen()
+            K = NumberField([x**2 - d for d in primes_needed], names=names)
+            sqrt_diag_inv = diagonal_matrix([~sqrt(K(d)) for d in diag])
+            return sqrt_diag_inv * dft_matrix
+
+        # positive characteristic case
+        assert F.characteristic() > 0
+        if not (F.is_field() and F.is_finite() and F.order().is_square()):
+            raise ValueError("the base ring must be a finite field of square order")
+        if F.characteristic().divides(G.cardinality()):
+            raise NotImplementedError("not implemented when p|n!; dimension of invariant forms may be greater than one")
+        q = F.order().sqrt()
+
+        def conj_square_root(u):
+            if not u:
+                return F.zero()
+            z = F.multiplicative_generator()
+            k = u.log(z)
+            if k % (q+1) != 0:
+                raise ValueError(f"Unable to factor as {u} is not in base field GF({q})")
+            return z ** ((k//(q+1)) % (q-1))
+
+        dft_matrix = self.dft()
+        sign_diag = (dft_matrix * dft_matrix.H).diagonal()
+        conj_sqrt_diag_inv = diagonal_matrix([~conj_square_root(d) for d in sign_diag])
+        return conj_sqrt_diag_inv * dft_matrix
+
     def _dft_seminormal(self, mult='l2r'):
         """
         Return the seminormal form of the discrete Fourier transform for ``self``.
 
         INPUT:
 
-        - ``mult`` -- string (default: `l2r`). If set to `r2l`,
+        - ``mult`` -- string (default: `l2r`); if set to `r2l`,
           this causes the method to use the antipodes
           (:meth:`antipode`) of the seminormal basis instead of
           the seminormal basis.
@@ -2133,11 +2197,11 @@ def epsilon_ik(self, itab, ktab, star=0, mult='l2r'):
 
         INPUT:
 
-        - ``itab``, ``ktab`` -- two standard tableaux of size `n`.
+        - ``itab``, ``ktab`` -- two standard tableaux of size `n`
 
-        - ``star`` -- integer (default: `0`).
+        - ``star`` -- integer (default: `0`)
 
-        - ``mult`` -- string (default: `l2r`). If set to `r2l`,
+        - ``mult`` -- string (default: `l2r`); if set to `r2l`,
           this causes the method to return the antipode
           (:meth:`antipode`) of `\epsilon(I, K)` instead of
           `\epsilon(I, K)` itself.
@@ -2507,9 +2571,9 @@ def epsilon_ik(itab, ktab, star=0):
 
     INPUT:
 
-    - ``itab``, ``ktab`` -- two standard tableaux of same size.
+    - ``itab``, ``ktab`` -- two standard tableaux of same size
 
-    - ``star`` -- integer (default: `0`).
+    - ``star`` -- integer (default: `0`)
 
     OUTPUT:
 
@@ -2647,7 +2711,7 @@ def kappa(alpha):
 
     INPUT:
 
-    - ``alpha`` -- integer partition (can be encoded as a list).
+    - ``alpha`` -- integer partition (can be encoded as a list)
 
     OUTPUT:
 
@@ -2684,15 +2748,15 @@ def a(tableau, star=0, base_ring=QQ):
     INPUT:
 
     - ``tableau`` -- Young tableau which contains every integer
-      from `1` to its size precisely once.
+      from `1` to its size precisely once
 
-    - ``star`` -- nonnegative integer (default: `0`). When this
+    - ``star`` -- nonnegative integer (default: `0`); when this
       optional variable is set, the method computes not the row
       projection operator of ``tableau``, but the row projection
       operator of the restriction of ``tableau`` to the entries
       ``1, 2, ..., tableau.size() - star`` instead.
 
-    - ``base_ring`` -- commutative ring (default: ``QQ``). When this
+    - ``base_ring`` -- commutative ring (default: ``QQ``); when this
       optional variable is set, the row projection operator is
       computed over a user-determined base ring instead of `\QQ`.
       (Note that symmetric group algebras currently don't preserve
@@ -2757,7 +2821,7 @@ def b(tableau, star=0, base_ring=QQ):
     INPUT:
 
     - ``tableau`` -- Young tableau which contains every integer
-      from `1` to its size precisely once.
+      from `1` to its size precisely once
 
     - ``star`` -- nonnegative integer (default: `0`). When this
       optional variable is set, the method computes not the column
@@ -3242,7 +3306,7 @@ def _to_sga(self, ind):
         from sage.combinat.rsk import RSK_inverse
         S = ind[1]
         T = ind[2]
-        w = RSK_inverse(T, S, output="permutation")
+        w = RSK_inverse(T, S, output='permutation')
         return self._algebra.kazhdan_lusztig_basis_element(w)
 
 
@@ -3435,7 +3499,7 @@ def __init__(self, R, n, q=None):
         """
         HeckeAlgebraSymmetricGroup_generic.__init__(self, R, n, q)
         self._name += " on the T basis"
-        self.print_options(prefix="T")
+        self.print_options(prefix='T')
 
     def t_action_on_basis(self, perm, i):
         r"""