+ sage: J = QuaternionHermitianSimpleEJA(3)
+ sage: J.one()
+ e0 + e9 + e14
+ sage: J.one().natural_representation()
+ [1 0 0 0 0 0 0 0 0 0 0 0]
+ [0 1 0 0 0 0 0 0 0 0 0 0]
+ [0 0 1 0 0 0 0 0 0 0 0 0]
+ [0 0 0 1 0 0 0 0 0 0 0 0]
+ [0 0 0 0 1 0 0 0 0 0 0 0]
+ [0 0 0 0 0 1 0 0 0 0 0 0]
+ [0 0 0 0 0 0 1 0 0 0 0 0]
+ [0 0 0 0 0 0 0 1 0 0 0 0]
+ [0 0 0 0 0 0 0 0 1 0 0 0]
+ [0 0 0 0 0 0 0 0 0 1 0 0]
+ [0 0 0 0 0 0 0 0 0 0 1 0]
+ [0 0 0 0 0 0 0 0 0 0 0 1]
+
+ """
+ B = self.parent().natural_basis()
+ W = B[0].matrix_space()
+ return W.linear_combination(zip(self.vector(), B))
+
+
+ def operator_matrix(self):
+ """
+ Return the matrix that represents left- (or right-)
+ multiplication by this element in the parent algebra.
+
+ We have to override this because the superclass method
+ returns a matrix that acts on row vectors (that is, on
+ the right).
+
+ EXAMPLES:
+
+ Test the first polarization identity from my notes, Koecher Chapter
+ III, or from Baes (2.3)::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: x = J.random_element()
+ sage: y = J.random_element()
+ sage: Lx = x.operator_matrix()
+ sage: Ly = y.operator_matrix()
+ sage: Lxx = (x*x).operator_matrix()
+ sage: Lxy = (x*y).operator_matrix()
+ sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
+ True
+
+ Test the second polarization identity from my notes or from
+ Baes (2.4)::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: x = J.random_element()
+ sage: y = J.random_element()
+ sage: z = J.random_element()
+ sage: Lx = x.operator_matrix()
+ sage: Ly = y.operator_matrix()
+ sage: Lz = z.operator_matrix()
+ sage: Lzy = (z*y).operator_matrix()
+ sage: Lxy = (x*y).operator_matrix()
+ sage: Lxz = (x*z).operator_matrix()
+ sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
+ True
+
+ Test the third polarization identity from my notes or from
+ Baes (2.5)::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: u = J.random_element()
+ sage: y = J.random_element()
+ sage: z = J.random_element()
+ sage: Lu = u.operator_matrix()
+ sage: Ly = y.operator_matrix()
+ sage: Lz = z.operator_matrix()
+ sage: Lzy = (z*y).operator_matrix()
+ sage: Luy = (u*y).operator_matrix()
+ sage: Luz = (u*z).operator_matrix()
+ sage: Luyz = (u*(y*z)).operator_matrix()
+ sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
+ sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
+ sage: bool(lhs == rhs)
+ True
+
+ """
+ fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
+ return fda_elt.matrix().transpose()
+
+
+ def quadratic_representation(self, other=None):
+ """
+ Return the quadratic representation of this element.
+
+ EXAMPLES:
+
+ The explicit form in the spin factor algebra is given by
+ Alizadeh's Example 11.12::
+
+ sage: set_random_seed()
+ sage: n = ZZ.random_element(1,10)
+ sage: J = JordanSpinSimpleEJA(n)
+ sage: x = J.random_element()
+ sage: x_vec = x.vector()
+ sage: x0 = x_vec[0]
+ sage: x_bar = x_vec[1:]
+ sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
+ sage: B = 2*x0*x_bar.row()
+ sage: C = 2*x0*x_bar.column()
+ sage: D = identity_matrix(QQ, n-1)
+ sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
+ sage: D = D + 2*x_bar.tensor_product(x_bar)
+ sage: Q = block_matrix(2,2,[A,B,C,D])
+ sage: Q == x.quadratic_representation()
+ True
+
+ Test all of the properties from Theorem 11.2 in Alizadeh::
+
+ sage: set_random_seed()
+ sage: J = random_eja()
+ sage: x = J.random_element()
+ sage: y = J.random_element()
+
+ Property 1:
+
+ sage: actual = x.quadratic_representation(y)
+ sage: expected = ( (x+y).quadratic_representation()
+ ....: -x.quadratic_representation()
+ ....: -y.quadratic_representation() ) / 2
+ sage: actual == expected
+ True
+
+ Property 2:
+
+ sage: alpha = QQ.random_element()
+ sage: actual = (alpha*x).quadratic_representation()
+ sage: expected = (alpha^2)*x.quadratic_representation()
+ sage: actual == expected
+ True
+
+ Property 5:
+
+ sage: Qy = y.quadratic_representation()
+ sage: actual = J(Qy*x.vector()).quadratic_representation()
+ sage: expected = Qy*x.quadratic_representation()*Qy
+ sage: actual == expected
+ True
+
+ Property 6:
+
+ sage: k = ZZ.random_element(1,10)
+ sage: actual = (x^k).quadratic_representation()
+ sage: expected = (x.quadratic_representation())^k
+ sage: actual == expected
+ True
+
+ """
+ if other is None:
+ other=self
+ elif not other in self.parent():
+ raise TypeError("'other' must live in the same algebra")
+
+ L = self.operator_matrix()
+ M = other.operator_matrix()
+ return ( L*M + M*L - (self*other).operator_matrix() )
+
+
+ def span_of_powers(self):
+ """
+ Return the vector space spanned by successive powers of
+ this element.
+ """
+ # The dimension of the subalgebra can't be greater than
+ # the big algebra, so just put everything into a list
+ # and let span() get rid of the excess.
+ V = self.vector().parent()
+ return V.span( (self**d).vector() for d in xrange(V.dimension()) )
+
+
+ def subalgebra_generated_by(self):
+ """
+ Return the associative subalgebra of the parent EJA generated
+ by this element.
+
+ TESTS::
+
+ sage: set_random_seed()
+ sage: x = random_eja().random_element()
+ sage: x.subalgebra_generated_by().is_associative()
+ True
+
+ Squaring in the subalgebra should be the same thing as
+ squaring in the superalgebra::
+
+ sage: set_random_seed()
+ sage: x = random_eja().random_element()
+ sage: u = x.subalgebra_generated_by().random_element()
+ sage: u.operator_matrix()*u.vector() == (u**2).vector()
+ True
+
+ """
+ # First get the subspace spanned by the powers of myself...
+ V = self.span_of_powers()
+ F = self.base_ring()
+
+ # Now figure out the entries of the right-multiplication
+ # matrix for the successive basis elements b0, b1,... of
+ # that subspace.
+ mats = []
+ for b_right in V.basis():
+ eja_b_right = self.parent()(b_right)
+ b_right_rows = []
+ # The first row of the right-multiplication matrix by
+ # b1 is what we get if we apply that matrix to b1. The
+ # second row of the right multiplication matrix by b1
+ # is what we get when we apply that matrix to b2...
+ #
+ # IMPORTANT: this assumes that all vectors are COLUMN
+ # vectors, unlike our superclass (which uses row vectors).
+ for b_left in V.basis():
+ eja_b_left = self.parent()(b_left)
+ # Multiply in the original EJA, but then get the
+ # coordinates from the subalgebra in terms of its
+ # basis.
+ this_row = V.coordinates((eja_b_left*eja_b_right).vector())
+ b_right_rows.append(this_row)
+ b_right_matrix = matrix(F, b_right_rows)
+ mats.append(b_right_matrix)
+
+ # It's an algebra of polynomials in one element, and EJAs
+ # are power-associative.
+ #
+ # TODO: choose generator names intelligently.
+ return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
+
+
+ def subalgebra_idempotent(self):
+ """
+ Find an idempotent in the associative subalgebra I generate
+ using Proposition 2.3.5 in Baes.
+
+ TESTS::
+
+ sage: set_random_seed()
+ sage: J = eja_rn(5)
+ sage: c = J.random_element().subalgebra_idempotent()
+ sage: c^2 == c
+ True
+ sage: J = JordanSpinSimpleEJA(5)
+ sage: c = J.random_element().subalgebra_idempotent()
+ sage: c^2 == c
+ True
+
+ """
+ if self.is_nilpotent():
+ raise ValueError("this only works with non-nilpotent elements!")
+
+ V = self.span_of_powers()
+ J = self.subalgebra_generated_by()
+ # Mis-design warning: the basis used for span_of_powers()
+ # and subalgebra_generated_by() must be the same, and in
+ # the same order!
+ u = J(V.coordinates(self.vector()))
+
+ # The image of the matrix of left-u^m-multiplication
+ # will be minimal for some natural number s...
+ s = 0
+ minimal_dim = V.dimension()
+ for i in xrange(1, V.dimension()):
+ this_dim = (u**i).operator_matrix().image().dimension()
+ if this_dim < minimal_dim:
+ minimal_dim = this_dim
+ s = i
+
+ # Now minimal_matrix should correspond to the smallest
+ # non-zero subspace in Baes's (or really, Koecher's)
+ # proposition.
+ #
+ # However, we need to restrict the matrix to work on the
+ # subspace... or do we? Can't we just solve, knowing that
+ # A(c) = u^(s+1) should have a solution in the big space,
+ # too?
+ #
+ # Beware, solve_right() means that we're using COLUMN vectors.
+ # Our FiniteDimensionalAlgebraElement superclass uses rows.
+ u_next = u**(s+1)
+ A = u_next.operator_matrix()
+ c_coordinates = A.solve_right(u_next.vector())
+
+ # Now c_coordinates is the idempotent we want, but it's in
+ # the coordinate system of the subalgebra.
+ #
+ # We need the basis for J, but as elements of the parent algebra.
+ #
+ basis = [self.parent(v) for v in V.basis()]
+ return self.parent().linear_combination(zip(c_coordinates, basis))
+
+
+ def trace(self):
+ """
+ Return my trace, the sum of my eigenvalues.
+
+ EXAMPLES::
+
+ sage: J = JordanSpinSimpleEJA(3)
+ sage: e0,e1,e2 = J.gens()
+ sage: x = e0 + e1 + e2
+ sage: x.trace()
+ 2
+
+ """
+ cs = self.characteristic_polynomial().coefficients(sparse=False)
+ if len(cs) >= 2:
+ return -1*cs[-2]
+ else:
+ raise ValueError('charpoly had fewer than 2 coefficients')
+
+
+ def trace_inner_product(self, other):
+ """
+ Return the trace inner product of myself and ``other``.
+ """
+ if not other in self.parent():
+ raise TypeError("'other' must live in the same algebra")
+
+ return (self*other).trace()
+
+
+def eja_rn(dimension, field=QQ):
+ """
+ Return the Euclidean Jordan Algebra corresponding to the set
+ `R^n` under the Hadamard product.
+
+ EXAMPLES:
+
+ This multiplication table can be verified by hand::
+
+ sage: J = eja_rn(3)
+ sage: e0,e1,e2 = J.gens()
+ sage: e0*e0
+ e0