]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: add an ungodly hack to get fast charpolys back.
[sage.d.git] / mjo / eja / eja_algebra.py
index 8c66de69ee4a451ce59b0319d1680b8474e7439d..698aa37e2d6698c680ec51d3db1b272e16b6df4d 100644 (file)
@@ -9,14 +9,16 @@ from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
 from sage.combinat.free_module import CombinatorialFreeModule
 from sage.matrix.constructor import matrix
+from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method
 from sage.misc.prandom import choice
 from sage.misc.table import table
-from sage.modules.free_module import VectorSpace
+from sage.modules.free_module import FreeModule, VectorSpace
 from sage.rings.integer_ring import ZZ
-from sage.rings.number_field.number_field import QuadraticField
+from sage.rings.number_field.number_field import NumberField, QuadraticField
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.rings.rational_field import QQ
+from sage.rings.real_lazy import CLF, RLF
 from sage.structure.element import is_Matrix
 
 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
@@ -58,6 +60,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         self._rank = rank
         self._natural_basis = natural_basis
 
+        # TODO: HACK for the charpoly.. needs redesign badly.
+        self._basis_normalizers = None
+
         if category is None:
             category = MagmaticAlgebras(field).FiniteDimensional()
             category = category.WithBasis().Unital()
@@ -134,13 +139,17 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             return self.zero()
 
         natural_basis = self.natural_basis()
-        if elt not in natural_basis[0].matrix_space():
+        basis_space = natural_basis[0].matrix_space()
+        if elt not in basis_space:
             raise ValueError("not a naturally-represented algebra element")
 
-        # Thanks for nothing! Matrix spaces aren't vector
-        # spaces in Sage, so we have to figure out its
-        # natural-basis coordinates ourselves.
-        V = VectorSpace(elt.base_ring(), elt.nrows()*elt.ncols())
+        # Thanks for nothing! Matrix spaces aren't vector spaces in
+        # Sage, so we have to figure out its natural-basis coordinates
+        # ourselves. We use the basis space's ring instead of the
+        # element's ring because the basis space might be an algebraic
+        # closure whereas the base ring of the 3-by-3 identity matrix
+        # could be QQ instead of QQbar.
+        V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols())
         W = V.span_of_basis( _mat2vec(s) for s in natural_basis )
         coords =  W.coordinate_vector(_mat2vec(elt))
         return self.from_vector(coords)
@@ -159,13 +168,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         Ensure that it says what we think it says::
 
             sage: JordanSpinEJA(2, field=QQ)
-            Euclidean Jordan algebra of degree 2 over Rational Field
+            Euclidean Jordan algebra of dimension 2 over Rational Field
             sage: JordanSpinEJA(3, field=RDF)
-            Euclidean Jordan algebra of degree 3 over Real Double Field
+            Euclidean Jordan algebra of dimension 3 over Real Double Field
 
         """
-        # TODO: change this to say "dimension" and fix all the tests.
-        fmt = "Euclidean Jordan algebra of degree {} over {}"
+        fmt = "Euclidean Jordan algebra of dimension {} over {}"
         return fmt.format(self.dimension(), self.base_ring())
 
     def product_on_basis(self, i, j):
@@ -208,12 +216,18 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         determinant).
         """
         z = self._a_regular_element()
-        V = self.vector_space()
-        V1 = V.span_of_basis( (z**k).to_vector() for k in range(self.rank()) )
+        # Don't use the parent vector space directly here in case this
+        # happens to be a subalgebra. In that case, we would be e.g.
+        # two-dimensional but span_of_basis() would expect three
+        # coordinates.
+        V = VectorSpace(self.base_ring(), self.vector_space().dimension())
+        basis = [ (z**k).to_vector() for k in range(self.rank()) ]
+        V1 = V.span_of_basis( basis )
         b =  (V1.basis() + V1.complement().basis())
         return V.span_of_basis(b)
 
 
+
     @cached_method
     def _charpoly_coeff(self, i):
         """
@@ -224,6 +238,19 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         store the trace/determinant (a_{r-1} and a_{0} respectively)
         separate from the entire characteristic polynomial.
         """
+        if self._basis_normalizers is not None:
+             # Must be a matrix class?
+             # WARNING/TODO: this whole mess is mis-designed.
+             n = self.natural_basis_space().nrows()
+             field = self.base_ring().base_ring() # yeeeeaaaahhh
+             J = self.__class__(n, field, False)
+             (_,x,_,_) = J._charpoly_matrix_system()
+             p = J._charpoly_coeff(i)
+             # p might be missing some vars, have to substitute "optionally"
+             pairs = zip(x.base_ring().gens(), self._basis_normalizers)
+             substitutions = { v: v*c for (v,c) in pairs }
+             return p.subs(substitutions)
+
         (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
         R = A_of_x.base_ring()
         if i >= self.rank():
@@ -260,16 +287,40 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         r = self.rank()
         n = self.dimension()
 
-        # Construct a new algebra over a multivariate polynomial ring...
+        # Turn my vector space into a module so that "vectors" can
+        # have multivatiate polynomial entries.
         names = tuple('X' + str(i) for i in range(1,n+1))
         R = PolynomialRing(self.base_ring(), names)
-        # Hack around the fact that our multiplication table is in terms of
-        # algebra elements but the constructor wants it in terms of vectors.
-        vmt = [ tuple(map(lambda x: x.to_vector(), ls))
-                for ls in self._multiplication_table ]
-        J = FiniteDimensionalEuclideanJordanAlgebra(R, tuple(vmt), r)
 
-        idmat = matrix.identity(J.base_ring(), n)
+        # Using change_ring() on the parent's vector space doesn't work
+        # here because, in a subalgebra, that vector space has a basis
+        # and change_ring() tries to bring the basis along with it. And
+        # that doesn't work unless the new ring is a PID, which it usually
+        # won't be.
+        V = FreeModule(R,n)
+
+        # Now let x = (X1,X2,...,Xn) be the vector whose entries are
+        # indeterminates...
+        x = V(names)
+
+        # And figure out the "left multiplication by x" matrix in
+        # that setting.
+        lmbx_cols = []
+        monomial_matrices = [ self.monomial(i).operator().matrix()
+                              for i in range(n) ] # don't recompute these!
+        for k in range(n):
+            ek = self.monomial(k).to_vector()
+            lmbx_cols.append(
+              sum( x[i]*(monomial_matrices[i]*ek)
+                   for i in range(n) ) )
+        Lx = matrix.column(R, lmbx_cols)
+
+        # Now we can compute powers of x "symbolically"
+        x_powers = [self.one().to_vector(), x]
+        for d in range(2, r+1):
+            x_powers.append( Lx*(x_powers[-1]) )
+
+        idmat = matrix.identity(R, n)
 
         W = self._charpoly_basis_space()
         W = W.change_ring(R.fraction_field())
@@ -289,18 +340,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # We want the middle equivalent thing in our matrix, but use
         # the first equivalent thing instead so that we can pass in
         # standard coordinates.
-        x = J.from_vector(W(R.gens()))
-
-        # Handle the zeroth power separately, because computing
-        # the unit element in J is mathematically suspect.
-        x0 = W.coordinate_vector(self.one().to_vector())
-        l1  = [ x0.column() ]
-        l1 += [ W.coordinate_vector((x**k).to_vector()).column()
-                for k in range(1,r) ]
-        l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)]
-        A_of_x = matrix.block(R, 1, n, (l1 + l2))
-        xr = W.coordinate_vector((x**r).to_vector())
-        return (A_of_x, x, xr, A_of_x.det())
+        x_powers = [ W.coordinate_vector(xp) for xp in x_powers ]
+        l2 = [idmat.column(k-1) for k in range(r+1, n+1)]
+        A_of_x = matrix.column(R, n, (x_powers[:r] + l2))
+        return (A_of_x, x, x_powers[r], A_of_x.det())
 
 
     @cached_method
@@ -390,9 +433,32 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
-        if (not x in self) or (not y in self):
-            raise TypeError("arguments must live in this algebra")
-        return x.trace_inner_product(y)
+        X = x.natural_representation()
+        Y = y.natural_representation()
+        return self.__class__.natural_inner_product(X,Y)
+
+
+    def is_trivial(self):
+        """
+        Return whether or not this algebra is trivial.
+
+        A trivial algebra contains only the zero element.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
+
+        EXAMPLES::
+
+            sage: J = ComplexHermitianEJA(3)
+            sage: J.is_trivial()
+            False
+            sage: A = J.zero().subalgebra_generated_by()
+            sage: A.is_trivial()
+            True
+
+        """
+        return self.dimension() == 0
 
 
     def multiplication_table(self):
@@ -454,8 +520,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             Finite family {0: e0, 1: e1, 2: e2}
             sage: J.natural_basis()
             (
-            [1 0]  [0 1]  [0 0]
-            [0 0], [1 0], [0 1]
+            [1 0]  [        0 1/2*sqrt2]  [0 0]
+            [0 0], [1/2*sqrt2         0], [0 1]
             )
 
         ::
@@ -471,11 +537,37 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         """
         if self._natural_basis is None:
-            return tuple( b.to_vector().column() for b in self.basis() )
+            M = self.natural_basis_space()
+            return tuple( M(b.to_vector()) for b in self.basis() )
         else:
             return self._natural_basis
 
 
+    def natural_basis_space(self):
+        """
+        Return the matrix space in which this algebra's natural basis
+        elements live.
+        """
+        if self._natural_basis is None or len(self._natural_basis) == 0:
+            return MatrixSpace(self.base_ring(), self.dimension(), 1)
+        else:
+            return self._natural_basis[0].matrix_space()
+
+
+    @staticmethod
+    def natural_inner_product(X,Y):
+        """
+        Compute the inner product of two naturally-represented elements.
+
+        For example in the real symmetric matrix EJA, this will compute
+        the trace inner-product of two n-by-n symmetric matrices. The
+        default should work for the real cartesian product EJA, the
+        Jordan spin EJA, and the real symmetric matrices. The others
+        will have to be overridden.
+        """
+        return (X.conjugate_transpose()*Y).trace()
+
+
     @cached_method
     def one(self):
         """
@@ -537,6 +629,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return self.linear_combination(zip(self.gens(), coeffs))
 
 
+    def random_element(self):
+        # Temporary workaround for https://trac.sagemath.org/ticket/28327
+        if self.is_trivial():
+            return self.zero()
+        else:
+            s = super(FiniteDimensionalEuclideanJordanAlgebra, self)
+            return s.random_element()
+
+
     def rank(self):
         """
         Return the rank of this EJA.
@@ -609,7 +710,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
             sage: J = RealSymmetricEJA(2)
             sage: J.vector_space()
-            Vector space of dimension 3 over Rational Field
+            Vector space of dimension 3 over...
 
         """
         return self.zero().to_vector().parent().ambient_vector_space()
@@ -650,17 +751,58 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: e2*e2
         e2
 
+    TESTS:
+
+    We can change the generator prefix::
+
+        sage: RealCartesianProductEJA(3, prefix='r').gens()
+        (r0, r1, r2)
+
+    Our inner product satisfies the Jordan axiom::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = RealCartesianProductEJA(n)
+        sage: x = J.random_element()
+        sage: y = J.random_element()
+        sage: z = J.random_element()
+        sage: (x*y).inner_product(z) == y.inner_product(x*z)
+        True
+
     """
-    def __init__(self, n, field=QQ):
+    def __init__(self, n, field=QQ, **kwargs):
         V = VectorSpace(field, n)
         mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
                        for i in range(n) ]
 
         fdeja = super(RealCartesianProductEJA, self)
-        return fdeja.__init__(field, mult_table, rank=n)
+        return fdeja.__init__(field, mult_table, rank=n, **kwargs)
 
     def inner_product(self, x, y):
-        return _usual_ip(x,y)
+        """
+        Faster to reimplement than to use natural representations.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
+
+        TESTS:
+
+        Ensure that this is the usual inner product for the algebras
+        over `R^n`::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: J = RealCartesianProductEJA(n)
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: X = x.natural_representation()
+            sage: Y = y.natural_representation()
+            sage: x.inner_product(y) == J.__class__.natural_inner_product(X,Y)
+            True
+
+        """
+        return x.to_vector().inner_product(y.to_vector())
 
 
 def random_eja():
@@ -696,7 +838,7 @@ def random_eja():
     TESTS::
 
         sage: random_eja()
-        Euclidean Jordan algebra of degree...
+        Euclidean Jordan algebra of dimension...
 
     """
 
@@ -714,9 +856,22 @@ def random_eja():
 
 
 
-def _real_symmetric_basis(n, field=QQ):
+def _real_symmetric_basis(n, field):
     """
     Return a basis for the space of real symmetric n-by-n matrices.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import _real_symmetric_basis
+
+    TESTS::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: B = _real_symmetric_basis(n, QQ)
+        sage: all( M.is_symmetric() for M in  B)
+        True
+
     """
     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
     # coordinates.
@@ -727,16 +882,21 @@ def _real_symmetric_basis(n, field=QQ):
             if i == j:
                 Sij = Eij
             else:
-                # Beware, orthogonal but not normalized!
                 Sij = Eij + Eij.transpose()
             S.append(Sij)
     return tuple(S)
 
 
-def _complex_hermitian_basis(n, field=QQ):
+def _complex_hermitian_basis(n, field):
     """
     Returns a basis for the space of complex Hermitian n-by-n matrices.
 
+    Why do we embed these? Basically, because all of numerical linear
+    algebra assumes that you're working with vectors consisting of `n`
+    entries from a field and scalars from the same field. There's no way
+    to tell SageMath that (for example) the vectors contain complex
+    numbers, while the scalar field is real.
+
     SETUP::
 
         sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
@@ -745,11 +905,15 @@ def _complex_hermitian_basis(n, field=QQ):
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
-        sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
+        sage: field = QuadraticField(2, 'sqrt2')
+        sage: B = _complex_hermitian_basis(n, field)
+        sage: all( M.is_symmetric() for M in  B)
         True
 
     """
-    F = QuadraticField(-1, 'I')
+    R = PolynomialRing(field, 'z')
+    z = R.gen()
+    F = NumberField(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
     I = F.gen()
 
     # This is like the symmetric case, but we need to be careful:
@@ -760,24 +924,33 @@ def _complex_hermitian_basis(n, field=QQ):
     S = []
     for i in xrange(n):
         for j in xrange(i+1):
-            Eij = matrix(field, n, lambda k,l: k==i and l==j)
+            Eij = matrix(F, n, lambda k,l: k==i and l==j)
             if i == j:
                 Sij = _embed_complex_matrix(Eij)
                 S.append(Sij)
             else:
-                # Beware, orthogonal but not normalized! The second one
-                # has a minus because it's conjugated.
+                # The second one has a minus because it's conjugated.
                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
                 S.append(Sij_real)
                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
                 S.append(Sij_imag)
-    return tuple(S)
+
+    # Since we embedded these, we can drop back to the "field" that we
+    # started with instead of the complex extension "F".
+    return tuple( s.change_ring(field) for s in S )
 
 
-def _quaternion_hermitian_basis(n, field=QQ):
+
+def _quaternion_hermitian_basis(n, field, normalize):
     """
     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
 
+    Why do we embed these? Basically, because all of numerical linear
+    algebra assumes that you're working with vectors consisting of `n`
+    entries from a field and scalars from the same field. There's no way
+    to tell SageMath that (for example) the vectors contain complex
+    numbers, while the scalar field is real.
+
     SETUP::
 
         sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
@@ -786,7 +959,8 @@ def _quaternion_hermitian_basis(n, field=QQ):
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
-        sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
+        sage: B = _quaternion_hermitian_basis(n, QQ, False)
+        sage: all( M.is_symmetric() for M in B )
         True
 
     """
@@ -862,7 +1036,7 @@ def _embed_complex_matrix(M):
 
     EXAMPLES::
 
-        sage: F = QuadraticField(-1,'i')
+        sage: F = QuadraticField(-1, 'i')
         sage: x1 = F(4 - 2*i)
         sage: x2 = F(1 + 2*i)
         sage: x3 = F(-i)
@@ -896,8 +1070,8 @@ def _embed_complex_matrix(M):
     field = M.base_ring()
     blocks = []
     for z in M.list():
-        a = z.real()
-        b = z.imag()
+        a = z.vector()[0] # real part, I guess
+        b = z.vector()[1] # imag part, I guess
         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
 
     # We can drop the imaginaries here.
@@ -940,7 +1114,10 @@ def _unembed_complex_matrix(M):
     if not n.mod(2).is_zero():
         raise ValueError("the matrix 'M' must be a complex embedding")
 
-    F = QuadraticField(-1, 'i')
+    field = M.base_ring() # This should already have sqrt2
+    R = PolynomialRing(field, 'z')
+    z = R.gen()
+    F = NumberField(z**2 + 1,'i', embedding=CLF(-1).sqrt())
     i = F.gen()
 
     # Go top-left to bottom-right (reading order), converting every
@@ -1076,10 +1253,6 @@ def _unembed_quaternion_matrix(M):
     return matrix(Q, n/4, elements)
 
 
-# The usual inner product on R^n.
-def _usual_ip(x,y):
-    return x.to_vector().inner_product(y.to_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.
@@ -1106,7 +1279,7 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: e0*e0
         e0
         sage: e1*e1
-        e0 + e2
+        1/2*e0 + 1/2*e2
         sage: e2*e2
         e2
 
@@ -1136,19 +1309,67 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: J(expected) == x*y
         True
 
+    We can change the generator prefix::
+
+        sage: RealSymmetricEJA(3, prefix='q').gens()
+        (q0, q1, q2, q3, q4, q5)
+
+    Our inner product satisfies the Jordan axiom::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = RealSymmetricEJA(n)
+        sage: x = J.random_element()
+        sage: y = J.random_element()
+        sage: z = J.random_element()
+        sage: (x*y).inner_product(z) == y.inner_product(x*z)
+        True
+
+    Our basis is normalized with respect to the natural inner product::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = RealSymmetricEJA(n)
+        sage: all( b.norm() == 1 for b in J.gens() )
+        True
+
+    Left-multiplication operators are symmetric because they satisfy
+    the Jordan axiom::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: x = RealSymmetricEJA(n).random_element()
+        sage: x.operator().matrix().is_symmetric()
+        True
+
     """
-    def __init__(self, n, field=QQ):
-        S = _real_symmetric_basis(n, field=field)
+    def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
+        S = _real_symmetric_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(RealSymmetricEJA, self)
         return fdeja.__init__(field,
                               Qs,
                               rank=n,
-                              natural_basis=S)
+                              natural_basis=S,
+                              **kwargs)
 
-    def inner_product(self, x, y):
-        return _matrix_ip(x,y)
 
 
 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
@@ -1188,28 +1409,74 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: J(expected) == x*y
         True
 
+    We can change the generator prefix::
+
+        sage: ComplexHermitianEJA(2, prefix='z').gens()
+        (z0, z1, z2, z3)
+
+    Our inner product satisfies the Jordan axiom::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = ComplexHermitianEJA(n)
+        sage: x = J.random_element()
+        sage: y = J.random_element()
+        sage: z = J.random_element()
+        sage: (x*y).inner_product(z) == y.inner_product(x*z)
+        True
+
+    Our basis is normalized with respect to the natural inner product::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,4)
+        sage: J = ComplexHermitianEJA(n)
+        sage: all( b.norm() == 1 for b in J.gens() )
+        True
+
+    Left-multiplication operators are symmetric because they satisfy
+    the Jordan axiom::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: x = ComplexHermitianEJA(n).random_element()
+        sage: x.operator().matrix().is_symmetric()
+        True
+
     """
-    def __init__(self, n, field=QQ):
-        S = _complex_hermitian_basis(n)
+    def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
+        S = _complex_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(ComplexHermitianEJA, self)
         return fdeja.__init__(field,
                               Qs,
                               rank=n,
-                              natural_basis=S)
+                              natural_basis=S,
+                              **kwargs)
 
 
-    def inner_product(self, x, y):
-        # 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.
-        return _matrix_ip(x,y)/2
-
+    @staticmethod
+    def natural_inner_product(X,Y):
+        Xu = _unembed_complex_matrix(X)
+        Yu = _unembed_complex_matrix(Y)
+        # 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):
     """
@@ -1248,16 +1515,33 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: J(expected) == x*y
         True
 
+    We can change the generator prefix::
+
+        sage: QuaternionHermitianEJA(2, prefix='a').gens()
+        (a0, a1, a2, a3, a4, a5)
+
+    Our inner product satisfies the Jordan axiom::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = QuaternionHermitianEJA(n)
+        sage: x = J.random_element()
+        sage: y = J.random_element()
+        sage: z = J.random_element()
+        sage: (x*y).inner_product(z) == y.inner_product(x*z)
+        True
+
     """
-    def __init__(self, n, field=QQ):
-        S = _quaternion_hermitian_basis(n)
+    def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
+        S = _quaternion_hermitian_basis(n, field, normalize_basis)
         Qs = _multiplication_table_from_matrix_basis(S)
 
         fdeja = super(QuaternionHermitianEJA, self)
         return fdeja.__init__(field,
                               Qs,
                               rank=n,
-                              natural_basis=S)
+                              natural_basis=S,
+                              **kwargs)
 
     def inner_product(self, x, y):
         # Since a+bi+cj+dk on the diagonal is represented as
@@ -1304,8 +1588,24 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: e2*e3
         0
 
+    We can change the generator prefix::
+
+        sage: JordanSpinEJA(2, prefix='B').gens()
+        (B0, B1)
+
+    Our inner product satisfies the Jordan axiom::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = JordanSpinEJA(n)
+        sage: x = J.random_element()
+        sage: y = J.random_element()
+        sage: z = J.random_element()
+        sage: (x*y).inner_product(z) == y.inner_product(x*z)
+        True
+
     """
-    def __init__(self, n, field=QQ):
+    def __init__(self, n, field=QQ, **kwargs):
         V = VectorSpace(field, n)
         mult_table = [[V.zero() for j in range(n)] for i in range(n)]
         for i in range(n):
@@ -1326,7 +1626,30 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
         # one-dimensional ambient space (because the rank is bounded by
         # the ambient dimension).
         fdeja = super(JordanSpinEJA, self)
-        return fdeja.__init__(field, mult_table, rank=min(n,2))
+        return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
 
     def inner_product(self, x, y):
-        return _usual_ip(x,y)
+        """
+        Faster to reimplement than to use natural representations.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+
+        TESTS:
+
+        Ensure that this is the usual inner product for the algebras
+        over `R^n`::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: J = JordanSpinEJA(n)
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: X = x.natural_representation()
+            sage: Y = y.natural_representation()
+            sage: x.inner_product(y) == J.__class__.natural_inner_product(X,Y)
+            True
+
+        """
+        return x.to_vector().inner_product(y.to_vector())