eja: factor out the ugly bits of constructing simple matrix EJAs.

diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py
index 7820db524ac2f5858947f48411c00b1d9b4d3eed..8dedc4f68d319123948c36437b68abca2ddb1eef 100644 (file)
--- a/mjo/eja/euclidean_jordan_algebra.py
+++ b/mjo/eja/euclidean_jordan_algebra.py
@@ -620,27 +620,79 @@ def eja_sn(dimension, field=QQ):
         e2
 
     """
-    Qs = []
+    S = _real_symmetric_basis(dimension, field=field)
+    Qs = _multiplication_table_from_matrix_basis(S)
 
-    # 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)
+    return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
+
+
+def random_eja():
+    """
+    Return a "random" finite-dimensional Euclidean Jordan Algebra.
+
+    ALGORITHM:
+
+    For now, we choose a random natural number ``n`` (greater than zero)
+    and then give you back one of the following:
+
+      * The cartesian product of the rational numbers ``n`` times; this is
+        ``QQ^n`` with the Hadamard product.
+
+      * The Jordan spin algebra on ``QQ^n``.
+
+      * The ``n``-by-``n`` rational symmetric matrices with the symmetric
+        product.
+
+    Later this might be extended to return Cartesian products of the
+    EJAs above.
 
+    TESTS::
+
+        sage: random_eja()
+        Euclidean Jordan algebra of degree...
+
+    """
+    n = ZZ.random_element(1,10).abs()
+    constructor = choice([eja_rn, eja_ln, eja_sn])
+    return constructor(dimension=n, field=QQ)
+
+
+
+def _real_symmetric_basis(n, field=QQ):
+    """
+    Return a basis for the space of real symmetric n-by-n matrices.
+    """
     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
     # coordinates.
     S = []
-
-    for i in xrange(dimension):
+    for i in xrange(n):
         for j in xrange(i+1):
-            Eij = matrix(field, dimension, lambda k,l: k==i and l==j)
+            Eij = matrix(field, n, lambda k,l: k==i and l==j)
             if i == j:
                 Sij = Eij
             else:
+                # Beware, orthogonal but not normalized!
                 Sij = Eij + Eij.transpose()
             S.append(Sij)
+    return S
+
+
+def _multiplication_table_from_matrix_basis(basis):
+    """
+    At least three of the five simple Euclidean Jordan algebras have the
+    symmetric multiplication (A,B) |-> (AB + BA)/2, where the
+    multiplication on the right is matrix multiplication. Given a basis
+    for the underlying matrix space, this function returns a
+    multiplication table (obtained by looping through the basis
+    elements) for an algebra of those matrices.
+    """
+    # 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.
+    field = basis[0].base_ring()
+    dimension = basis[0].nrows()
 
     def mat2vec(m):
         return vector(field, m.list())
@@ -648,14 +700,16 @@ def eja_sn(dimension, field=QQ):
     def vec2mat(v):
         return matrix(field, dimension, v.list())
 
-    W = V.span( mat2vec(s) for s in S )
+    V = VectorSpace(field, dimension**2)
+    W = V.span( mat2vec(s) for s in basis )
 
     # Taking the span above reorders our basis (thanks, jerk!) so we
     # need to put our "matrix basis" in the same order as the
     # (reordered) vector basis.
     S = [ vec2mat(b) for b in W.basis() ]
 
-    for s in S:
+    Qs = []
+    for s in basis:
         # Brute force the 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. BEWARE:
@@ -664,13 +718,13 @@ def eja_sn(dimension, field=QQ):
         # constructor uses ROW vectors and not COLUMN vectors. That's
         # why we're computing rows here and not columns.
         Q_rows = []
-        for t in S:
+        for t in basis:
             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)
+    return Qs
 
 
 def random_eja():