From c27983908a6e5349251894fa6e9c716d88eeaa3b Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Wed, 3 Jul 2019 19:36:13 -0400
Subject: [PATCH] eja: add the simple algebra S^n.

---
 mjo/eja/euclidean_jordan_algebra.py | 58 +++++++++++++++++++++++++++++
 1 file changed, 58 insertions(+)

diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py
index 2404af0..f71f75d 100644
--- a/mjo/eja/euclidean_jordan_algebra.py
+++ b/mjo/eja/euclidean_jordan_algebra.py
@@ -596,3 +596,61 @@ def eja_ln(dimension, field=QQ):
     # ambient dimension).
     rank = min(dimension,2)
     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=rank)
+
+
+def eja_sn(dimension, field=QQ):
+    """
+    Return the simple Jordan algebra of ``dimension``-by-``dimension``
+    symmetric matrices over ``field``.
+
+    EXAMPLES::
+
+        sage: J = eja_sn(2)
+        sage: e0, e1, e2 = J.gens()
+        sage: e0*e0
+        e0
+        sage: e1*e1
+        e0 + e2
+        sage: e2*e2
+        e2
+
+    """
+    Qs = []
+
+    # In S^2, for example, we nominally have four coordinates even
+    # though the space is of dimension three only. The vector space V
+    # is supposed to hold the entire long vector, and the subspace W
+    # of V will be spanned by the vectors that arise from symmetric
+    # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
+    V = VectorSpace(field, dimension**2)
+
+    # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
+    # coordinates.
+    S = []
+
+    for i in xrange(dimension):
+        for j in xrange(i+1):
+            Eij = matrix(field, dimension, lambda k,l: k==i and l==j)
+            if i == j:
+                Sij = Eij
+            else:
+                Sij = Eij + Eij.transpose()
+            S.append(Sij)
+
+    def mat2vec(m):
+        return vector(field, m.list())
+
+    W = V.span( mat2vec(s) for s in S )
+
+    for s in S:
+        # Brute force the right-multiplication-by-s matrix by looping
+        # through all elements of the basis and doing the computation
+        # to find out what the corresponding row should be.
+        Q_rows = []
+        for t in S:
+            this_row = mat2vec((s*t + t*s)/2)
+            Q_rows.append(W.coordinates(this_row))
+        Q = matrix(field,Q_rows)
+        Qs.append(Q)
+
+    return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
-- 
2.44.2