eja: add matrix inner products, and factor out the two used in R^n.

diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py
index ca344edc709352dea2d1940a7646e058066d4b8d..a73cfaf689c7206f7644e10021006e4bcf6f08f7 100644 (file)
--- a/mjo/eja/euclidean_jordan_algebra.py
+++ b/mjo/eja/euclidean_jordan_algebra.py
@@ -244,6 +244,24 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 sage: J(x).inner_product(J(y))
                 32
 
+            The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
+            multiplication is the usual matrix multiplication in `S^n`,
+            so the inner product of the identity matrix with itself
+            should be the `n`::
+
+                sage: J = RealSymmetricSimpleEJA(3)
+                sage: J.one().inner_product(J.one())
+                3
+
+            Likewise, the inner product on `C^n` is `<X,Y> =
+            Re(trace(X*Y))`, where we must necessarily take the real
+            part because the product of Hermitian matrices may not be
+            Hermitian::
+
+                sage: J = ComplexHermitianSimpleEJA(3)
+                sage: J.one().inner_product(J.one())
+                3
+
             TESTS:
 
             Ensure that we can always compute an inner product, and that
@@ -941,13 +959,10 @@ def eja_rn(dimension, field=QQ):
     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
            for i in xrange(dimension) ]
 
-    # The usual inner product on R^n.
-    ip = lambda x, y: x.vector().inner_product(y.vector())
-
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=dimension,
-                                                   inner_product=ip)
+                                                   inner_product=_usual_ip)
 
 
 
@@ -1161,14 +1176,26 @@ def _unembed_complex_matrix(M):
         for j in xrange(n/2):
             submat = M[2*k:2*k+2,2*j:2*j+2]
             if submat[0,0] != submat[1,1]:
-                raise ArgumentError('bad real submatrix')
+                raise ValueError('bad real submatrix')
             if submat[0,1] != -submat[1,0]:
-                raise ArgumentError('bad imag submatrix')
+                raise ValueError('bad imag submatrix')
             z = submat[0,0] + submat[1,0]*i
             elements.append(z)
 
     return matrix(F, n/2, elements)
 
+# The usual inner product on R^n.
+def _usual_ip(x,y):
+    return x.vector().inner_product(y.vector())
+
+# The inner product used for the real symmetric simple EJA.
+# We keep it as a separate function because e.g. the complex
+# algebra uses the same inner product, except divided by 2.
+def _matrix_ip(X,Y):
+    X_mat = X.natural_representation()
+    Y_mat = Y.natural_representation()
+    return (X_mat*Y_mat).trace()
+
 
 def RealSymmetricSimpleEJA(n, field=QQ):
     """
@@ -1204,7 +1231,8 @@ def RealSymmetricSimpleEJA(n, field=QQ):
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=n,
-                                                   natural_basis=T)
+                                                   natural_basis=T,
+                                                   inner_product=_matrix_ip)
 
 
 def ComplexHermitianSimpleEJA(n, field=QQ):
@@ -1227,10 +1255,21 @@ def ComplexHermitianSimpleEJA(n, field=QQ):
     """
     S = _complex_hermitian_basis(n)
     (Qs, T) = _multiplication_table_from_matrix_basis(S)
+
+    # Since a+bi on the diagonal is represented as
+    #
+    #   a + bi  = [  a  b  ]
+    #             [ -b  a  ],
+    #
+    # we'll double-count the "a" entries if we take the trace of
+    # the embedding.
+    ip = lambda X,Y: _matrix_ip(X,Y)/2
+
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=n,
-                                                   natural_basis=T)
+                                                   natural_basis=T,
+                                                   inner_product=ip)
 
 
 def QuaternionHermitianSimpleEJA(n):
@@ -1277,6 +1316,13 @@ def JordanSpinSimpleEJA(n, field=QQ):
         sage: e2*e3
         0
 
+    In one dimension, this is the reals under multiplication::
+
+        sage: J1 = JordanSpinSimpleEJA(1)
+        sage: J2 = eja_rn(1)
+        sage: J1 == J2
+        True
+
     """
     Qs = []
     id_matrix = identity_matrix(field, n)
@@ -1291,13 +1337,10 @@ def JordanSpinSimpleEJA(n, field=QQ):
         Qi[0,0] = Qi[0,0] * ~field(2)
         Qs.append(Qi)
 
-    # The usual inner product on R^n.
-    ip = lambda x, y: x.vector().inner_product(y.vector())
-
     # The rank of the spin factor algebra is two, UNLESS we're in a
     # one-dimensional ambient space (the rank is bounded by the
     # ambient dimension).
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=min(n,2),
-                                                   inner_product=ip)
+                                                   inner_product=_usual_ip)