2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
9 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
10 from sage
.combinat
.free_module
import CombinatorialFreeModule
11 from sage
.matrix
.constructor
import matrix
12 from sage
.matrix
.matrix_space
import MatrixSpace
13 from sage
.misc
.cachefunc
import cached_method
14 from sage
.misc
.prandom
import choice
15 from sage
.misc
.table
import table
16 from sage
.modules
.free_module
import FreeModule
, VectorSpace
17 from sage
.rings
.integer_ring
import ZZ
18 from sage
.rings
.number_field
.number_field
import NumberField
, QuadraticField
19 from sage
.rings
.polynomial
.polynomial_ring_constructor
import PolynomialRing
20 from sage
.rings
.rational_field
import QQ
21 from sage
.rings
.real_lazy
import CLF
, RLF
22 from sage
.structure
.element
import is_Matrix
24 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
25 from mjo
.eja
.eja_utils
import _mat2vec
27 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
28 # This is an ugly hack needed to prevent the category framework
29 # from implementing a coercion from our base ring (e.g. the
30 # rationals) into the algebra. First of all -- such a coercion is
31 # nonsense to begin with. But more importantly, it tries to do so
32 # in the category of rings, and since our algebras aren't
33 # associative they generally won't be rings.
34 _no_generic_basering_coercion
= True
46 sage: from mjo.eja.eja_algebra import random_eja
50 By definition, Jordan multiplication commutes::
52 sage: set_random_seed()
53 sage: J = random_eja()
54 sage: x = J.random_element()
55 sage: y = J.random_element()
61 self
._natural
_basis
= natural_basis
63 # TODO: HACK for the charpoly.. needs redesign badly.
64 self
._basis
_normalizers
= None
67 category
= MagmaticAlgebras(field
).FiniteDimensional()
68 category
= category
.WithBasis().Unital()
70 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
72 range(len(mult_table
)),
75 self
.print_options(bracket
='')
77 # The multiplication table we're given is necessarily in terms
78 # of vectors, because we don't have an algebra yet for
79 # anything to be an element of. However, it's faster in the
80 # long run to have the multiplication table be in terms of
81 # algebra elements. We do this after calling the superclass
82 # constructor so that from_vector() knows what to do.
83 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
84 for ls
in mult_table
]
87 def _element_constructor_(self
, elt
):
89 Construct an element of this algebra from its natural
92 This gets called only after the parent element _call_ method
93 fails to find a coercion for the argument.
97 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
98 ....: RealCartesianProductEJA,
99 ....: RealSymmetricEJA)
103 The identity in `S^n` is converted to the identity in the EJA::
105 sage: J = RealSymmetricEJA(3)
106 sage: I = matrix.identity(QQ,3)
107 sage: J(I) == J.one()
110 This skew-symmetric matrix can't be represented in the EJA::
112 sage: J = RealSymmetricEJA(3)
113 sage: A = matrix(QQ,3, lambda i,j: i-j)
115 Traceback (most recent call last):
117 ArithmeticError: vector is not in free module
121 Ensure that we can convert any element of the two non-matrix
122 simple algebras (whose natural representations are their usual
123 vector representations) back and forth faithfully::
125 sage: set_random_seed()
126 sage: J = RealCartesianProductEJA(5)
127 sage: x = J.random_element()
128 sage: J(x.to_vector().column()) == x
130 sage: J = JordanSpinEJA(5)
131 sage: x = J.random_element()
132 sage: J(x.to_vector().column()) == x
137 # The superclass implementation of random_element()
138 # needs to be able to coerce "0" into the algebra.
141 natural_basis
= self
.natural_basis()
142 basis_space
= natural_basis
[0].matrix_space()
143 if elt
not in basis_space
:
144 raise ValueError("not a naturally-represented algebra element")
146 # Thanks for nothing! Matrix spaces aren't vector spaces in
147 # Sage, so we have to figure out its natural-basis coordinates
148 # ourselves. We use the basis space's ring instead of the
149 # element's ring because the basis space might be an algebraic
150 # closure whereas the base ring of the 3-by-3 identity matrix
151 # could be QQ instead of QQbar.
152 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
153 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
154 coords
= W
.coordinate_vector(_mat2vec(elt
))
155 return self
.from_vector(coords
)
160 Return a string representation of ``self``.
164 sage: from mjo.eja.eja_algebra import JordanSpinEJA
168 Ensure that it says what we think it says::
170 sage: JordanSpinEJA(2, field=QQ)
171 Euclidean Jordan algebra of dimension 2 over Rational Field
172 sage: JordanSpinEJA(3, field=RDF)
173 Euclidean Jordan algebra of dimension 3 over Real Double Field
176 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
177 return fmt
.format(self
.dimension(), self
.base_ring())
179 def product_on_basis(self
, i
, j
):
180 return self
._multiplication
_table
[i
][j
]
182 def _a_regular_element(self
):
184 Guess a regular element. Needed to compute the basis for our
185 characteristic polynomial coefficients.
189 sage: from mjo.eja.eja_algebra import random_eja
193 Ensure that this hacky method succeeds for every algebra that we
194 know how to construct::
196 sage: set_random_seed()
197 sage: J = random_eja()
198 sage: J._a_regular_element().is_regular()
203 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
204 if not z
.is_regular():
205 raise ValueError("don't know a regular element")
210 def _charpoly_basis_space(self
):
212 Return the vector space spanned by the basis used in our
213 characteristic polynomial coefficients. This is used not only to
214 compute those coefficients, but also any time we need to
215 evaluate the coefficients (like when we compute the trace or
218 z
= self
._a
_regular
_element
()
219 # Don't use the parent vector space directly here in case this
220 # happens to be a subalgebra. In that case, we would be e.g.
221 # two-dimensional but span_of_basis() would expect three
223 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
224 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
225 V1
= V
.span_of_basis( basis
)
226 b
= (V1
.basis() + V1
.complement().basis())
227 return V
.span_of_basis(b
)
232 def _charpoly_coeff(self
, i
):
234 Return the coefficient polynomial "a_{i}" of this algebra's
235 general characteristic polynomial.
237 Having this be a separate cached method lets us compute and
238 store the trace/determinant (a_{r-1} and a_{0} respectively)
239 separate from the entire characteristic polynomial.
241 if self
._basis
_normalizers
is not None:
242 # Must be a matrix class?
243 # WARNING/TODO: this whole mess is mis-designed.
244 n
= self
.natural_basis_space().nrows()
245 field
= self
.base_ring().base_ring() # yeeeeaaaahhh
246 J
= self
.__class
__(n
, field
, False)
247 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
248 p
= J
._charpoly
_coeff
(i
)
249 # p might be missing some vars, have to substitute "optionally"
250 pairs
= zip(x
.base_ring().gens(), self
._basis
_normalizers
)
251 substitutions
= { v: v*c for (v,c) in pairs }
252 return p
.subs(substitutions
)
254 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
255 R
= A_of_x
.base_ring()
257 # Guaranteed by theory
260 # Danger: the in-place modification is done for performance
261 # reasons (reconstructing a matrix with huge polynomial
262 # entries is slow), but I don't know how cached_method works,
263 # so it's highly possible that we're modifying some global
264 # list variable by reference, here. In other words, you
265 # probably shouldn't call this method twice on the same
266 # algebra, at the same time, in two threads
267 Ai_orig
= A_of_x
.column(i
)
268 A_of_x
.set_column(i
,xr
)
269 numerator
= A_of_x
.det()
270 A_of_x
.set_column(i
,Ai_orig
)
272 # We're relying on the theory here to ensure that each a_i is
273 # indeed back in R, and the added negative signs are to make
274 # the whole charpoly expression sum to zero.
275 return R(-numerator
/detA
)
279 def _charpoly_matrix_system(self
):
281 Compute the matrix whose entries A_ij are polynomials in
282 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
283 corresponding to `x^r` and the determinent of the matrix A =
284 [A_ij]. In other words, all of the fixed (cachable) data needed
285 to compute the coefficients of the characteristic polynomial.
290 # Turn my vector space into a module so that "vectors" can
291 # have multivatiate polynomial entries.
292 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
293 R
= PolynomialRing(self
.base_ring(), names
)
295 # Using change_ring() on the parent's vector space doesn't work
296 # here because, in a subalgebra, that vector space has a basis
297 # and change_ring() tries to bring the basis along with it. And
298 # that doesn't work unless the new ring is a PID, which it usually
302 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
306 # And figure out the "left multiplication by x" matrix in
309 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
310 for i
in range(n
) ] # don't recompute these!
312 ek
= self
.monomial(k
).to_vector()
314 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
315 for i
in range(n
) ) )
316 Lx
= matrix
.column(R
, lmbx_cols
)
318 # Now we can compute powers of x "symbolically"
319 x_powers
= [self
.one().to_vector(), x
]
320 for d
in range(2, r
+1):
321 x_powers
.append( Lx
*(x_powers
[-1]) )
323 idmat
= matrix
.identity(R
, n
)
325 W
= self
._charpoly
_basis
_space
()
326 W
= W
.change_ring(R
.fraction_field())
328 # Starting with the standard coordinates x = (X1,X2,...,Xn)
329 # and then converting the entries to W-coordinates allows us
330 # to pass in the standard coordinates to the charpoly and get
331 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
334 # W.coordinates(x^2) eval'd at (standard z-coords)
338 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
340 # We want the middle equivalent thing in our matrix, but use
341 # the first equivalent thing instead so that we can pass in
342 # standard coordinates.
343 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
344 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
345 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
346 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
350 def characteristic_polynomial(self
):
352 Return a characteristic polynomial that works for all elements
355 The resulting polynomial has `n+1` variables, where `n` is the
356 dimension of this algebra. The first `n` variables correspond to
357 the coordinates of an algebra element: when evaluated at the
358 coordinates of an algebra element with respect to a certain
359 basis, the result is a univariate polynomial (in the one
360 remaining variable ``t``), namely the characteristic polynomial
365 sage: from mjo.eja.eja_algebra import JordanSpinEJA
369 The characteristic polynomial in the spin algebra is given in
370 Alizadeh, Example 11.11::
372 sage: J = JordanSpinEJA(3)
373 sage: p = J.characteristic_polynomial(); p
374 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
375 sage: xvec = J.one().to_vector()
383 # The list of coefficient polynomials a_1, a_2, ..., a_n.
384 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
386 # We go to a bit of trouble here to reorder the
387 # indeterminates, so that it's easier to evaluate the
388 # characteristic polynomial at x's coordinates and get back
389 # something in terms of t, which is what we want.
391 S
= PolynomialRing(self
.base_ring(),'t')
393 S
= PolynomialRing(S
, R
.variable_names())
396 # Note: all entries past the rth should be zero. The
397 # coefficient of the highest power (x^r) is 1, but it doesn't
398 # appear in the solution vector which contains coefficients
399 # for the other powers (to make them sum to x^r).
401 a
[r
] = 1 # corresponds to x^r
403 # When the rank is equal to the dimension, trying to
404 # assign a[r] goes out-of-bounds.
405 a
.append(1) # corresponds to x^r
407 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
410 def inner_product(self
, x
, y
):
412 The inner product associated with this Euclidean Jordan algebra.
414 Defaults to the trace inner product, but can be overridden by
415 subclasses if they are sure that the necessary properties are
420 sage: from mjo.eja.eja_algebra import random_eja
424 The inner product must satisfy its axiom for this algebra to truly
425 be a Euclidean Jordan Algebra::
427 sage: set_random_seed()
428 sage: J = random_eja()
429 sage: x = J.random_element()
430 sage: y = J.random_element()
431 sage: z = J.random_element()
432 sage: (x*y).inner_product(z) == y.inner_product(x*z)
436 X
= x
.natural_representation()
437 Y
= y
.natural_representation()
438 return self
.natural_inner_product(X
,Y
)
441 def is_trivial(self
):
443 Return whether or not this algebra is trivial.
445 A trivial algebra contains only the zero element.
449 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
453 sage: J = ComplexHermitianEJA(3)
456 sage: A = J.zero().subalgebra_generated_by()
461 return self
.dimension() == 0
464 def multiplication_table(self
):
466 Return a visual representation of this algebra's multiplication
467 table (on basis elements).
471 sage: from mjo.eja.eja_algebra import JordanSpinEJA
475 sage: J = JordanSpinEJA(4)
476 sage: J.multiplication_table()
477 +----++----+----+----+----+
478 | * || e0 | e1 | e2 | e3 |
479 +====++====+====+====+====+
480 | e0 || e0 | e1 | e2 | e3 |
481 +----++----+----+----+----+
482 | e1 || e1 | e0 | 0 | 0 |
483 +----++----+----+----+----+
484 | e2 || e2 | 0 | e0 | 0 |
485 +----++----+----+----+----+
486 | e3 || e3 | 0 | 0 | e0 |
487 +----++----+----+----+----+
490 M
= list(self
._multiplication
_table
) # copy
491 for i
in range(len(M
)):
492 # M had better be "square"
493 M
[i
] = [self
.monomial(i
)] + M
[i
]
494 M
= [["*"] + list(self
.gens())] + M
495 return table(M
, header_row
=True, header_column
=True, frame
=True)
498 def natural_basis(self
):
500 Return a more-natural representation of this algebra's basis.
502 Every finite-dimensional Euclidean Jordan Algebra is a direct
503 sum of five simple algebras, four of which comprise Hermitian
504 matrices. This method returns the original "natural" basis
505 for our underlying vector space. (Typically, the natural basis
506 is used to construct the multiplication table in the first place.)
508 Note that this will always return a matrix. The standard basis
509 in `R^n` will be returned as `n`-by-`1` column matrices.
513 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
514 ....: RealSymmetricEJA)
518 sage: J = RealSymmetricEJA(2)
520 Finite family {0: e0, 1: e1, 2: e2}
521 sage: J.natural_basis()
523 [1 0] [ 0 1/2*sqrt2] [0 0]
524 [0 0], [1/2*sqrt2 0], [0 1]
529 sage: J = JordanSpinEJA(2)
531 Finite family {0: e0, 1: e1}
532 sage: J.natural_basis()
539 if self
._natural
_basis
is None:
540 M
= self
.natural_basis_space()
541 return tuple( M(b
.to_vector()) for b
in self
.basis() )
543 return self
._natural
_basis
546 def natural_basis_space(self
):
548 Return the matrix space in which this algebra's natural basis
551 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
552 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
554 return self
._natural
_basis
[0].matrix_space()
558 def natural_inner_product(X
,Y
):
560 Compute the inner product of two naturally-represented elements.
562 For example in the real symmetric matrix EJA, this will compute
563 the trace inner-product of two n-by-n symmetric matrices. The
564 default should work for the real cartesian product EJA, the
565 Jordan spin EJA, and the real symmetric matrices. The others
566 will have to be overridden.
568 return (X
.conjugate_transpose()*Y
).trace()
574 Return the unit element of this algebra.
578 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
583 sage: J = RealCartesianProductEJA(5)
585 e0 + e1 + e2 + e3 + e4
589 The identity element acts like the identity::
591 sage: set_random_seed()
592 sage: J = random_eja()
593 sage: x = J.random_element()
594 sage: J.one()*x == x and x*J.one() == x
597 The matrix of the unit element's operator is the identity::
599 sage: set_random_seed()
600 sage: J = random_eja()
601 sage: actual = J.one().operator().matrix()
602 sage: expected = matrix.identity(J.base_ring(), J.dimension())
603 sage: actual == expected
607 # We can brute-force compute the matrices of the operators
608 # that correspond to the basis elements of this algebra.
609 # If some linear combination of those basis elements is the
610 # algebra identity, then the same linear combination of
611 # their matrices has to be the identity matrix.
613 # Of course, matrices aren't vectors in sage, so we have to
614 # appeal to the "long vectors" isometry.
615 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
617 # Now we use basis linear algebra to find the coefficients,
618 # of the matrices-as-vectors-linear-combination, which should
619 # work for the original algebra basis too.
620 A
= matrix
.column(self
.base_ring(), oper_vecs
)
622 # We used the isometry on the left-hand side already, but we
623 # still need to do it for the right-hand side. Recall that we
624 # wanted something that summed to the identity matrix.
625 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
627 # Now if there's an identity element in the algebra, this should work.
628 coeffs
= A
.solve_right(b
)
629 return self
.linear_combination(zip(self
.gens(), coeffs
))
632 def random_element(self
):
633 # Temporary workaround for https://trac.sagemath.org/ticket/28327
634 if self
.is_trivial():
637 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
638 return s
.random_element()
643 Return the rank of this EJA.
647 The author knows of no algorithm to compute the rank of an EJA
648 where only the multiplication table is known. In lieu of one, we
649 require the rank to be specified when the algebra is created,
650 and simply pass along that number here.
654 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
655 ....: RealSymmetricEJA,
656 ....: ComplexHermitianEJA,
657 ....: QuaternionHermitianEJA,
662 The rank of the Jordan spin algebra is always two::
664 sage: JordanSpinEJA(2).rank()
666 sage: JordanSpinEJA(3).rank()
668 sage: JordanSpinEJA(4).rank()
671 The rank of the `n`-by-`n` Hermitian real, complex, or
672 quaternion matrices is `n`::
674 sage: RealSymmetricEJA(2).rank()
676 sage: ComplexHermitianEJA(2).rank()
678 sage: QuaternionHermitianEJA(2).rank()
680 sage: RealSymmetricEJA(5).rank()
682 sage: ComplexHermitianEJA(5).rank()
684 sage: QuaternionHermitianEJA(5).rank()
689 Ensure that every EJA that we know how to construct has a
690 positive integer rank::
692 sage: set_random_seed()
693 sage: r = random_eja().rank()
694 sage: r in ZZ and r > 0
701 def vector_space(self
):
703 Return the vector space that underlies this algebra.
707 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
711 sage: J = RealSymmetricEJA(2)
712 sage: J.vector_space()
713 Vector space of dimension 3 over...
716 return self
.zero().to_vector().parent().ambient_vector_space()
719 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
722 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
724 Return the Euclidean Jordan Algebra corresponding to the set
725 `R^n` under the Hadamard product.
727 Note: this is nothing more than the Cartesian product of ``n``
728 copies of the spin algebra. Once Cartesian product algebras
729 are implemented, this can go.
733 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
737 This multiplication table can be verified by hand::
739 sage: J = RealCartesianProductEJA(3)
740 sage: e0,e1,e2 = J.gens()
756 We can change the generator prefix::
758 sage: RealCartesianProductEJA(3, prefix='r').gens()
761 Our inner product satisfies the Jordan axiom::
763 sage: set_random_seed()
764 sage: n = ZZ.random_element(1,5)
765 sage: J = RealCartesianProductEJA(n)
766 sage: x = J.random_element()
767 sage: y = J.random_element()
768 sage: z = J.random_element()
769 sage: (x*y).inner_product(z) == y.inner_product(x*z)
773 def __init__(self
, n
, field
=QQ
, **kwargs
):
774 V
= VectorSpace(field
, n
)
775 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
778 fdeja
= super(RealCartesianProductEJA
, self
)
779 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
781 def inner_product(self
, x
, y
):
783 Faster to reimplement than to use natural representations.
787 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
791 Ensure that this is the usual inner product for the algebras
794 sage: set_random_seed()
795 sage: n = ZZ.random_element(1,5)
796 sage: J = RealCartesianProductEJA(n)
797 sage: x = J.random_element()
798 sage: y = J.random_element()
799 sage: X = x.natural_representation()
800 sage: Y = y.natural_representation()
801 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
805 return x
.to_vector().inner_product(y
.to_vector())
810 Return a "random" finite-dimensional Euclidean Jordan Algebra.
814 For now, we choose a random natural number ``n`` (greater than zero)
815 and then give you back one of the following:
817 * The cartesian product of the rational numbers ``n`` times; this is
818 ``QQ^n`` with the Hadamard product.
820 * The Jordan spin algebra on ``QQ^n``.
822 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
825 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
826 in the space of ``2n``-by-``2n`` real symmetric matrices.
828 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
829 in the space of ``4n``-by-``4n`` real symmetric matrices.
831 Later this might be extended to return Cartesian products of the
836 sage: from mjo.eja.eja_algebra import random_eja
841 Euclidean Jordan algebra of dimension...
845 # The max_n component lets us choose different upper bounds on the
846 # value "n" that gets passed to the constructor. This is needed
847 # because e.g. R^{10} is reasonable to test, while the Hermitian
848 # 10-by-10 quaternion matrices are not.
849 (constructor
, max_n
) = choice([(RealCartesianProductEJA
, 6),
851 (RealSymmetricEJA
, 5),
852 (ComplexHermitianEJA
, 4),
853 (QuaternionHermitianEJA
, 3)])
854 n
= ZZ
.random_element(1, max_n
)
855 return constructor(n
, field
=QQ
)
859 def _real_symmetric_basis(n
, field
):
861 Return a basis for the space of real symmetric n-by-n matrices.
865 sage: from mjo.eja.eja_algebra import _real_symmetric_basis
869 sage: set_random_seed()
870 sage: n = ZZ.random_element(1,5)
871 sage: B = _real_symmetric_basis(n, QQ)
872 sage: all( M.is_symmetric() for M in B)
876 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
880 for j
in xrange(i
+1):
881 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
885 Sij
= Eij
+ Eij
.transpose()
890 def _complex_hermitian_basis(n
, field
):
892 Returns a basis for the space of complex Hermitian n-by-n matrices.
894 Why do we embed these? Basically, because all of numerical linear
895 algebra assumes that you're working with vectors consisting of `n`
896 entries from a field and scalars from the same field. There's no way
897 to tell SageMath that (for example) the vectors contain complex
898 numbers, while the scalar field is real.
902 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
906 sage: set_random_seed()
907 sage: n = ZZ.random_element(1,5)
908 sage: field = QuadraticField(2, 'sqrt2')
909 sage: B = _complex_hermitian_basis(n, field)
910 sage: all( M.is_symmetric() for M in B)
914 R
= PolynomialRing(field
, 'z')
916 F
= NumberField(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
919 # This is like the symmetric case, but we need to be careful:
921 # * We want conjugate-symmetry, not just symmetry.
922 # * The diagonal will (as a result) be real.
926 for j
in xrange(i
+1):
927 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
929 Sij
= _embed_complex_matrix(Eij
)
932 # The second one has a minus because it's conjugated.
933 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
935 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
938 # Since we embedded these, we can drop back to the "field" that we
939 # started with instead of the complex extension "F".
940 return tuple( s
.change_ring(field
) for s
in S
)
944 def _quaternion_hermitian_basis(n
, field
):
946 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
948 Why do we embed these? Basically, because all of numerical linear
949 algebra assumes that you're working with vectors consisting of `n`
950 entries from a field and scalars from the same field. There's no way
951 to tell SageMath that (for example) the vectors contain complex
952 numbers, while the scalar field is real.
956 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
960 sage: set_random_seed()
961 sage: n = ZZ.random_element(1,5)
962 sage: B = _quaternion_hermitian_basis(n, QQ)
963 sage: all( M.is_symmetric() for M in B )
967 Q
= QuaternionAlgebra(QQ
,-1,-1)
970 # This is like the symmetric case, but we need to be careful:
972 # * We want conjugate-symmetry, not just symmetry.
973 # * The diagonal will (as a result) be real.
977 for j
in xrange(i
+1):
978 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
980 Sij
= _embed_quaternion_matrix(Eij
)
983 # Beware, orthogonal but not normalized! The second,
984 # third, and fourth ones have a minus because they're
986 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
988 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
990 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
992 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
998 def _multiplication_table_from_matrix_basis(basis
):
1000 At least three of the five simple Euclidean Jordan algebras have the
1001 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1002 multiplication on the right is matrix multiplication. Given a basis
1003 for the underlying matrix space, this function returns a
1004 multiplication table (obtained by looping through the basis
1005 elements) for an algebra of those matrices.
1007 # In S^2, for example, we nominally have four coordinates even
1008 # though the space is of dimension three only. The vector space V
1009 # is supposed to hold the entire long vector, and the subspace W
1010 # of V will be spanned by the vectors that arise from symmetric
1011 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1012 field
= basis
[0].base_ring()
1013 dimension
= basis
[0].nrows()
1015 V
= VectorSpace(field
, dimension
**2)
1016 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1018 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1021 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1022 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1027 def _embed_complex_matrix(M
):
1029 Embed the n-by-n complex matrix ``M`` into the space of real
1030 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1031 bi` to the block matrix ``[[a,b],[-b,a]]``.
1035 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
1039 sage: F = QuadraticField(-1, 'i')
1040 sage: x1 = F(4 - 2*i)
1041 sage: x2 = F(1 + 2*i)
1044 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1045 sage: _embed_complex_matrix(M)
1054 Embedding is a homomorphism (isomorphism, in fact)::
1056 sage: set_random_seed()
1057 sage: n = ZZ.random_element(5)
1058 sage: F = QuadraticField(-1, 'i')
1059 sage: X = random_matrix(F, n)
1060 sage: Y = random_matrix(F, n)
1061 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1062 sage: expected = _embed_complex_matrix(X*Y)
1063 sage: actual == expected
1069 raise ValueError("the matrix 'M' must be square")
1070 field
= M
.base_ring()
1073 a
= z
.vector()[0] # real part, I guess
1074 b
= z
.vector()[1] # imag part, I guess
1075 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1077 # We can drop the imaginaries here.
1078 return matrix
.block(field
.base_ring(), n
, blocks
)
1081 def _unembed_complex_matrix(M
):
1083 The inverse of _embed_complex_matrix().
1087 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
1088 ....: _unembed_complex_matrix)
1092 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1093 ....: [-2, 1, -4, 3],
1094 ....: [ 9, 10, 11, 12],
1095 ....: [-10, 9, -12, 11] ])
1096 sage: _unembed_complex_matrix(A)
1098 [ 10*i + 9 12*i + 11]
1102 Unembedding is the inverse of embedding::
1104 sage: set_random_seed()
1105 sage: F = QuadraticField(-1, 'i')
1106 sage: M = random_matrix(F, 3)
1107 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1113 raise ValueError("the matrix 'M' must be square")
1114 if not n
.mod(2).is_zero():
1115 raise ValueError("the matrix 'M' must be a complex embedding")
1117 field
= M
.base_ring() # This should already have sqrt2
1118 R
= PolynomialRing(field
, 'z')
1120 F
= NumberField(z
**2 + 1,'i', embedding
=CLF(-1).sqrt())
1123 # Go top-left to bottom-right (reading order), converting every
1124 # 2-by-2 block we see to a single complex element.
1126 for k
in xrange(n
/2):
1127 for j
in xrange(n
/2):
1128 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1129 if submat
[0,0] != submat
[1,1]:
1130 raise ValueError('bad on-diagonal submatrix')
1131 if submat
[0,1] != -submat
[1,0]:
1132 raise ValueError('bad off-diagonal submatrix')
1133 z
= submat
[0,0] + submat
[0,1]*i
1136 return matrix(F
, n
/2, elements
)
1139 def _embed_quaternion_matrix(M
):
1141 Embed the n-by-n quaternion matrix ``M`` into the space of real
1142 matrices of size 4n-by-4n by first sending each quaternion entry
1143 `z = a + bi + cj + dk` to the block-complex matrix
1144 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1149 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
1153 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1154 sage: i,j,k = Q.gens()
1155 sage: x = 1 + 2*i + 3*j + 4*k
1156 sage: M = matrix(Q, 1, [[x]])
1157 sage: _embed_quaternion_matrix(M)
1163 Embedding is a homomorphism (isomorphism, in fact)::
1165 sage: set_random_seed()
1166 sage: n = ZZ.random_element(5)
1167 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1168 sage: X = random_matrix(Q, n)
1169 sage: Y = random_matrix(Q, n)
1170 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1171 sage: expected = _embed_quaternion_matrix(X*Y)
1172 sage: actual == expected
1176 quaternions
= M
.base_ring()
1179 raise ValueError("the matrix 'M' must be square")
1181 F
= QuadraticField(-1, 'i')
1186 t
= z
.coefficient_tuple()
1191 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1192 [-c
+ d
*i
, a
- b
*i
]])
1193 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1195 # We should have real entries by now, so use the realest field
1196 # we've got for the return value.
1197 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1200 def _unembed_quaternion_matrix(M
):
1202 The inverse of _embed_quaternion_matrix().
1206 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1207 ....: _unembed_quaternion_matrix)
1211 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1212 ....: [-2, 1, -4, 3],
1213 ....: [-3, 4, 1, -2],
1214 ....: [-4, -3, 2, 1]])
1215 sage: _unembed_quaternion_matrix(M)
1216 [1 + 2*i + 3*j + 4*k]
1220 Unembedding is the inverse of embedding::
1222 sage: set_random_seed()
1223 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1224 sage: M = random_matrix(Q, 3)
1225 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1231 raise ValueError("the matrix 'M' must be square")
1232 if not n
.mod(4).is_zero():
1233 raise ValueError("the matrix 'M' must be a complex embedding")
1235 # Use the base ring of the matrix to ensure that its entries can be
1236 # multiplied by elements of the quaternion algebra.
1237 field
= M
.base_ring()
1238 Q
= QuaternionAlgebra(field
,-1,-1)
1241 # Go top-left to bottom-right (reading order), converting every
1242 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1245 for l
in xrange(n
/4):
1246 for m
in xrange(n
/4):
1247 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1248 if submat
[0,0] != submat
[1,1].conjugate():
1249 raise ValueError('bad on-diagonal submatrix')
1250 if submat
[0,1] != -submat
[1,0].conjugate():
1251 raise ValueError('bad off-diagonal submatrix')
1252 z
= submat
[0,0].vector()[0] # real part
1253 z
+= submat
[0,0].vector()[1]*i
# imag part
1254 z
+= submat
[0,1].vector()[0]*j
# real part
1255 z
+= submat
[0,1].vector()[1]*k
# imag part
1258 return matrix(Q
, n
/4, elements
)
1261 # The inner product used for the real symmetric simple EJA.
1262 # We keep it as a separate function because e.g. the complex
1263 # algebra uses the same inner product, except divided by 2.
1264 def _matrix_ip(X
,Y
):
1265 X_mat
= X
.natural_representation()
1266 Y_mat
= Y
.natural_representation()
1267 return (X_mat
*Y_mat
).trace()
1270 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1272 The rank-n simple EJA consisting of real symmetric n-by-n
1273 matrices, the usual symmetric Jordan product, and the trace inner
1274 product. It has dimension `(n^2 + n)/2` over the reals.
1278 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1282 sage: J = RealSymmetricEJA(2)
1283 sage: e0, e1, e2 = J.gens()
1293 The dimension of this algebra is `(n^2 + n) / 2`::
1295 sage: set_random_seed()
1296 sage: n = ZZ.random_element(1,5)
1297 sage: J = RealSymmetricEJA(n)
1298 sage: J.dimension() == (n^2 + n)/2
1301 The Jordan multiplication is what we think it is::
1303 sage: set_random_seed()
1304 sage: n = ZZ.random_element(1,5)
1305 sage: J = RealSymmetricEJA(n)
1306 sage: x = J.random_element()
1307 sage: y = J.random_element()
1308 sage: actual = (x*y).natural_representation()
1309 sage: X = x.natural_representation()
1310 sage: Y = y.natural_representation()
1311 sage: expected = (X*Y + Y*X)/2
1312 sage: actual == expected
1314 sage: J(expected) == x*y
1317 We can change the generator prefix::
1319 sage: RealSymmetricEJA(3, prefix='q').gens()
1320 (q0, q1, q2, q3, q4, q5)
1322 Our inner product satisfies the Jordan axiom::
1324 sage: set_random_seed()
1325 sage: n = ZZ.random_element(1,5)
1326 sage: J = RealSymmetricEJA(n)
1327 sage: x = J.random_element()
1328 sage: y = J.random_element()
1329 sage: z = J.random_element()
1330 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1333 Our natural basis is normalized with respect to the natural inner
1334 product unless we specify otherwise::
1336 sage: set_random_seed()
1337 sage: n = ZZ.random_element(1,5)
1338 sage: J = RealSymmetricEJA(n)
1339 sage: all( b.norm() == 1 for b in J.gens() )
1342 Since our natural basis is normalized with respect to the natural
1343 inner product, and since we know that this algebra is an EJA, any
1344 left-multiplication operator's matrix will be symmetric because
1345 natural->EJA basis representation is an isometry and within the EJA
1346 the operator is self-adjoint by the Jordan axiom::
1348 sage: set_random_seed()
1349 sage: n = ZZ.random_element(1,5)
1350 sage: x = RealSymmetricEJA(n).random_element()
1351 sage: x.operator().matrix().is_symmetric()
1355 def __init__(self
, n
, field
=QQ
, normalize_basis
=True, **kwargs
):
1356 S
= _real_symmetric_basis(n
, field
)
1358 if n
> 1 and normalize_basis
:
1359 # We'll need sqrt(2) to normalize the basis, and this
1360 # winds up in the multiplication table, so the whole
1361 # algebra needs to be over the field extension.
1362 R
= PolynomialRing(field
, 'z')
1365 if p
.is_irreducible():
1366 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1367 S
= [ s
.change_ring(field
) for s
in S
]
1368 self
._basis
_normalizers
= tuple(
1369 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in S
)
1370 S
= tuple( s
*c
for (s
,c
) in zip(S
,self
._basis
_normalizers
) )
1372 Qs
= _multiplication_table_from_matrix_basis(S
)
1374 fdeja
= super(RealSymmetricEJA
, self
)
1375 return fdeja
.__init
__(field
,
1383 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1385 The rank-n simple EJA consisting of complex Hermitian n-by-n
1386 matrices over the real numbers, the usual symmetric Jordan product,
1387 and the real-part-of-trace inner product. It has dimension `n^2` over
1392 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1396 The dimension of this algebra is `n^2`::
1398 sage: set_random_seed()
1399 sage: n = ZZ.random_element(1,5)
1400 sage: J = ComplexHermitianEJA(n)
1401 sage: J.dimension() == n^2
1404 The Jordan multiplication is what we think it is::
1406 sage: set_random_seed()
1407 sage: n = ZZ.random_element(1,5)
1408 sage: J = ComplexHermitianEJA(n)
1409 sage: x = J.random_element()
1410 sage: y = J.random_element()
1411 sage: actual = (x*y).natural_representation()
1412 sage: X = x.natural_representation()
1413 sage: Y = y.natural_representation()
1414 sage: expected = (X*Y + Y*X)/2
1415 sage: actual == expected
1417 sage: J(expected) == x*y
1420 We can change the generator prefix::
1422 sage: ComplexHermitianEJA(2, prefix='z').gens()
1425 Our inner product satisfies the Jordan axiom::
1427 sage: set_random_seed()
1428 sage: n = ZZ.random_element(1,5)
1429 sage: J = ComplexHermitianEJA(n)
1430 sage: x = J.random_element()
1431 sage: y = J.random_element()
1432 sage: z = J.random_element()
1433 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1436 Our natural basis is normalized with respect to the natural inner
1437 product unless we specify otherwise::
1439 sage: set_random_seed()
1440 sage: n = ZZ.random_element(1,4)
1441 sage: J = ComplexHermitianEJA(n)
1442 sage: all( b.norm() == 1 for b in J.gens() )
1445 Since our natural basis is normalized with respect to the natural
1446 inner product, and since we know that this algebra is an EJA, any
1447 left-multiplication operator's matrix will be symmetric because
1448 natural->EJA basis representation is an isometry and within the EJA
1449 the operator is self-adjoint by the Jordan axiom::
1451 sage: set_random_seed()
1452 sage: n = ZZ.random_element(1,5)
1453 sage: x = ComplexHermitianEJA(n).random_element()
1454 sage: x.operator().matrix().is_symmetric()
1458 def __init__(self
, n
, field
=QQ
, normalize_basis
=True, **kwargs
):
1459 S
= _complex_hermitian_basis(n
, field
)
1461 if n
> 1 and normalize_basis
:
1462 # We'll need sqrt(2) to normalize the basis, and this
1463 # winds up in the multiplication table, so the whole
1464 # algebra needs to be over the field extension.
1465 R
= PolynomialRing(field
, 'z')
1468 if p
.is_irreducible():
1469 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1470 S
= [ s
.change_ring(field
) for s
in S
]
1471 self
._basis
_normalizers
= tuple(
1472 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in S
)
1473 S
= tuple( s
*c
for (s
,c
) in zip(S
,self
._basis
_normalizers
) )
1475 Qs
= _multiplication_table_from_matrix_basis(S
)
1477 fdeja
= super(ComplexHermitianEJA
, self
)
1478 return fdeja
.__init
__(field
,
1486 def natural_inner_product(X
,Y
):
1487 Xu
= _unembed_complex_matrix(X
)
1488 Yu
= _unembed_complex_matrix(Y
)
1489 # The trace need not be real; consider Xu = (i*I) and Yu = I.
1490 return ((Xu
*Yu
).trace()).vector()[0] # real part, I guess
1494 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1496 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1497 matrices, the usual symmetric Jordan product, and the
1498 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1503 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1507 The dimension of this algebra is `n^2`::
1509 sage: set_random_seed()
1510 sage: n = ZZ.random_element(1,4)
1511 sage: J = QuaternionHermitianEJA(n)
1512 sage: J.dimension() == 2*(n^2) - n
1515 The Jordan multiplication is what we think it is::
1517 sage: set_random_seed()
1518 sage: n = ZZ.random_element(1,4)
1519 sage: J = QuaternionHermitianEJA(n)
1520 sage: x = J.random_element()
1521 sage: y = J.random_element()
1522 sage: actual = (x*y).natural_representation()
1523 sage: X = x.natural_representation()
1524 sage: Y = y.natural_representation()
1525 sage: expected = (X*Y + Y*X)/2
1526 sage: actual == expected
1528 sage: J(expected) == x*y
1531 We can change the generator prefix::
1533 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1534 (a0, a1, a2, a3, a4, a5)
1536 Our inner product satisfies the Jordan axiom::
1538 sage: set_random_seed()
1539 sage: n = ZZ.random_element(1,4)
1540 sage: J = QuaternionHermitianEJA(n)
1541 sage: x = J.random_element()
1542 sage: y = J.random_element()
1543 sage: z = J.random_element()
1544 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1547 Our natural basis is normalized with respect to the natural inner
1548 product unless we specify otherwise::
1550 sage: set_random_seed()
1551 sage: n = ZZ.random_element(1,4)
1552 sage: J = QuaternionHermitianEJA(n)
1553 sage: all( b.norm() == 1 for b in J.gens() )
1556 Since our natural basis is normalized with respect to the natural
1557 inner product, and since we know that this algebra is an EJA, any
1558 left-multiplication operator's matrix will be symmetric because
1559 natural->EJA basis representation is an isometry and within the EJA
1560 the operator is self-adjoint by the Jordan axiom::
1562 sage: set_random_seed()
1563 sage: n = ZZ.random_element(1,5)
1564 sage: x = QuaternionHermitianEJA(n).random_element()
1565 sage: x.operator().matrix().is_symmetric()
1569 def __init__(self
, n
, field
=QQ
, normalize_basis
=True, **kwargs
):
1570 S
= _quaternion_hermitian_basis(n
, field
)
1572 if n
> 1 and normalize_basis
:
1573 # We'll need sqrt(2) to normalize the basis, and this
1574 # winds up in the multiplication table, so the whole
1575 # algebra needs to be over the field extension.
1576 R
= PolynomialRing(field
, 'z')
1579 if p
.is_irreducible():
1580 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1581 S
= [ s
.change_ring(field
) for s
in S
]
1582 self
._basis
_normalizers
= tuple(
1583 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in S
)
1584 S
= tuple( s
*c
for (s
,c
) in zip(S
,self
._basis
_normalizers
) )
1586 Qs
= _multiplication_table_from_matrix_basis(S
)
1588 fdeja
= super(QuaternionHermitianEJA
, self
)
1589 return fdeja
.__init
__(field
,
1596 def natural_inner_product(X
,Y
):
1597 Xu
= _unembed_quaternion_matrix(X
)
1598 Yu
= _unembed_quaternion_matrix(Y
)
1599 # The trace need not be real; consider Xu = (i*I) and Yu = I.
1600 # The result will be a quaternion algebra element, which doesn't
1601 # have a "vector" method, but does have coefficient_tuple() method
1602 # that returns the coefficients of 1, i, j, and k -- in that order.
1603 return ((Xu
*Yu
).trace()).coefficient_tuple()[0]
1607 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1609 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1610 with the usual inner product and jordan product ``x*y =
1611 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1616 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1620 This multiplication table can be verified by hand::
1622 sage: J = JordanSpinEJA(4)
1623 sage: e0,e1,e2,e3 = J.gens()
1639 We can change the generator prefix::
1641 sage: JordanSpinEJA(2, prefix='B').gens()
1644 Our inner product satisfies the Jordan axiom::
1646 sage: set_random_seed()
1647 sage: n = ZZ.random_element(1,5)
1648 sage: J = JordanSpinEJA(n)
1649 sage: x = J.random_element()
1650 sage: y = J.random_element()
1651 sage: z = J.random_element()
1652 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1656 def __init__(self
, n
, field
=QQ
, **kwargs
):
1657 V
= VectorSpace(field
, n
)
1658 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1668 z0
= x
.inner_product(y
)
1669 zbar
= y0
*xbar
+ x0
*ybar
1670 z
= V([z0
] + zbar
.list())
1671 mult_table
[i
][j
] = z
1673 # The rank of the spin algebra is two, unless we're in a
1674 # one-dimensional ambient space (because the rank is bounded by
1675 # the ambient dimension).
1676 fdeja
= super(JordanSpinEJA
, self
)
1677 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1679 def inner_product(self
, x
, y
):
1681 Faster to reimplement than to use natural representations.
1685 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1689 Ensure that this is the usual inner product for the algebras
1692 sage: set_random_seed()
1693 sage: n = ZZ.random_element(1,5)
1694 sage: J = JordanSpinEJA(n)
1695 sage: x = J.random_element()
1696 sage: y = J.random_element()
1697 sage: X = x.natural_representation()
1698 sage: Y = y.natural_representation()
1699 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1703 return x
.to_vector().inner_product(y
.to_vector())