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.random_instance()
127 sage: x = J.random_element()
128 sage: J(x.to_vector().column()) == x
130 sage: J = JordanSpinEJA.random_instance()
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
)
159 def _max_test_case_size():
161 Return an integer "size" that is an upper bound on the size of
162 this algebra when it is used in a random test
163 case. Unfortunately, the term "size" is quite vague -- when
164 dealing with `R^n` under either the Hadamard or Jordan spin
165 product, the "size" refers to the dimension `n`. When dealing
166 with a matrix algebra (real symmetric or complex/quaternion
167 Hermitian), it refers to the size of the matrix, which is
168 far less than the dimension of the underlying vector space.
170 We default to five in this class, which is safe in `R^n`. The
171 matrix algebra subclasses (or any class where the "size" is
172 interpreted to be far less than the dimension) should override
173 with a smaller number.
180 Return a string representation of ``self``.
184 sage: from mjo.eja.eja_algebra import JordanSpinEJA
188 Ensure that it says what we think it says::
190 sage: JordanSpinEJA(2, field=QQ)
191 Euclidean Jordan algebra of dimension 2 over Rational Field
192 sage: JordanSpinEJA(3, field=RDF)
193 Euclidean Jordan algebra of dimension 3 over Real Double Field
196 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
197 return fmt
.format(self
.dimension(), self
.base_ring())
199 def product_on_basis(self
, i
, j
):
200 return self
._multiplication
_table
[i
][j
]
202 def _a_regular_element(self
):
204 Guess a regular element. Needed to compute the basis for our
205 characteristic polynomial coefficients.
209 sage: from mjo.eja.eja_algebra import random_eja
213 Ensure that this hacky method succeeds for every algebra that we
214 know how to construct::
216 sage: set_random_seed()
217 sage: J = random_eja()
218 sage: J._a_regular_element().is_regular()
223 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
224 if not z
.is_regular():
225 raise ValueError("don't know a regular element")
230 def _charpoly_basis_space(self
):
232 Return the vector space spanned by the basis used in our
233 characteristic polynomial coefficients. This is used not only to
234 compute those coefficients, but also any time we need to
235 evaluate the coefficients (like when we compute the trace or
238 z
= self
._a
_regular
_element
()
239 # Don't use the parent vector space directly here in case this
240 # happens to be a subalgebra. In that case, we would be e.g.
241 # two-dimensional but span_of_basis() would expect three
243 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
244 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
245 V1
= V
.span_of_basis( basis
)
246 b
= (V1
.basis() + V1
.complement().basis())
247 return V
.span_of_basis(b
)
252 def _charpoly_coeff(self
, i
):
254 Return the coefficient polynomial "a_{i}" of this algebra's
255 general characteristic polynomial.
257 Having this be a separate cached method lets us compute and
258 store the trace/determinant (a_{r-1} and a_{0} respectively)
259 separate from the entire characteristic polynomial.
261 if self
._basis
_normalizers
is not None:
262 # Must be a matrix class?
263 # WARNING/TODO: this whole mess is mis-designed.
264 n
= self
.natural_basis_space().nrows()
265 field
= self
.base_ring().base_ring() # yeeeeaaaahhh
266 J
= self
.__class
__(n
, field
, False)
267 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
268 p
= J
._charpoly
_coeff
(i
)
269 # p might be missing some vars, have to substitute "optionally"
270 pairs
= zip(x
.base_ring().gens(), self
._basis
_normalizers
)
271 substitutions
= { v: v*c for (v,c) in pairs }
272 return p
.subs(substitutions
)
274 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
275 R
= A_of_x
.base_ring()
277 # Guaranteed by theory
280 # Danger: the in-place modification is done for performance
281 # reasons (reconstructing a matrix with huge polynomial
282 # entries is slow), but I don't know how cached_method works,
283 # so it's highly possible that we're modifying some global
284 # list variable by reference, here. In other words, you
285 # probably shouldn't call this method twice on the same
286 # algebra, at the same time, in two threads
287 Ai_orig
= A_of_x
.column(i
)
288 A_of_x
.set_column(i
,xr
)
289 numerator
= A_of_x
.det()
290 A_of_x
.set_column(i
,Ai_orig
)
292 # We're relying on the theory here to ensure that each a_i is
293 # indeed back in R, and the added negative signs are to make
294 # the whole charpoly expression sum to zero.
295 return R(-numerator
/detA
)
299 def _charpoly_matrix_system(self
):
301 Compute the matrix whose entries A_ij are polynomials in
302 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
303 corresponding to `x^r` and the determinent of the matrix A =
304 [A_ij]. In other words, all of the fixed (cachable) data needed
305 to compute the coefficients of the characteristic polynomial.
310 # Turn my vector space into a module so that "vectors" can
311 # have multivatiate polynomial entries.
312 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
313 R
= PolynomialRing(self
.base_ring(), names
)
315 # Using change_ring() on the parent's vector space doesn't work
316 # here because, in a subalgebra, that vector space has a basis
317 # and change_ring() tries to bring the basis along with it. And
318 # that doesn't work unless the new ring is a PID, which it usually
322 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
326 # And figure out the "left multiplication by x" matrix in
329 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
330 for i
in range(n
) ] # don't recompute these!
332 ek
= self
.monomial(k
).to_vector()
334 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
335 for i
in range(n
) ) )
336 Lx
= matrix
.column(R
, lmbx_cols
)
338 # Now we can compute powers of x "symbolically"
339 x_powers
= [self
.one().to_vector(), x
]
340 for d
in range(2, r
+1):
341 x_powers
.append( Lx
*(x_powers
[-1]) )
343 idmat
= matrix
.identity(R
, n
)
345 W
= self
._charpoly
_basis
_space
()
346 W
= W
.change_ring(R
.fraction_field())
348 # Starting with the standard coordinates x = (X1,X2,...,Xn)
349 # and then converting the entries to W-coordinates allows us
350 # to pass in the standard coordinates to the charpoly and get
351 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
354 # W.coordinates(x^2) eval'd at (standard z-coords)
358 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
360 # We want the middle equivalent thing in our matrix, but use
361 # the first equivalent thing instead so that we can pass in
362 # standard coordinates.
363 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
364 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
365 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
366 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
370 def characteristic_polynomial(self
):
372 Return a characteristic polynomial that works for all elements
375 The resulting polynomial has `n+1` variables, where `n` is the
376 dimension of this algebra. The first `n` variables correspond to
377 the coordinates of an algebra element: when evaluated at the
378 coordinates of an algebra element with respect to a certain
379 basis, the result is a univariate polynomial (in the one
380 remaining variable ``t``), namely the characteristic polynomial
385 sage: from mjo.eja.eja_algebra import JordanSpinEJA
389 The characteristic polynomial in the spin algebra is given in
390 Alizadeh, Example 11.11::
392 sage: J = JordanSpinEJA(3)
393 sage: p = J.characteristic_polynomial(); p
394 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
395 sage: xvec = J.one().to_vector()
403 # The list of coefficient polynomials a_1, a_2, ..., a_n.
404 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
406 # We go to a bit of trouble here to reorder the
407 # indeterminates, so that it's easier to evaluate the
408 # characteristic polynomial at x's coordinates and get back
409 # something in terms of t, which is what we want.
411 S
= PolynomialRing(self
.base_ring(),'t')
413 S
= PolynomialRing(S
, R
.variable_names())
416 # Note: all entries past the rth should be zero. The
417 # coefficient of the highest power (x^r) is 1, but it doesn't
418 # appear in the solution vector which contains coefficients
419 # for the other powers (to make them sum to x^r).
421 a
[r
] = 1 # corresponds to x^r
423 # When the rank is equal to the dimension, trying to
424 # assign a[r] goes out-of-bounds.
425 a
.append(1) # corresponds to x^r
427 return sum( a
[k
]*(t
**k
) for k
in range(len(a
)) )
430 def inner_product(self
, x
, y
):
432 The inner product associated with this Euclidean Jordan algebra.
434 Defaults to the trace inner product, but can be overridden by
435 subclasses if they are sure that the necessary properties are
440 sage: from mjo.eja.eja_algebra import random_eja
444 The inner product must satisfy its axiom for this algebra to truly
445 be a Euclidean Jordan Algebra::
447 sage: set_random_seed()
448 sage: J = random_eja()
449 sage: x = J.random_element()
450 sage: y = J.random_element()
451 sage: z = J.random_element()
452 sage: (x*y).inner_product(z) == y.inner_product(x*z)
456 X
= x
.natural_representation()
457 Y
= y
.natural_representation()
458 return self
.natural_inner_product(X
,Y
)
461 def is_trivial(self
):
463 Return whether or not this algebra is trivial.
465 A trivial algebra contains only the zero element.
469 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
473 sage: J = ComplexHermitianEJA(3)
476 sage: A = J.zero().subalgebra_generated_by()
481 return self
.dimension() == 0
484 def multiplication_table(self
):
486 Return a visual representation of this algebra's multiplication
487 table (on basis elements).
491 sage: from mjo.eja.eja_algebra import JordanSpinEJA
495 sage: J = JordanSpinEJA(4)
496 sage: J.multiplication_table()
497 +----++----+----+----+----+
498 | * || e0 | e1 | e2 | e3 |
499 +====++====+====+====+====+
500 | e0 || e0 | e1 | e2 | e3 |
501 +----++----+----+----+----+
502 | e1 || e1 | e0 | 0 | 0 |
503 +----++----+----+----+----+
504 | e2 || e2 | 0 | e0 | 0 |
505 +----++----+----+----+----+
506 | e3 || e3 | 0 | 0 | e0 |
507 +----++----+----+----+----+
510 M
= list(self
._multiplication
_table
) # copy
511 for i
in range(len(M
)):
512 # M had better be "square"
513 M
[i
] = [self
.monomial(i
)] + M
[i
]
514 M
= [["*"] + list(self
.gens())] + M
515 return table(M
, header_row
=True, header_column
=True, frame
=True)
518 def natural_basis(self
):
520 Return a more-natural representation of this algebra's basis.
522 Every finite-dimensional Euclidean Jordan Algebra is a direct
523 sum of five simple algebras, four of which comprise Hermitian
524 matrices. This method returns the original "natural" basis
525 for our underlying vector space. (Typically, the natural basis
526 is used to construct the multiplication table in the first place.)
528 Note that this will always return a matrix. The standard basis
529 in `R^n` will be returned as `n`-by-`1` column matrices.
533 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
534 ....: RealSymmetricEJA)
538 sage: J = RealSymmetricEJA(2)
540 Finite family {0: e0, 1: e1, 2: e2}
541 sage: J.natural_basis()
543 [1 0] [ 0 1/2*sqrt2] [0 0]
544 [0 0], [1/2*sqrt2 0], [0 1]
549 sage: J = JordanSpinEJA(2)
551 Finite family {0: e0, 1: e1}
552 sage: J.natural_basis()
559 if self
._natural
_basis
is None:
560 M
= self
.natural_basis_space()
561 return tuple( M(b
.to_vector()) for b
in self
.basis() )
563 return self
._natural
_basis
566 def natural_basis_space(self
):
568 Return the matrix space in which this algebra's natural basis
571 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
572 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
574 return self
._natural
_basis
[0].matrix_space()
578 def natural_inner_product(X
,Y
):
580 Compute the inner product of two naturally-represented elements.
582 For example in the real symmetric matrix EJA, this will compute
583 the trace inner-product of two n-by-n symmetric matrices. The
584 default should work for the real cartesian product EJA, the
585 Jordan spin EJA, and the real symmetric matrices. The others
586 will have to be overridden.
588 return (X
.conjugate_transpose()*Y
).trace()
594 Return the unit element of this algebra.
598 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
603 sage: J = RealCartesianProductEJA(5)
605 e0 + e1 + e2 + e3 + e4
609 The identity element acts like the identity::
611 sage: set_random_seed()
612 sage: J = random_eja()
613 sage: x = J.random_element()
614 sage: J.one()*x == x and x*J.one() == x
617 The matrix of the unit element's operator is the identity::
619 sage: set_random_seed()
620 sage: J = random_eja()
621 sage: actual = J.one().operator().matrix()
622 sage: expected = matrix.identity(J.base_ring(), J.dimension())
623 sage: actual == expected
627 # We can brute-force compute the matrices of the operators
628 # that correspond to the basis elements of this algebra.
629 # If some linear combination of those basis elements is the
630 # algebra identity, then the same linear combination of
631 # their matrices has to be the identity matrix.
633 # Of course, matrices aren't vectors in sage, so we have to
634 # appeal to the "long vectors" isometry.
635 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
637 # Now we use basis linear algebra to find the coefficients,
638 # of the matrices-as-vectors-linear-combination, which should
639 # work for the original algebra basis too.
640 A
= matrix
.column(self
.base_ring(), oper_vecs
)
642 # We used the isometry on the left-hand side already, but we
643 # still need to do it for the right-hand side. Recall that we
644 # wanted something that summed to the identity matrix.
645 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
647 # Now if there's an identity element in the algebra, this should work.
648 coeffs
= A
.solve_right(b
)
649 return self
.linear_combination(zip(self
.gens(), coeffs
))
652 def random_element(self
):
653 # Temporary workaround for https://trac.sagemath.org/ticket/28327
654 if self
.is_trivial():
657 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
658 return s
.random_element()
662 def random_instance(cls
, field
=QQ
, **kwargs
):
664 Return a random instance of this type of algebra.
666 In subclasses for algebras that we know how to construct, this
667 is a shortcut for constructing test cases and examples.
669 if cls
is FiniteDimensionalEuclideanJordanAlgebra
:
670 # Red flag! But in theory we could do this I guess. The
671 # only finite-dimensional exceptional EJA is the
672 # octononions. So, we could just create an EJA from an
673 # associative matrix algebra (generated by a subset of
674 # elements) with the symmetric product. Or, we could punt
675 # to random_eja() here, override it in our subclasses, and
676 # not worry about it.
677 raise NotImplementedError
679 n
= ZZ
.random_element(1, cls
._max
_test
_case
_size
())
680 return cls(n
, field
, **kwargs
)
685 Return the rank of this EJA.
689 The author knows of no algorithm to compute the rank of an EJA
690 where only the multiplication table is known. In lieu of one, we
691 require the rank to be specified when the algebra is created,
692 and simply pass along that number here.
696 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
697 ....: RealSymmetricEJA,
698 ....: ComplexHermitianEJA,
699 ....: QuaternionHermitianEJA,
704 The rank of the Jordan spin algebra is always two::
706 sage: JordanSpinEJA(2).rank()
708 sage: JordanSpinEJA(3).rank()
710 sage: JordanSpinEJA(4).rank()
713 The rank of the `n`-by-`n` Hermitian real, complex, or
714 quaternion matrices is `n`::
716 sage: RealSymmetricEJA(2).rank()
718 sage: ComplexHermitianEJA(2).rank()
720 sage: QuaternionHermitianEJA(2).rank()
722 sage: RealSymmetricEJA(5).rank()
724 sage: ComplexHermitianEJA(5).rank()
726 sage: QuaternionHermitianEJA(5).rank()
731 Ensure that every EJA that we know how to construct has a
732 positive integer rank::
734 sage: set_random_seed()
735 sage: r = random_eja().rank()
736 sage: r in ZZ and r > 0
743 def vector_space(self
):
745 Return the vector space that underlies this algebra.
749 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
753 sage: J = RealSymmetricEJA(2)
754 sage: J.vector_space()
755 Vector space of dimension 3 over...
758 return self
.zero().to_vector().parent().ambient_vector_space()
761 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
764 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
766 Return the Euclidean Jordan Algebra corresponding to the set
767 `R^n` under the Hadamard product.
769 Note: this is nothing more than the Cartesian product of ``n``
770 copies of the spin algebra. Once Cartesian product algebras
771 are implemented, this can go.
775 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
779 This multiplication table can be verified by hand::
781 sage: J = RealCartesianProductEJA(3)
782 sage: e0,e1,e2 = J.gens()
798 We can change the generator prefix::
800 sage: RealCartesianProductEJA(3, prefix='r').gens()
803 Our inner product satisfies the Jordan axiom::
805 sage: set_random_seed()
806 sage: J = RealCartesianProductEJA.random_instance()
807 sage: x = J.random_element()
808 sage: y = J.random_element()
809 sage: z = J.random_element()
810 sage: (x*y).inner_product(z) == y.inner_product(x*z)
814 def __init__(self
, n
, field
=QQ
, **kwargs
):
815 V
= VectorSpace(field
, n
)
816 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in range(n
) ]
819 fdeja
= super(RealCartesianProductEJA
, self
)
820 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
822 def inner_product(self
, x
, y
):
824 Faster to reimplement than to use natural representations.
828 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
832 Ensure that this is the usual inner product for the algebras
835 sage: set_random_seed()
836 sage: J = RealCartesianProductEJA.random_instance()
837 sage: x = J.random_element()
838 sage: y = J.random_element()
839 sage: X = x.natural_representation()
840 sage: Y = y.natural_representation()
841 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
845 return x
.to_vector().inner_product(y
.to_vector())
850 Return a "random" finite-dimensional Euclidean Jordan Algebra.
854 For now, we choose a random natural number ``n`` (greater than zero)
855 and then give you back one of the following:
857 * The cartesian product of the rational numbers ``n`` times; this is
858 ``QQ^n`` with the Hadamard product.
860 * The Jordan spin algebra on ``QQ^n``.
862 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
865 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
866 in the space of ``2n``-by-``2n`` real symmetric matrices.
868 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
869 in the space of ``4n``-by-``4n`` real symmetric matrices.
871 Later this might be extended to return Cartesian products of the
876 sage: from mjo.eja.eja_algebra import random_eja
881 Euclidean Jordan algebra of dimension...
884 classname
= choice([RealCartesianProductEJA
,
888 QuaternionHermitianEJA
])
889 return classname
.random_instance()
893 def _real_symmetric_basis(n
, field
):
895 Return a basis for the space of real symmetric n-by-n matrices.
899 sage: from mjo.eja.eja_algebra import _real_symmetric_basis
903 sage: set_random_seed()
904 sage: n = ZZ.random_element(1,5)
905 sage: B = _real_symmetric_basis(n, QQ)
906 sage: all( M.is_symmetric() for M in B)
910 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
914 for j
in xrange(i
+1):
915 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
919 Sij
= Eij
+ Eij
.transpose()
924 def _complex_hermitian_basis(n
, field
):
926 Returns a basis for the space of complex Hermitian n-by-n matrices.
928 Why do we embed these? Basically, because all of numerical linear
929 algebra assumes that you're working with vectors consisting of `n`
930 entries from a field and scalars from the same field. There's no way
931 to tell SageMath that (for example) the vectors contain complex
932 numbers, while the scalar field is real.
936 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
940 sage: set_random_seed()
941 sage: n = ZZ.random_element(1,5)
942 sage: field = QuadraticField(2, 'sqrt2')
943 sage: B = _complex_hermitian_basis(n, field)
944 sage: all( M.is_symmetric() for M in B)
948 R
= PolynomialRing(field
, 'z')
950 F
= NumberField(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
953 # This is like the symmetric case, but we need to be careful:
955 # * We want conjugate-symmetry, not just symmetry.
956 # * The diagonal will (as a result) be real.
960 for j
in xrange(i
+1):
961 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
963 Sij
= _embed_complex_matrix(Eij
)
966 # The second one has a minus because it's conjugated.
967 Sij_real
= _embed_complex_matrix(Eij
+ Eij
.transpose())
969 Sij_imag
= _embed_complex_matrix(I
*Eij
- I
*Eij
.transpose())
972 # Since we embedded these, we can drop back to the "field" that we
973 # started with instead of the complex extension "F".
974 return tuple( s
.change_ring(field
) for s
in S
)
978 def _quaternion_hermitian_basis(n
, field
):
980 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
982 Why do we embed these? Basically, because all of numerical linear
983 algebra assumes that you're working with vectors consisting of `n`
984 entries from a field and scalars from the same field. There's no way
985 to tell SageMath that (for example) the vectors contain complex
986 numbers, while the scalar field is real.
990 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
994 sage: set_random_seed()
995 sage: n = ZZ.random_element(1,5)
996 sage: B = _quaternion_hermitian_basis(n, QQ)
997 sage: all( M.is_symmetric() for M in B )
1001 Q
= QuaternionAlgebra(QQ
,-1,-1)
1004 # This is like the symmetric case, but we need to be careful:
1006 # * We want conjugate-symmetry, not just symmetry.
1007 # * The diagonal will (as a result) be real.
1011 for j
in xrange(i
+1):
1012 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1014 Sij
= _embed_quaternion_matrix(Eij
)
1017 # Beware, orthogonal but not normalized! The second,
1018 # third, and fourth ones have a minus because they're
1020 Sij_real
= _embed_quaternion_matrix(Eij
+ Eij
.transpose())
1022 Sij_I
= _embed_quaternion_matrix(I
*Eij
- I
*Eij
.transpose())
1024 Sij_J
= _embed_quaternion_matrix(J
*Eij
- J
*Eij
.transpose())
1026 Sij_K
= _embed_quaternion_matrix(K
*Eij
- K
*Eij
.transpose())
1032 def _multiplication_table_from_matrix_basis(basis
):
1034 At least three of the five simple Euclidean Jordan algebras have the
1035 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1036 multiplication on the right is matrix multiplication. Given a basis
1037 for the underlying matrix space, this function returns a
1038 multiplication table (obtained by looping through the basis
1039 elements) for an algebra of those matrices.
1041 # In S^2, for example, we nominally have four coordinates even
1042 # though the space is of dimension three only. The vector space V
1043 # is supposed to hold the entire long vector, and the subspace W
1044 # of V will be spanned by the vectors that arise from symmetric
1045 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1046 field
= basis
[0].base_ring()
1047 dimension
= basis
[0].nrows()
1049 V
= VectorSpace(field
, dimension
**2)
1050 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
1052 mult_table
= [[W
.zero() for j
in range(n
)] for i
in range(n
)]
1055 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
1056 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
1061 def _embed_complex_matrix(M
):
1063 Embed the n-by-n complex matrix ``M`` into the space of real
1064 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1065 bi` to the block matrix ``[[a,b],[-b,a]]``.
1069 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
1070 ....: ComplexHermitianEJA)
1074 sage: F = QuadraticField(-1, 'i')
1075 sage: x1 = F(4 - 2*i)
1076 sage: x2 = F(1 + 2*i)
1079 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1080 sage: _embed_complex_matrix(M)
1089 Embedding is a homomorphism (isomorphism, in fact)::
1091 sage: set_random_seed()
1092 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1093 sage: n = ZZ.random_element(n_max)
1094 sage: F = QuadraticField(-1, 'i')
1095 sage: X = random_matrix(F, n)
1096 sage: Y = random_matrix(F, n)
1097 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1098 sage: expected = _embed_complex_matrix(X*Y)
1099 sage: actual == expected
1105 raise ValueError("the matrix 'M' must be square")
1106 field
= M
.base_ring()
1109 a
= z
.vector()[0] # real part, I guess
1110 b
= z
.vector()[1] # imag part, I guess
1111 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1113 # We can drop the imaginaries here.
1114 return matrix
.block(field
.base_ring(), n
, blocks
)
1117 def _unembed_complex_matrix(M
):
1119 The inverse of _embed_complex_matrix().
1123 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
1124 ....: _unembed_complex_matrix)
1128 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1129 ....: [-2, 1, -4, 3],
1130 ....: [ 9, 10, 11, 12],
1131 ....: [-10, 9, -12, 11] ])
1132 sage: _unembed_complex_matrix(A)
1134 [ 10*i + 9 12*i + 11]
1138 Unembedding is the inverse of embedding::
1140 sage: set_random_seed()
1141 sage: F = QuadraticField(-1, 'i')
1142 sage: M = random_matrix(F, 3)
1143 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1149 raise ValueError("the matrix 'M' must be square")
1150 if not n
.mod(2).is_zero():
1151 raise ValueError("the matrix 'M' must be a complex embedding")
1153 field
= M
.base_ring() # This should already have sqrt2
1154 R
= PolynomialRing(field
, 'z')
1156 F
= NumberField(z
**2 + 1,'i', embedding
=CLF(-1).sqrt())
1159 # Go top-left to bottom-right (reading order), converting every
1160 # 2-by-2 block we see to a single complex element.
1162 for k
in xrange(n
/2):
1163 for j
in xrange(n
/2):
1164 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1165 if submat
[0,0] != submat
[1,1]:
1166 raise ValueError('bad on-diagonal submatrix')
1167 if submat
[0,1] != -submat
[1,0]:
1168 raise ValueError('bad off-diagonal submatrix')
1169 z
= submat
[0,0] + submat
[0,1]*i
1172 return matrix(F
, n
/2, elements
)
1175 def _embed_quaternion_matrix(M
):
1177 Embed the n-by-n quaternion matrix ``M`` into the space of real
1178 matrices of size 4n-by-4n by first sending each quaternion entry
1179 `z = a + bi + cj + dk` to the block-complex matrix
1180 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1185 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1186 ....: QuaternionHermitianEJA)
1190 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1191 sage: i,j,k = Q.gens()
1192 sage: x = 1 + 2*i + 3*j + 4*k
1193 sage: M = matrix(Q, 1, [[x]])
1194 sage: _embed_quaternion_matrix(M)
1200 Embedding is a homomorphism (isomorphism, in fact)::
1202 sage: set_random_seed()
1203 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1204 sage: n = ZZ.random_element(n_max)
1205 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1206 sage: X = random_matrix(Q, n)
1207 sage: Y = random_matrix(Q, n)
1208 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1209 sage: expected = _embed_quaternion_matrix(X*Y)
1210 sage: actual == expected
1214 quaternions
= M
.base_ring()
1217 raise ValueError("the matrix 'M' must be square")
1219 F
= QuadraticField(-1, 'i')
1224 t
= z
.coefficient_tuple()
1229 cplx_matrix
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1230 [-c
+ d
*i
, a
- b
*i
]])
1231 blocks
.append(_embed_complex_matrix(cplx_matrix
))
1233 # We should have real entries by now, so use the realest field
1234 # we've got for the return value.
1235 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1238 def _unembed_quaternion_matrix(M
):
1240 The inverse of _embed_quaternion_matrix().
1244 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1245 ....: _unembed_quaternion_matrix)
1249 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1250 ....: [-2, 1, -4, 3],
1251 ....: [-3, 4, 1, -2],
1252 ....: [-4, -3, 2, 1]])
1253 sage: _unembed_quaternion_matrix(M)
1254 [1 + 2*i + 3*j + 4*k]
1258 Unembedding is the inverse of embedding::
1260 sage: set_random_seed()
1261 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1262 sage: M = random_matrix(Q, 3)
1263 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1269 raise ValueError("the matrix 'M' must be square")
1270 if not n
.mod(4).is_zero():
1271 raise ValueError("the matrix 'M' must be a complex embedding")
1273 # Use the base ring of the matrix to ensure that its entries can be
1274 # multiplied by elements of the quaternion algebra.
1275 field
= M
.base_ring()
1276 Q
= QuaternionAlgebra(field
,-1,-1)
1279 # Go top-left to bottom-right (reading order), converting every
1280 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1283 for l
in xrange(n
/4):
1284 for m
in xrange(n
/4):
1285 submat
= _unembed_complex_matrix(M
[4*l
:4*l
+4,4*m
:4*m
+4])
1286 if submat
[0,0] != submat
[1,1].conjugate():
1287 raise ValueError('bad on-diagonal submatrix')
1288 if submat
[0,1] != -submat
[1,0].conjugate():
1289 raise ValueError('bad off-diagonal submatrix')
1290 z
= submat
[0,0].vector()[0] # real part
1291 z
+= submat
[0,0].vector()[1]*i
# imag part
1292 z
+= submat
[0,1].vector()[0]*j
# real part
1293 z
+= submat
[0,1].vector()[1]*k
# imag part
1296 return matrix(Q
, n
/4, elements
)
1299 # The inner product used for the real symmetric simple EJA.
1300 # We keep it as a separate function because e.g. the complex
1301 # algebra uses the same inner product, except divided by 2.
1302 def _matrix_ip(X
,Y
):
1303 X_mat
= X
.natural_representation()
1304 Y_mat
= Y
.natural_representation()
1305 return (X_mat
*Y_mat
).trace()
1308 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1310 The rank-n simple EJA consisting of real symmetric n-by-n
1311 matrices, the usual symmetric Jordan product, and the trace inner
1312 product. It has dimension `(n^2 + n)/2` over the reals.
1316 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1320 sage: J = RealSymmetricEJA(2)
1321 sage: e0, e1, e2 = J.gens()
1331 The dimension of this algebra is `(n^2 + n) / 2`::
1333 sage: set_random_seed()
1334 sage: n_max = RealSymmetricEJA._max_test_case_size()
1335 sage: n = ZZ.random_element(1, n_max)
1336 sage: J = RealSymmetricEJA(n)
1337 sage: J.dimension() == (n^2 + n)/2
1340 The Jordan multiplication is what we think it is::
1342 sage: set_random_seed()
1343 sage: J = RealSymmetricEJA.random_instance()
1344 sage: x = J.random_element()
1345 sage: y = J.random_element()
1346 sage: actual = (x*y).natural_representation()
1347 sage: X = x.natural_representation()
1348 sage: Y = y.natural_representation()
1349 sage: expected = (X*Y + Y*X)/2
1350 sage: actual == expected
1352 sage: J(expected) == x*y
1355 We can change the generator prefix::
1357 sage: RealSymmetricEJA(3, prefix='q').gens()
1358 (q0, q1, q2, q3, q4, q5)
1360 Our inner product satisfies the Jordan axiom::
1362 sage: set_random_seed()
1363 sage: J = RealSymmetricEJA.random_instance()
1364 sage: x = J.random_element()
1365 sage: y = J.random_element()
1366 sage: z = J.random_element()
1367 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1370 Our natural basis is normalized with respect to the natural inner
1371 product unless we specify otherwise::
1373 sage: set_random_seed()
1374 sage: J = RealSymmetricEJA.random_instance()
1375 sage: all( b.norm() == 1 for b in J.gens() )
1378 Since our natural basis is normalized with respect to the natural
1379 inner product, and since we know that this algebra is an EJA, any
1380 left-multiplication operator's matrix will be symmetric because
1381 natural->EJA basis representation is an isometry and within the EJA
1382 the operator is self-adjoint by the Jordan axiom::
1384 sage: set_random_seed()
1385 sage: x = RealSymmetricEJA.random_instance().random_element()
1386 sage: x.operator().matrix().is_symmetric()
1390 def __init__(self
, n
, field
=QQ
, normalize_basis
=True, **kwargs
):
1391 S
= _real_symmetric_basis(n
, field
)
1393 if n
> 1 and normalize_basis
:
1394 # We'll need sqrt(2) to normalize the basis, and this
1395 # winds up in the multiplication table, so the whole
1396 # algebra needs to be over the field extension.
1397 R
= PolynomialRing(field
, 'z')
1400 if p
.is_irreducible():
1401 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1402 S
= [ s
.change_ring(field
) for s
in S
]
1403 self
._basis
_normalizers
= tuple(
1404 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in S
)
1405 S
= tuple( s
*c
for (s
,c
) in zip(S
,self
._basis
_normalizers
) )
1407 Qs
= _multiplication_table_from_matrix_basis(S
)
1409 fdeja
= super(RealSymmetricEJA
, self
)
1410 return fdeja
.__init
__(field
,
1417 def _max_test_case_size():
1421 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1423 The rank-n simple EJA consisting of complex Hermitian n-by-n
1424 matrices over the real numbers, the usual symmetric Jordan product,
1425 and the real-part-of-trace inner product. It has dimension `n^2` over
1430 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1434 The dimension of this algebra is `n^2`::
1436 sage: set_random_seed()
1437 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1438 sage: n = ZZ.random_element(1, n_max)
1439 sage: J = ComplexHermitianEJA(n)
1440 sage: J.dimension() == n^2
1443 The Jordan multiplication is what we think it is::
1445 sage: set_random_seed()
1446 sage: J = ComplexHermitianEJA.random_instance()
1447 sage: x = J.random_element()
1448 sage: y = J.random_element()
1449 sage: actual = (x*y).natural_representation()
1450 sage: X = x.natural_representation()
1451 sage: Y = y.natural_representation()
1452 sage: expected = (X*Y + Y*X)/2
1453 sage: actual == expected
1455 sage: J(expected) == x*y
1458 We can change the generator prefix::
1460 sage: ComplexHermitianEJA(2, prefix='z').gens()
1463 Our inner product satisfies the Jordan axiom::
1465 sage: set_random_seed()
1466 sage: J = ComplexHermitianEJA.random_instance()
1467 sage: x = J.random_element()
1468 sage: y = J.random_element()
1469 sage: z = J.random_element()
1470 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1473 Our natural basis is normalized with respect to the natural inner
1474 product unless we specify otherwise::
1476 sage: set_random_seed()
1477 sage: J = ComplexHermitianEJA.random_instance()
1478 sage: all( b.norm() == 1 for b in J.gens() )
1481 Since our natural basis is normalized with respect to the natural
1482 inner product, and since we know that this algebra is an EJA, any
1483 left-multiplication operator's matrix will be symmetric because
1484 natural->EJA basis representation is an isometry and within the EJA
1485 the operator is self-adjoint by the Jordan axiom::
1487 sage: set_random_seed()
1488 sage: x = ComplexHermitianEJA.random_instance().random_element()
1489 sage: x.operator().matrix().is_symmetric()
1493 def __init__(self
, n
, field
=QQ
, normalize_basis
=True, **kwargs
):
1494 S
= _complex_hermitian_basis(n
, field
)
1496 if n
> 1 and normalize_basis
:
1497 # We'll need sqrt(2) to normalize the basis, and this
1498 # winds up in the multiplication table, so the whole
1499 # algebra needs to be over the field extension.
1500 R
= PolynomialRing(field
, 'z')
1503 if p
.is_irreducible():
1504 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1505 S
= [ s
.change_ring(field
) for s
in S
]
1506 self
._basis
_normalizers
= tuple(
1507 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in S
)
1508 S
= tuple( s
*c
for (s
,c
) in zip(S
,self
._basis
_normalizers
) )
1510 Qs
= _multiplication_table_from_matrix_basis(S
)
1512 fdeja
= super(ComplexHermitianEJA
, self
)
1513 return fdeja
.__init
__(field
,
1521 def _max_test_case_size():
1525 def natural_inner_product(X
,Y
):
1526 Xu
= _unembed_complex_matrix(X
)
1527 Yu
= _unembed_complex_matrix(Y
)
1528 # The trace need not be real; consider Xu = (i*I) and Yu = I.
1529 return ((Xu
*Yu
).trace()).vector()[0] # real part, I guess
1533 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1535 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1536 matrices, the usual symmetric Jordan product, and the
1537 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1542 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1546 The dimension of this algebra is `2*n^2 - n`::
1548 sage: set_random_seed()
1549 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1550 sage: n = ZZ.random_element(1, n_max)
1551 sage: J = QuaternionHermitianEJA(n)
1552 sage: J.dimension() == 2*(n^2) - n
1555 The Jordan multiplication is what we think it is::
1557 sage: set_random_seed()
1558 sage: J = QuaternionHermitianEJA.random_instance()
1559 sage: x = J.random_element()
1560 sage: y = J.random_element()
1561 sage: actual = (x*y).natural_representation()
1562 sage: X = x.natural_representation()
1563 sage: Y = y.natural_representation()
1564 sage: expected = (X*Y + Y*X)/2
1565 sage: actual == expected
1567 sage: J(expected) == x*y
1570 We can change the generator prefix::
1572 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1573 (a0, a1, a2, a3, a4, a5)
1575 Our inner product satisfies the Jordan axiom::
1577 sage: set_random_seed()
1578 sage: J = QuaternionHermitianEJA.random_instance()
1579 sage: x = J.random_element()
1580 sage: y = J.random_element()
1581 sage: z = J.random_element()
1582 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1585 Our natural basis is normalized with respect to the natural inner
1586 product unless we specify otherwise::
1588 sage: set_random_seed()
1589 sage: J = QuaternionHermitianEJA.random_instance()
1590 sage: all( b.norm() == 1 for b in J.gens() )
1593 Since our natural basis is normalized with respect to the natural
1594 inner product, and since we know that this algebra is an EJA, any
1595 left-multiplication operator's matrix will be symmetric because
1596 natural->EJA basis representation is an isometry and within the EJA
1597 the operator is self-adjoint by the Jordan axiom::
1599 sage: set_random_seed()
1600 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1601 sage: x.operator().matrix().is_symmetric()
1605 def __init__(self
, n
, field
=QQ
, normalize_basis
=True, **kwargs
):
1606 S
= _quaternion_hermitian_basis(n
, field
)
1608 if n
> 1 and normalize_basis
:
1609 # We'll need sqrt(2) to normalize the basis, and this
1610 # winds up in the multiplication table, so the whole
1611 # algebra needs to be over the field extension.
1612 R
= PolynomialRing(field
, 'z')
1615 if p
.is_irreducible():
1616 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
1617 S
= [ s
.change_ring(field
) for s
in S
]
1618 self
._basis
_normalizers
= tuple(
1619 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in S
)
1620 S
= tuple( s
*c
for (s
,c
) in zip(S
,self
._basis
_normalizers
) )
1622 Qs
= _multiplication_table_from_matrix_basis(S
)
1624 fdeja
= super(QuaternionHermitianEJA
, self
)
1625 return fdeja
.__init
__(field
,
1632 def _max_test_case_size():
1636 def natural_inner_product(X
,Y
):
1637 Xu
= _unembed_quaternion_matrix(X
)
1638 Yu
= _unembed_quaternion_matrix(Y
)
1639 # The trace need not be real; consider Xu = (i*I) and Yu = I.
1640 # The result will be a quaternion algebra element, which doesn't
1641 # have a "vector" method, but does have coefficient_tuple() method
1642 # that returns the coefficients of 1, i, j, and k -- in that order.
1643 return ((Xu
*Yu
).trace()).coefficient_tuple()[0]
1647 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1649 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1650 with the usual inner product and jordan product ``x*y =
1651 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1656 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1660 This multiplication table can be verified by hand::
1662 sage: J = JordanSpinEJA(4)
1663 sage: e0,e1,e2,e3 = J.gens()
1679 We can change the generator prefix::
1681 sage: JordanSpinEJA(2, prefix='B').gens()
1684 Our inner product satisfies the Jordan axiom::
1686 sage: set_random_seed()
1687 sage: J = JordanSpinEJA.random_instance()
1688 sage: x = J.random_element()
1689 sage: y = J.random_element()
1690 sage: z = J.random_element()
1691 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1695 def __init__(self
, n
, field
=QQ
, **kwargs
):
1696 V
= VectorSpace(field
, n
)
1697 mult_table
= [[V
.zero() for j
in range(n
)] for i
in range(n
)]
1707 z0
= x
.inner_product(y
)
1708 zbar
= y0
*xbar
+ x0
*ybar
1709 z
= V([z0
] + zbar
.list())
1710 mult_table
[i
][j
] = z
1712 # The rank of the spin algebra is two, unless we're in a
1713 # one-dimensional ambient space (because the rank is bounded by
1714 # the ambient dimension).
1715 fdeja
= super(JordanSpinEJA
, self
)
1716 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1718 def inner_product(self
, x
, y
):
1720 Faster to reimplement than to use natural representations.
1724 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1728 Ensure that this is the usual inner product for the algebras
1731 sage: set_random_seed()
1732 sage: J = JordanSpinEJA.random_instance()
1733 sage: x = J.random_element()
1734 sage: y = J.random_element()
1735 sage: X = x.natural_representation()
1736 sage: Y = y.natural_representation()
1737 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1741 return x
.to_vector().inner_product(y
.to_vector())