From b558f4e44d9f2c542a7b988f457abccc82a3641c Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Wed, 26 Jan 2022 11:12:20 -0500
Subject: [PATCH] eja: begin work on the ComplexSkewHermitianEJA.

---
 mjo/eja/all.py         |   1 +
 mjo/eja/eja_algebra.py | 126 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 127 insertions(+)

diff --git a/mjo/eja/all.py b/mjo/eja/all.py
index ef020fa..9fe9cb7 100644
--- a/mjo/eja/all.py
+++ b/mjo/eja/all.py
@@ -5,6 +5,7 @@ All imports from mjo.eja modules.
 from mjo.eja.eja_algebra import (AlbertEJA,
                                  BilinearFormEJA,
                                  ComplexHermitianEJA,
+                                 ComplexSkewHermitianEJA,
                                  HadamardEJA,
                                  JordanSpinEJA,
                                  OctonionHermitianEJA,
diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index b76c5fc..fb01646 100644
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -3542,3 +3542,129 @@ def random_eja(max_dimension=None, *args, **kwargs):
         # if the sub-call also Decides on a cartesian product.
         J2 = random_eja(new_max_dimension, *args, **kwargs)
         return cartesian_product([J1,J2])
+
+
+class ComplexSkewHermitianEJA(RationalBasisEJA):
+    r"""
+    The EJA described in Faraut and Koranyi's Exercise III.1.b.
+    """
+    @staticmethod
+    def _denormalized_basis(A):
+        """
+        SETUP::
+
+            sage: from mjo.hurwitz import ComplexMatrixAlgebra
+            sage: from mjo.eja.eja_algebra import ComplexSkewHermitianEJA
+
+        TESTS::
+
+            sage: n = ZZ.random_element(1,2)
+            sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ)
+            sage: B = ComplexSkewHermitianEJA._denormalized_basis(A)
+            sage: all( M.is_skew_hermitian() for M in  B)
+            True
+
+        """
+        es = A.entry_algebra_gens()
+        gen = lambda A,m: A.monomial(m)
+
+        basis = []
+
+        # The size of the blocks. We're going to treat these thing as
+        # 2x2 block matrices,
+        #
+        #   [  x1      x2      ]
+        #   [ -x2^*    x1-conj ]
+        #
+        # where x1 is skew-Hermitian and x2 is symmetric.
+        #
+        m = A.nrows()/2
+
+        # We only loop through the top half of the matrix, because the
+        # bottom can be constructed from the top.
+        for i in range(m):
+
+            # First do the top-left block, which is skew-Hermitian.
+            # We can compute the bottom-right block in the process.
+            for j in range(i+1):
+                if i == j:
+                    # Top-left block's entry.
+                    E_ii = gen(A, (i,j,es[1]))
+
+                    # Bottom-right block's entry.
+                    E_ii += gen(A, (i+m,j+m,es[1])).conjugate()
+                    basis.append(E_ii)
+                else:
+                    for e in es:
+                        # Top-left block's entry.
+                        E_ij  = gen(A, (i,j,e))
+                        E_ij -= E_ij.conjugate_transpose()
+
+                        # Bottom-right block's entry.
+                        F_ij  = gen(A, (i+m,j+m,e)).conjugate()
+                        F_ij -= F_ij.conjugate_transpose()
+
+                        basis.append(E_ij + F_ij)
+
+            # Now do the top-right block, which is symmetric, and compute
+            # the bottom-left block along the way.
+            for j in range(m,i+m+1):
+                if (i+m) == j:
+                    # A symmetric (not Hermitian!) complex matrix can
+                    # have both real and complex entries on its
+                    # diagonal.
+                    for e in es:
+                        # Top-right block's entry.
+                        E_ii = gen(A, (i,j,e))
+
+                        # Bottom-left block's entry.
+                        E_ii -= gen(A, (i-m,j-m,e)).conjugate()
+                        basis.append(E_ii)
+                else:
+                    for e in es:
+                        # Top-right block's entry. BEWARE! We're not
+                        # reflecting across the main diagonal as in
+                        # (i,j)~(j,i). We're only reflecting across
+                        # the diagonal for the top-right block.
+                        E_ij  = gen(A, (i,j,e))
+
+                        # Shift it back to non-offset coords, transpose,
+                        # and put it back:
+                        #
+                        # (i,j) -> (i,j-m) -> (j-m, i) -> (j-m, i+m)
+                        E_ij += gen(A, (j-m,i+m,e))
+
+                        # Bottom-left's block's below-diagonal entry.
+                        # Just shift the top-right coords down m and
+                        # left m.
+                        F_ij  = -gen(A, (i+m,j-m,e)).conjugate()
+                        F_ij += -gen(A, (j,i,e)).conjugate()
+
+                        basis.append(E_ij + F_ij)
+
+        return tuple( basis )
+
+
+    def __init__(self, n, field=AA, **kwargs):
+        # New code; always check the axioms.
+        if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
+
+        from mjo.hurwitz import ComplexMatrixAlgebra
+        A = ComplexMatrixAlgebra(2*n, scalars=field)
+
+        I_n = matrix.identity(ZZ, n)
+        J = matrix.block(ZZ, 2, 2, (0, I_n, -I_n, 0), subdivide=False)
+        J = A.from_list(J.rows())
+
+        def jordan_product(X,Y):
+            return (X*J*Y + Y*J*X)/2
+
+        def inner_product(X,Y):
+            return (X*Y.conjugate_transpose()).trace().real()
+
+        super().__init__(self._denormalized_basis(A),
+                         jordan_product,
+                         inner_product,
+                         field=field,
+                         matrix_space=A,
+                         **kwargs)
-- 
2.44.2