eja: enable normalization of the natural quaternion basis.

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 2bfa3714420eff322b6b2fb19e177107a191ef29..b439ff9a46a578d22c9444d96862a69ea5c4a95d 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -1232,7 +1232,10 @@ def _unembed_quaternion_matrix(M):
     if not n.mod(4).is_zero():
         raise ValueError("the matrix 'M' must be a complex embedding")
 
-    Q = QuaternionAlgebra(QQ,-1,-1)
+    # Use the base ring of the matrix to ensure that its entries can be
+    # multiplied by elements of the quaternion algebra.
+    field = M.base_ring()
+    Q = QuaternionAlgebra(field,-1,-1)
     i,j,k = Q.gens()
 
     # Go top-left to bottom-right (reading order), converting every
@@ -1246,8 +1249,10 @@ def _unembed_quaternion_matrix(M):
                 raise ValueError('bad on-diagonal submatrix')
             if submat[0,1] != -submat[1,0].conjugate():
                 raise ValueError('bad off-diagonal submatrix')
-            z  = submat[0,0].real() + submat[0,0].imag()*i
-            z += submat[0,1].real()*j + submat[0,1].imag()*k
+            z  = submat[0,0].vector()[0]   # real part
+            z += submat[0,0].vector()[1]*i # imag part
+            z += submat[0,1].vector()[0]*j # real part
+            z += submat[0,1].vector()[1]*k # imag part
             elements.append(z)
 
     return matrix(Q, n/4, elements)
@@ -1486,6 +1491,8 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
         # The trace need not be real; consider Xu = (i*I) and Yu = I.
         return ((Xu*Yu).trace()).vector()[0] # real part, I guess
 
+
+
 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
@@ -1502,7 +1509,7 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     The dimension of this algebra is `n^2`::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
+        sage: n = ZZ.random_element(1,4)
         sage: J = QuaternionHermitianEJA(n)
         sage: J.dimension() == 2*(n^2) - n
         True
@@ -1510,7 +1517,7 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     The Jordan multiplication is what we think it is::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
+        sage: n = ZZ.random_element(1,4)
         sage: J = QuaternionHermitianEJA(n)
         sage: x = J.random_element()
         sage: y = J.random_element()
@@ -1531,7 +1538,7 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     Our inner product satisfies the Jordan axiom::
 
         sage: set_random_seed()
-        sage: n = ZZ.random_element(1,5)
+        sage: n = ZZ.random_element(1,4)
         sage: J = QuaternionHermitianEJA(n)
         sage: x = J.random_element()
         sage: y = J.random_element()
@@ -1539,9 +1546,46 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: (x*y).inner_product(z) == y.inner_product(x*z)
         True
 
+    Our natural basis is normalized with respect to the natural inner
+    product unless we specify otherwise::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,4)
+        sage: J = QuaternionHermitianEJA(n)
+        sage: all( b.norm() == 1 for b in J.gens() )
+        True
+
+    Since our natural basis is normalized with respect to the natural
+    inner product, and since we know that this algebra is an EJA, any
+    left-multiplication operator's matrix will be symmetric because
+    natural->EJA basis representation is an isometry and within the EJA
+    the operator is self-adjoint by the Jordan axiom::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: x = QuaternionHermitianEJA(n).random_element()
+        sage: x.operator().matrix().is_symmetric()
+        True
+
     """
     def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
-        S = _quaternion_hermitian_basis(n, field, normalize_basis)
+        S = _quaternion_hermitian_basis(n, field)
+
+        if n > 1 and normalize_basis:
+            # We'll need sqrt(2) to normalize the basis, and this
+            # winds up in the multiplication table, so the whole
+            # algebra needs to be over the field extension.
+            R = PolynomialRing(field, 'z')
+            z = R.gen()
+            p = z**2 - 2
+            if p.is_irreducible():
+                field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt())
+                S = [ s.change_ring(field) for s in S ]
+            self._basis_normalizers = tuple(
+                ~(self.__class__.natural_inner_product(s,s).sqrt())
+                for s in S )
+            S = tuple( s*c for (s,c) in zip(S,self._basis_normalizers) )
+
         Qs = _multiplication_table_from_matrix_basis(S)
 
         fdeja = super(QuaternionHermitianEJA, self)
@@ -1551,17 +1595,16 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
                               natural_basis=S,
                               **kwargs)
 
-    def inner_product(self, x, y):
-        # Since a+bi+cj+dk on the diagonal is represented as
-        #
-        #   a + bi +cj + dk = [  a  b  c  d]
-        #                     [ -b  a -d  c]
-        #                     [ -c  d  a -b]
-        #                     [ -d -c  b  a],
-        #
-        # we'll quadruple-count the "a" entries if we take the trace of
-        # the embedding.
-        return _matrix_ip(x,y)/4
+    @staticmethod
+    def natural_inner_product(X,Y):
+        Xu = _unembed_quaternion_matrix(X)
+        Yu = _unembed_quaternion_matrix(Y)
+        # The trace need not be real; consider Xu = (i*I) and Yu = I.
+        # The result will be a quaternion algebra element, which doesn't
+        # have a "vector" method, but does have coefficient_tuple() method
+        # that returns the coefficients of 1, i, j, and k -- in that order.
+        return ((Xu*Yu).trace()).coefficient_tuple()[0]
+
 
 
 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):