unitary DFT for symmetric group algebra

src/sage/combinat/symmetric_group_algebra.py

@@ -2066,8 +2066,74 @@
             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``.

src/sage/combinat/symmetric_group_representations.py

@@ -253,6 +253,8 @@
         return YoungRepresentations_Orthogonal(n, ring=ring, cache_matrices=cache_matrices)
     elif implementation == "specht":
         return SpechtRepresentations(n, ring=ring, cache_matrices=cache_matrices)
+    elif implementation == "unitary":
+        return UnitaryRepresentations(n, ring=ring, cache_matrices=cache_matrices)
     else:
         raise NotImplementedError("only seminormal, orthogonal and specht are implemented")
 
@@ -337,7 +339,7 @@
             sage: spc == SymmetricGroupRepresentation([3])
             True
         """
-        if not isinstance(other, type(other)):
+        if not isinstance(other, type(self)):
             return False
         return (self._ring, self._partition) == (other._ring, other._partition)
 
@@ -659,7 +661,7 @@
                 iu = index_lookup[u]
                 iv = index_lookup[v]
                 M[iu, iu], M[iu, iv], M[iv, iu], M[iv, iv] = \
-                    self._2x2_matrix_entries(beta)
+                    self._2x2_matrix_entries(self._ring(beta))
         return M
 
     def _representation_matrix_uncached(self, permutation):
@@ -914,7 +916,7 @@
         """
         if permutation is None:
             permutation = Permutation(range(1, 1 + self._n))
-        Q = matrix(QQ, len(self._yang_baxter_graph))
+        Q = matrix(self._ring, len(self._yang_baxter_graph))
         for i, v in enumerate(self._dual_vertices):
             for j, u in enumerate(self._yang_baxter_graph):
                 Q[i, j] = self.scalar_product(tuple(permutation.action(v)), u)
@@ -999,6 +1001,170 @@
         """
         return "Specht representations of the symmetric group of order %s! over %s" % (self._n, self._ring)
 
+# #### Unitary Representation ###############################################
+
+
+class UnitaryRepresentation(SymmetricGroupRepresentation_generic_class):
+    def __init__(self, parent, partition):
+        r"""
+        Initialize ``self``.
+
+        EXAMPLES::
+
+            sage: U = SymmetricGroupRepresentation([2,1], "unitary")
+            sage: TestSuite(U).run()
+            sage: U = SymmetricGroupRepresentation([2,1], "unitary", GF(7))
+            Traceback (most recent call last):
+            ...
+            ValueError: the base ring must be a finite field of square order
+            sage: U = SymmetricGroupRepresentation([2,1], "unitary", GF(7**2))
+            sage: TestSuite(U).run()
+        """
+        if parent._ring.characteristic() == 0:
+            orth = SymmetricGroupRepresentation(partition, 'orthogonal')
+            self.representation_matrix = orth.representation_matrix
+            self._representation_matrix_uncached = orth._representation_matrix_uncached
+        else:
+            if not (parent._ring.is_field()
+                    and parent._ring.is_finite()
+                    and parent._ring.order().is_square()):
+                raise ValueError("the base ring must be a finite field of square order")
+            from sage.arith.misc import factorial
+            if parent._ring.characteristic().divides(factorial(parent._n)):
+                raise NotImplementedError("not implemented when p|n!; dimension of invariant forms may be greater than one")
+            self._q = parent._ring.order().sqrt()
+            self._specht = Permutations(sum(partition)).algebra(parent._ring).specht_module(partition)
+        super().__init__(parent, partition)
+
+    def _repr_(self):
+        r"""
+        String representation of ``self``.
+
+        EXAMPLES::
+
+            sage: SymmetricGroupRepresentation([2,1], "unitary")
+            Unitary representation of the symmetric group corresponding to [2, 1]
+        """
+        return f"Unitary representation of the symmetric group corresponding to {self._partition}"
+
+    @cached_method
+    def representation_matrix(self, permutation):
+        r"""
+        Return the matrix representing the ``permutation`` in this
+        irreducible representation.
+
+        .. NOTE::
+
+            This method caches the results.
+
+        EXAMPLES::
+
+            sage: unitary_specht = SymmetricGroupRepresentation([3,1], 'unitary', ring=GF(7**2))
+            sage: unitary_specht.representation_matrix(Permutation([2,1,3,4]))
+            [6 0 0]
+            [0 1 0]
+            [0 0 1]
+            sage: unitary_specht.representation_matrix(Permutation([3,2,1,4]))
+            [       4 2*z2 + 3        0]
+            [5*z2 + 5        3        0]
+            [       0        0        1]
+            sage: all(A * A.H == 1 for A in [unitary_specht.representation_matrix(g)
+            ....:                            for g in Permutations(4)])
+            True
+        """
+        ret = self._representation_matrix_uncached(permutation)
+        ret.set_immutable()
+        return ret
+
+    @lazy_attribute
+    def _unitary_change_basis_matrix(self):
+        """
+        Compute the change of basis matrix.
+
+        EXAMPLES::
+
+            sage: unitary_specht = SymmetricGroupRepresentation([3,1], 'unitary', ring=GF(7**2))
+            sage: unitary_specht._unitary_change_basis_matrix
+            [       1        4        4]
+            [       0 2*z2 + 2 3*z2 + 3]
+            [       0        0 6*z2 + 1]
+            sage: unitary_specht = SymmetricGroupRepresentation([2,2], 'unitary', ring=GF(7**2))
+            sage: unitary_specht._unitary_change_basis_matrix
+            [       1        4]
+            [       0 2*z2 + 2]
+        """
+        from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+        G = Permutations(self._n)
+        F = self._ring
+        rho = self._specht.representation_matrix
+
+        # compute the invariant symmetric bilinear matrix
+        d_rho = self._specht.dimension()
+        R = PolynomialRing(F, 'u', d_rho**2)
+        U_vars = R.gens()
+        Utemp = matrix(R, d_rho, d_rho, U_vars)
+
+        def augmented_matrix(g):
+            rho_g = rho(g)
+            equation_matrix = rho_g.transpose() * Utemp * rho_g.conjugate() - Utemp
+            augmented_system = []
+            for i in range(d_rho):
+                for j in range(d_rho):
+                    linear_expression = equation_matrix[i, j]
+                    row = [linear_expression.coefficient(u) for u in U_vars]
+                    augmented_system.append(row)
+            return augmented_system
+
+        total_system = sum((augmented_matrix(g) for g in G), [])
+        null_space = matrix(F, total_system).right_kernel()
+        U = matrix(F, d_rho, d_rho, null_space.basis()[0])
+        return U.cholesky(extended=True).H
+
+    def _representation_matrix_uncached(self, permutation):
+        r"""
+        Return the matrix representing the ``permutation`` in this
+        irreducible representation.
+
+        .. NOTE::
+
+            This method caches the results.
+
+        EXAMPLES::
+
+            sage: unitary_specht = SymmetricGroupRepresentation([2,2], 'unitary', ring=GF(7**2))
+            sage: unitary_specht._representation_matrix_uncached(Permutation([3,1,4,2]))
+            [       4 5*z2 + 4]
+            [2*z2 + 2        3]
+            sage: unitary_specht._representation_matrix_uncached(Permutation([1,2,4,3]))
+            [6 0]
+            [0 1]
+            sage: all(A * A.H == 1 for A in [unitary_specht._representation_matrix_uncached(g)
+            ....:                            for g in Permutations(4)])
+            True
+        """
+        assert self._ring.characteristic() > 0
+        rho = self._specht.representation_matrix
+        A = self._unitary_change_basis_matrix
+        return A * rho(permutation) * A.inverse()
+
+
+class UnitaryRepresentations(SymmetricGroupRepresentations_class):
+    _default_ring = SR
+
+    Element = UnitaryRepresentation
+
+    def _repr_(self):
+        r"""
+        String representation of ``self``.
+
+        EXAMPLES::
+
+            sage: spc = SymmetricGroupRepresentations(4)
+            sage: spc
+            Specht representations of the symmetric group of order 4! over Integer Ring
+        """
+        return "Unitary representations of the symmetric group of order %s! over %s" % (self._n, self._ring)
+
 # ##### Miscellaneous functions ############################################