]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/eja/eja_algebra.py
7436ed36ce676dd9bbf0eda9922ac2aecacc6e6e
[sage.d.git] / mjo / eja / eja_algebra.py
1 """
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.
6 """
7
8 from itertools import repeat
9
10 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
11 from sage.categories.magmatic_algebras import MagmaticAlgebras
12 from sage.combinat.free_module import CombinatorialFreeModule
13 from sage.matrix.constructor import matrix
14 from sage.matrix.matrix_space import MatrixSpace
15 from sage.misc.cachefunc import cached_method
16 from sage.misc.lazy_import import lazy_import
17 from sage.misc.prandom import choice
18 from sage.misc.table import table
19 from sage.modules.free_module import FreeModule, VectorSpace
20 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
21 PolynomialRing,
22 QuadraticField)
23 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
24 lazy_import('mjo.eja.eja_subalgebra',
25 'FiniteDimensionalEuclideanJordanSubalgebra')
26 from mjo.eja.eja_utils import _mat2vec
27
28 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
29
30 def _coerce_map_from_base_ring(self):
31 """
32 Disable the map from the base ring into the algebra.
33
34 Performing a nonsense conversion like this automatically
35 is counterpedagogical. The fallback is to try the usual
36 element constructor, which should also fail.
37
38 SETUP::
39
40 sage: from mjo.eja.eja_algebra import random_eja
41
42 TESTS::
43
44 sage: set_random_seed()
45 sage: J = random_eja()
46 sage: J(1)
47 Traceback (most recent call last):
48 ...
49 ValueError: not a naturally-represented algebra element
50
51 """
52 return None
53
54 def __init__(self,
55 field,
56 mult_table,
57 rank,
58 prefix='e',
59 category=None,
60 natural_basis=None,
61 check=True):
62 """
63 SETUP::
64
65 sage: from mjo.eja.eja_algebra import (JordanSpinEJA, random_eja)
66
67 EXAMPLES:
68
69 By definition, Jordan multiplication commutes::
70
71 sage: set_random_seed()
72 sage: J = random_eja()
73 sage: x,y = J.random_elements(2)
74 sage: x*y == y*x
75 True
76
77 TESTS:
78
79 The ``field`` we're given must be real::
80
81 sage: JordanSpinEJA(2,QQbar)
82 Traceback (most recent call last):
83 ...
84 ValueError: field is not real
85
86 """
87 if check:
88 if not field.is_subring(RR):
89 # Note: this does return true for the real algebraic
90 # field, and any quadratic field where we've specified
91 # a real embedding.
92 raise ValueError('field is not real')
93
94 self._rank = rank
95 self._natural_basis = natural_basis
96
97 if category is None:
98 category = MagmaticAlgebras(field).FiniteDimensional()
99 category = category.WithBasis().Unital()
100
101 fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
102 fda.__init__(field,
103 range(len(mult_table)),
104 prefix=prefix,
105 category=category)
106 self.print_options(bracket='')
107
108 # The multiplication table we're given is necessarily in terms
109 # of vectors, because we don't have an algebra yet for
110 # anything to be an element of. However, it's faster in the
111 # long run to have the multiplication table be in terms of
112 # algebra elements. We do this after calling the superclass
113 # constructor so that from_vector() knows what to do.
114 self._multiplication_table = [
115 list(map(lambda x: self.from_vector(x), ls))
116 for ls in mult_table
117 ]
118
119
120 def _element_constructor_(self, elt):
121 """
122 Construct an element of this algebra from its natural
123 representation.
124
125 This gets called only after the parent element _call_ method
126 fails to find a coercion for the argument.
127
128 SETUP::
129
130 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
131 ....: HadamardEJA,
132 ....: RealSymmetricEJA)
133
134 EXAMPLES:
135
136 The identity in `S^n` is converted to the identity in the EJA::
137
138 sage: J = RealSymmetricEJA(3)
139 sage: I = matrix.identity(QQ,3)
140 sage: J(I) == J.one()
141 True
142
143 This skew-symmetric matrix can't be represented in the EJA::
144
145 sage: J = RealSymmetricEJA(3)
146 sage: A = matrix(QQ,3, lambda i,j: i-j)
147 sage: J(A)
148 Traceback (most recent call last):
149 ...
150 ArithmeticError: vector is not in free module
151
152 TESTS:
153
154 Ensure that we can convert any element of the two non-matrix
155 simple algebras (whose natural representations are their usual
156 vector representations) back and forth faithfully::
157
158 sage: set_random_seed()
159 sage: J = HadamardEJA.random_instance()
160 sage: x = J.random_element()
161 sage: J(x.to_vector().column()) == x
162 True
163 sage: J = JordanSpinEJA.random_instance()
164 sage: x = J.random_element()
165 sage: J(x.to_vector().column()) == x
166 True
167
168 """
169 msg = "not a naturally-represented algebra element"
170 if elt == 0:
171 # The superclass implementation of random_element()
172 # needs to be able to coerce "0" into the algebra.
173 return self.zero()
174 elif elt in self.base_ring():
175 # Ensure that no base ring -> algebra coercion is performed
176 # by this method. There's some stupidity in sage that would
177 # otherwise propagate to this method; for example, sage thinks
178 # that the integer 3 belongs to the space of 2-by-2 matrices.
179 raise ValueError(msg)
180
181 natural_basis = self.natural_basis()
182 basis_space = natural_basis[0].matrix_space()
183 if elt not in basis_space:
184 raise ValueError(msg)
185
186 # Thanks for nothing! Matrix spaces aren't vector spaces in
187 # Sage, so we have to figure out its natural-basis coordinates
188 # ourselves. We use the basis space's ring instead of the
189 # element's ring because the basis space might be an algebraic
190 # closure whereas the base ring of the 3-by-3 identity matrix
191 # could be QQ instead of QQbar.
192 V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols())
193 W = V.span_of_basis( _mat2vec(s) for s in natural_basis )
194 coords = W.coordinate_vector(_mat2vec(elt))
195 return self.from_vector(coords)
196
197
198 def _repr_(self):
199 """
200 Return a string representation of ``self``.
201
202 SETUP::
203
204 sage: from mjo.eja.eja_algebra import JordanSpinEJA
205
206 TESTS:
207
208 Ensure that it says what we think it says::
209
210 sage: JordanSpinEJA(2, field=AA)
211 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
212 sage: JordanSpinEJA(3, field=RDF)
213 Euclidean Jordan algebra of dimension 3 over Real Double Field
214
215 """
216 fmt = "Euclidean Jordan algebra of dimension {} over {}"
217 return fmt.format(self.dimension(), self.base_ring())
218
219 def product_on_basis(self, i, j):
220 return self._multiplication_table[i][j]
221
222 def _a_regular_element(self):
223 """
224 Guess a regular element. Needed to compute the basis for our
225 characteristic polynomial coefficients.
226
227 SETUP::
228
229 sage: from mjo.eja.eja_algebra import random_eja
230
231 TESTS:
232
233 Ensure that this hacky method succeeds for every algebra that we
234 know how to construct::
235
236 sage: set_random_seed()
237 sage: J = random_eja()
238 sage: J._a_regular_element().is_regular()
239 True
240
241 """
242 gs = self.gens()
243 z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
244 if not z.is_regular():
245 raise ValueError("don't know a regular element")
246 return z
247
248
249 @cached_method
250 def _charpoly_basis_space(self):
251 """
252 Return the vector space spanned by the basis used in our
253 characteristic polynomial coefficients. This is used not only to
254 compute those coefficients, but also any time we need to
255 evaluate the coefficients (like when we compute the trace or
256 determinant).
257 """
258 z = self._a_regular_element()
259 # Don't use the parent vector space directly here in case this
260 # happens to be a subalgebra. In that case, we would be e.g.
261 # two-dimensional but span_of_basis() would expect three
262 # coordinates.
263 V = VectorSpace(self.base_ring(), self.vector_space().dimension())
264 basis = [ (z**k).to_vector() for k in range(self.rank()) ]
265 V1 = V.span_of_basis( basis )
266 b = (V1.basis() + V1.complement().basis())
267 return V.span_of_basis(b)
268
269
270
271 @cached_method
272 def _charpoly_coeff(self, i):
273 """
274 Return the coefficient polynomial "a_{i}" of this algebra's
275 general characteristic polynomial.
276
277 Having this be a separate cached method lets us compute and
278 store the trace/determinant (a_{r-1} and a_{0} respectively)
279 separate from the entire characteristic polynomial.
280 """
281 (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
282 R = A_of_x.base_ring()
283
284 if i == self.rank():
285 return R.one()
286 if i > self.rank():
287 # Guaranteed by theory
288 return R.zero()
289
290 # Danger: the in-place modification is done for performance
291 # reasons (reconstructing a matrix with huge polynomial
292 # entries is slow), but I don't know how cached_method works,
293 # so it's highly possible that we're modifying some global
294 # list variable by reference, here. In other words, you
295 # probably shouldn't call this method twice on the same
296 # algebra, at the same time, in two threads
297 Ai_orig = A_of_x.column(i)
298 A_of_x.set_column(i,xr)
299 numerator = A_of_x.det()
300 A_of_x.set_column(i,Ai_orig)
301
302 # We're relying on the theory here to ensure that each a_i is
303 # indeed back in R, and the added negative signs are to make
304 # the whole charpoly expression sum to zero.
305 return R(-numerator/detA)
306
307
308 @cached_method
309 def _charpoly_matrix_system(self):
310 """
311 Compute the matrix whose entries A_ij are polynomials in
312 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
313 corresponding to `x^r` and the determinent of the matrix A =
314 [A_ij]. In other words, all of the fixed (cachable) data needed
315 to compute the coefficients of the characteristic polynomial.
316 """
317 r = self.rank()
318 n = self.dimension()
319
320 # Turn my vector space into a module so that "vectors" can
321 # have multivatiate polynomial entries.
322 names = tuple('X' + str(i) for i in range(1,n+1))
323 R = PolynomialRing(self.base_ring(), names)
324
325 # Using change_ring() on the parent's vector space doesn't work
326 # here because, in a subalgebra, that vector space has a basis
327 # and change_ring() tries to bring the basis along with it. And
328 # that doesn't work unless the new ring is a PID, which it usually
329 # won't be.
330 V = FreeModule(R,n)
331
332 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
333 # indeterminates...
334 x = V(names)
335
336 # And figure out the "left multiplication by x" matrix in
337 # that setting.
338 lmbx_cols = []
339 monomial_matrices = [ self.monomial(i).operator().matrix()
340 for i in range(n) ] # don't recompute these!
341 for k in range(n):
342 ek = self.monomial(k).to_vector()
343 lmbx_cols.append(
344 sum( x[i]*(monomial_matrices[i]*ek)
345 for i in range(n) ) )
346 Lx = matrix.column(R, lmbx_cols)
347
348 # Now we can compute powers of x "symbolically"
349 x_powers = [self.one().to_vector(), x]
350 for d in range(2, r+1):
351 x_powers.append( Lx*(x_powers[-1]) )
352
353 idmat = matrix.identity(R, n)
354
355 W = self._charpoly_basis_space()
356 W = W.change_ring(R.fraction_field())
357
358 # Starting with the standard coordinates x = (X1,X2,...,Xn)
359 # and then converting the entries to W-coordinates allows us
360 # to pass in the standard coordinates to the charpoly and get
361 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
362 # we have
363 #
364 # W.coordinates(x^2) eval'd at (standard z-coords)
365 # =
366 # W-coords of (z^2)
367 # =
368 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
369 #
370 # We want the middle equivalent thing in our matrix, but use
371 # the first equivalent thing instead so that we can pass in
372 # standard coordinates.
373 x_powers = [ W.coordinate_vector(xp) for xp in x_powers ]
374 l2 = [idmat.column(k-1) for k in range(r+1, n+1)]
375 A_of_x = matrix.column(R, n, (x_powers[:r] + l2))
376 return (A_of_x, x, x_powers[r], A_of_x.det())
377
378
379 @cached_method
380 def characteristic_polynomial(self):
381 """
382 Return a characteristic polynomial that works for all elements
383 of this algebra.
384
385 The resulting polynomial has `n+1` variables, where `n` is the
386 dimension of this algebra. The first `n` variables correspond to
387 the coordinates of an algebra element: when evaluated at the
388 coordinates of an algebra element with respect to a certain
389 basis, the result is a univariate polynomial (in the one
390 remaining variable ``t``), namely the characteristic polynomial
391 of that element.
392
393 SETUP::
394
395 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
396
397 EXAMPLES:
398
399 The characteristic polynomial in the spin algebra is given in
400 Alizadeh, Example 11.11::
401
402 sage: J = JordanSpinEJA(3)
403 sage: p = J.characteristic_polynomial(); p
404 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
405 sage: xvec = J.one().to_vector()
406 sage: p(*xvec)
407 t^2 - 2*t + 1
408
409 By definition, the characteristic polynomial is a monic
410 degree-zero polynomial in a rank-zero algebra. Note that
411 Cayley-Hamilton is indeed satisfied since the polynomial
412 ``1`` evaluates to the identity element of the algebra on
413 any argument::
414
415 sage: J = TrivialEJA()
416 sage: J.characteristic_polynomial()
417 1
418
419 """
420 r = self.rank()
421 n = self.dimension()
422
423 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n.
424 a = [ self._charpoly_coeff(i) for i in range(r+1) ]
425
426 # We go to a bit of trouble here to reorder the
427 # indeterminates, so that it's easier to evaluate the
428 # characteristic polynomial at x's coordinates and get back
429 # something in terms of t, which is what we want.
430 R = a[0].parent()
431 S = PolynomialRing(self.base_ring(),'t')
432 t = S.gen(0)
433 S = PolynomialRing(S, R.variable_names())
434 t = S(t)
435
436 return sum( a[k]*(t**k) for k in range(len(a)) )
437
438
439 def inner_product(self, x, y):
440 """
441 The inner product associated with this Euclidean Jordan algebra.
442
443 Defaults to the trace inner product, but can be overridden by
444 subclasses if they are sure that the necessary properties are
445 satisfied.
446
447 SETUP::
448
449 sage: from mjo.eja.eja_algebra import random_eja
450
451 EXAMPLES:
452
453 Our inner product is "associative," which means the following for
454 a symmetric bilinear form::
455
456 sage: set_random_seed()
457 sage: J = random_eja()
458 sage: x,y,z = J.random_elements(3)
459 sage: (x*y).inner_product(z) == y.inner_product(x*z)
460 True
461
462 """
463 X = x.natural_representation()
464 Y = y.natural_representation()
465 return self.natural_inner_product(X,Y)
466
467
468 def is_trivial(self):
469 """
470 Return whether or not this algebra is trivial.
471
472 A trivial algebra contains only the zero element.
473
474 SETUP::
475
476 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
477 ....: TrivialEJA)
478
479 EXAMPLES::
480
481 sage: J = ComplexHermitianEJA(3)
482 sage: J.is_trivial()
483 False
484
485 ::
486
487 sage: J = TrivialEJA()
488 sage: J.is_trivial()
489 True
490
491 """
492 return self.dimension() == 0
493
494
495 def multiplication_table(self):
496 """
497 Return a visual representation of this algebra's multiplication
498 table (on basis elements).
499
500 SETUP::
501
502 sage: from mjo.eja.eja_algebra import JordanSpinEJA
503
504 EXAMPLES::
505
506 sage: J = JordanSpinEJA(4)
507 sage: J.multiplication_table()
508 +----++----+----+----+----+
509 | * || e0 | e1 | e2 | e3 |
510 +====++====+====+====+====+
511 | e0 || e0 | e1 | e2 | e3 |
512 +----++----+----+----+----+
513 | e1 || e1 | e0 | 0 | 0 |
514 +----++----+----+----+----+
515 | e2 || e2 | 0 | e0 | 0 |
516 +----++----+----+----+----+
517 | e3 || e3 | 0 | 0 | e0 |
518 +----++----+----+----+----+
519
520 """
521 M = list(self._multiplication_table) # copy
522 for i in range(len(M)):
523 # M had better be "square"
524 M[i] = [self.monomial(i)] + M[i]
525 M = [["*"] + list(self.gens())] + M
526 return table(M, header_row=True, header_column=True, frame=True)
527
528
529 def natural_basis(self):
530 """
531 Return a more-natural representation of this algebra's basis.
532
533 Every finite-dimensional Euclidean Jordan Algebra is a direct
534 sum of five simple algebras, four of which comprise Hermitian
535 matrices. This method returns the original "natural" basis
536 for our underlying vector space. (Typically, the natural basis
537 is used to construct the multiplication table in the first place.)
538
539 Note that this will always return a matrix. The standard basis
540 in `R^n` will be returned as `n`-by-`1` column matrices.
541
542 SETUP::
543
544 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
545 ....: RealSymmetricEJA)
546
547 EXAMPLES::
548
549 sage: J = RealSymmetricEJA(2)
550 sage: J.basis()
551 Finite family {0: e0, 1: e1, 2: e2}
552 sage: J.natural_basis()
553 (
554 [1 0] [ 0 0.7071067811865475?] [0 0]
555 [0 0], [0.7071067811865475? 0], [0 1]
556 )
557
558 ::
559
560 sage: J = JordanSpinEJA(2)
561 sage: J.basis()
562 Finite family {0: e0, 1: e1}
563 sage: J.natural_basis()
564 (
565 [1] [0]
566 [0], [1]
567 )
568
569 """
570 if self._natural_basis is None:
571 M = self.natural_basis_space()
572 return tuple( M(b.to_vector()) for b in self.basis() )
573 else:
574 return self._natural_basis
575
576
577 def natural_basis_space(self):
578 """
579 Return the matrix space in which this algebra's natural basis
580 elements live.
581 """
582 if self._natural_basis is None or len(self._natural_basis) == 0:
583 return MatrixSpace(self.base_ring(), self.dimension(), 1)
584 else:
585 return self._natural_basis[0].matrix_space()
586
587
588 @staticmethod
589 def natural_inner_product(X,Y):
590 """
591 Compute the inner product of two naturally-represented elements.
592
593 For example in the real symmetric matrix EJA, this will compute
594 the trace inner-product of two n-by-n symmetric matrices. The
595 default should work for the real cartesian product EJA, the
596 Jordan spin EJA, and the real symmetric matrices. The others
597 will have to be overridden.
598 """
599 return (X.conjugate_transpose()*Y).trace()
600
601
602 @cached_method
603 def one(self):
604 """
605 Return the unit element of this algebra.
606
607 SETUP::
608
609 sage: from mjo.eja.eja_algebra import (HadamardEJA,
610 ....: random_eja)
611
612 EXAMPLES::
613
614 sage: J = HadamardEJA(5)
615 sage: J.one()
616 e0 + e1 + e2 + e3 + e4
617
618 TESTS:
619
620 The identity element acts like the identity::
621
622 sage: set_random_seed()
623 sage: J = random_eja()
624 sage: x = J.random_element()
625 sage: J.one()*x == x and x*J.one() == x
626 True
627
628 The matrix of the unit element's operator is the identity::
629
630 sage: set_random_seed()
631 sage: J = random_eja()
632 sage: actual = J.one().operator().matrix()
633 sage: expected = matrix.identity(J.base_ring(), J.dimension())
634 sage: actual == expected
635 True
636
637 """
638 # We can brute-force compute the matrices of the operators
639 # that correspond to the basis elements of this algebra.
640 # If some linear combination of those basis elements is the
641 # algebra identity, then the same linear combination of
642 # their matrices has to be the identity matrix.
643 #
644 # Of course, matrices aren't vectors in sage, so we have to
645 # appeal to the "long vectors" isometry.
646 oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
647
648 # Now we use basis linear algebra to find the coefficients,
649 # of the matrices-as-vectors-linear-combination, which should
650 # work for the original algebra basis too.
651 A = matrix.column(self.base_ring(), oper_vecs)
652
653 # We used the isometry on the left-hand side already, but we
654 # still need to do it for the right-hand side. Recall that we
655 # wanted something that summed to the identity matrix.
656 b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
657
658 # Now if there's an identity element in the algebra, this should work.
659 coeffs = A.solve_right(b)
660 return self.linear_combination(zip(self.gens(), coeffs))
661
662
663 def peirce_decomposition(self, c):
664 """
665 The Peirce decomposition of this algebra relative to the
666 idempotent ``c``.
667
668 In the future, this can be extended to a complete system of
669 orthogonal idempotents.
670
671 INPUT:
672
673 - ``c`` -- an idempotent of this algebra.
674
675 OUTPUT:
676
677 A triple (J0, J5, J1) containing two subalgebras and one subspace
678 of this algebra,
679
680 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
681 corresponding to the eigenvalue zero.
682
683 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
684 corresponding to the eigenvalue one-half.
685
686 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
687 corresponding to the eigenvalue one.
688
689 These are the only possible eigenspaces for that operator, and this
690 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
691 orthogonal, and are subalgebras of this algebra with the appropriate
692 restrictions.
693
694 SETUP::
695
696 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
697
698 EXAMPLES:
699
700 The canonical example comes from the symmetric matrices, which
701 decompose into diagonal and off-diagonal parts::
702
703 sage: J = RealSymmetricEJA(3)
704 sage: C = matrix(QQ, [ [1,0,0],
705 ....: [0,1,0],
706 ....: [0,0,0] ])
707 sage: c = J(C)
708 sage: J0,J5,J1 = J.peirce_decomposition(c)
709 sage: J0
710 Euclidean Jordan algebra of dimension 1...
711 sage: J5
712 Vector space of degree 6 and dimension 2...
713 sage: J1
714 Euclidean Jordan algebra of dimension 3...
715
716 TESTS:
717
718 Every algebra decomposes trivially with respect to its identity
719 element::
720
721 sage: set_random_seed()
722 sage: J = random_eja()
723 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
724 sage: J0.dimension() == 0 and J5.dimension() == 0
725 True
726 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
727 True
728
729 The identity elements in the two subalgebras are the
730 projections onto their respective subspaces of the
731 superalgebra's identity element::
732
733 sage: set_random_seed()
734 sage: J = random_eja()
735 sage: x = J.random_element()
736 sage: if not J.is_trivial():
737 ....: while x.is_nilpotent():
738 ....: x = J.random_element()
739 sage: c = x.subalgebra_idempotent()
740 sage: J0,J5,J1 = J.peirce_decomposition(c)
741 sage: J1(c) == J1.one()
742 True
743 sage: J0(J.one() - c) == J0.one()
744 True
745
746 """
747 if not c.is_idempotent():
748 raise ValueError("element is not idempotent: %s" % c)
749
750 # Default these to what they should be if they turn out to be
751 # trivial, because eigenspaces_left() won't return eigenvalues
752 # corresponding to trivial spaces (e.g. it returns only the
753 # eigenspace corresponding to lambda=1 if you take the
754 # decomposition relative to the identity element).
755 trivial = FiniteDimensionalEuclideanJordanSubalgebra(self, ())
756 J0 = trivial # eigenvalue zero
757 J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
758 J1 = trivial # eigenvalue one
759
760 for (eigval, eigspace) in c.operator().matrix().right_eigenspaces():
761 if eigval == ~(self.base_ring()(2)):
762 J5 = eigspace
763 else:
764 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
765 subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, gens)
766 if eigval == 0:
767 J0 = subalg
768 elif eigval == 1:
769 J1 = subalg
770 else:
771 raise ValueError("unexpected eigenvalue: %s" % eigval)
772
773 return (J0, J5, J1)
774
775
776 def random_elements(self, count):
777 """
778 Return ``count`` random elements as a tuple.
779
780 SETUP::
781
782 sage: from mjo.eja.eja_algebra import JordanSpinEJA
783
784 EXAMPLES::
785
786 sage: J = JordanSpinEJA(3)
787 sage: x,y,z = J.random_elements(3)
788 sage: all( [ x in J, y in J, z in J ])
789 True
790 sage: len( J.random_elements(10) ) == 10
791 True
792
793 """
794 return tuple( self.random_element() for idx in range(count) )
795
796
797 def _rank_computation(self):
798 r"""
799 Compute the rank of this algebra.
800
801 SETUP::
802
803 sage: from mjo.eja.eja_algebra import (HadamardEJA,
804 ....: JordanSpinEJA,
805 ....: RealSymmetricEJA,
806 ....: ComplexHermitianEJA,
807 ....: QuaternionHermitianEJA)
808
809 EXAMPLES::
810
811 sage: J = HadamardEJA(4)
812 sage: J._rank_computation() == J.rank()
813 True
814 sage: J = JordanSpinEJA(4)
815 sage: J._rank_computation() == J.rank()
816 True
817 sage: J = RealSymmetricEJA(3)
818 sage: J._rank_computation() == J.rank()
819 True
820 sage: J = ComplexHermitianEJA(2)
821 sage: J._rank_computation() == J.rank()
822 True
823 sage: J = QuaternionHermitianEJA(2)
824 sage: J._rank_computation() == J.rank()
825 True
826
827 """
828 n = self.dimension()
829 if n == 0:
830 return 0
831 elif n == 1:
832 return 1
833
834 var_names = [ "X" + str(z) for z in range(1,n+1) ]
835 R = PolynomialRing(self.base_ring(), var_names)
836 vars = R.gens()
837
838 def L_x_i_j(i,j):
839 # From a result in my book, these are the entries of the
840 # basis representation of L_x.
841 return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
842 for k in range(n) )
843
844 L_x = matrix(R, n, n, L_x_i_j)
845 x_powers = [ vars[k]*(L_x**k)*self.one().to_vector()
846 for k in range(n) ]
847
848 # Can assume n >= 2
849 M = matrix([x_powers[0]])
850 old_rank = 1
851
852 for d in range(1,n):
853 M = matrix(M.rows() + [x_powers[d]])
854 M.echelonize()
855 new_rank = M.rank()
856 if new_rank == old_rank:
857 return new_rank
858 else:
859 old_rank = new_rank
860
861 return n
862
863 def rank(self):
864 """
865 Return the rank of this EJA.
866
867 ALGORITHM:
868
869 The author knows of no algorithm to compute the rank of an EJA
870 where only the multiplication table is known. In lieu of one, we
871 require the rank to be specified when the algebra is created,
872 and simply pass along that number here.
873
874 SETUP::
875
876 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
877 ....: RealSymmetricEJA,
878 ....: ComplexHermitianEJA,
879 ....: QuaternionHermitianEJA,
880 ....: random_eja)
881
882 EXAMPLES:
883
884 The rank of the Jordan spin algebra is always two::
885
886 sage: JordanSpinEJA(2).rank()
887 2
888 sage: JordanSpinEJA(3).rank()
889 2
890 sage: JordanSpinEJA(4).rank()
891 2
892
893 The rank of the `n`-by-`n` Hermitian real, complex, or
894 quaternion matrices is `n`::
895
896 sage: RealSymmetricEJA(4).rank()
897 4
898 sage: ComplexHermitianEJA(3).rank()
899 3
900 sage: QuaternionHermitianEJA(2).rank()
901 2
902
903 TESTS:
904
905 Ensure that every EJA that we know how to construct has a
906 positive integer rank, unless the algebra is trivial in
907 which case its rank will be zero::
908
909 sage: set_random_seed()
910 sage: J = random_eja()
911 sage: r = J.rank()
912 sage: r in ZZ
913 True
914 sage: r > 0 or (r == 0 and J.is_trivial())
915 True
916
917 """
918 return self._rank
919
920
921 def vector_space(self):
922 """
923 Return the vector space that underlies this algebra.
924
925 SETUP::
926
927 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
928
929 EXAMPLES::
930
931 sage: J = RealSymmetricEJA(2)
932 sage: J.vector_space()
933 Vector space of dimension 3 over...
934
935 """
936 return self.zero().to_vector().parent().ambient_vector_space()
937
938
939 Element = FiniteDimensionalEuclideanJordanAlgebraElement
940
941
942 class KnownRankEJA(object):
943 """
944 A class for algebras that we actually know we can construct. The
945 main issue is that, for most of our methods to make sense, we need
946 to know the rank of our algebra. Thus we can't simply generate a
947 "random" algebra, or even check that a given basis and product
948 satisfy the axioms; because even if everything looks OK, we wouldn't
949 know the rank we need to actuallty build the thing.
950
951 Not really a subclass of FDEJA because doing that causes method
952 resolution errors, e.g.
953
954 TypeError: Error when calling the metaclass bases
955 Cannot create a consistent method resolution
956 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
957 KnownRankEJA
958
959 """
960 @staticmethod
961 def _max_test_case_size():
962 """
963 Return an integer "size" that is an upper bound on the size of
964 this algebra when it is used in a random test
965 case. Unfortunately, the term "size" is quite vague -- when
966 dealing with `R^n` under either the Hadamard or Jordan spin
967 product, the "size" refers to the dimension `n`. When dealing
968 with a matrix algebra (real symmetric or complex/quaternion
969 Hermitian), it refers to the size of the matrix, which is
970 far less than the dimension of the underlying vector space.
971
972 We default to five in this class, which is safe in `R^n`. The
973 matrix algebra subclasses (or any class where the "size" is
974 interpreted to be far less than the dimension) should override
975 with a smaller number.
976 """
977 return 5
978
979 @classmethod
980 def random_instance(cls, field=AA, **kwargs):
981 """
982 Return a random instance of this type of algebra.
983
984 Beware, this will crash for "most instances" because the
985 constructor below looks wrong.
986 """
987 if cls is TrivialEJA:
988 # The TrivialEJA class doesn't take an "n" argument because
989 # there's only one.
990 return cls(field)
991
992 n = ZZ.random_element(cls._max_test_case_size()) + 1
993 return cls(n, field, **kwargs)
994
995
996 class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
997 """
998 Return the Euclidean Jordan Algebra corresponding to the set
999 `R^n` under the Hadamard product.
1000
1001 Note: this is nothing more than the Cartesian product of ``n``
1002 copies of the spin algebra. Once Cartesian product algebras
1003 are implemented, this can go.
1004
1005 SETUP::
1006
1007 sage: from mjo.eja.eja_algebra import HadamardEJA
1008
1009 EXAMPLES:
1010
1011 This multiplication table can be verified by hand::
1012
1013 sage: J = HadamardEJA(3)
1014 sage: e0,e1,e2 = J.gens()
1015 sage: e0*e0
1016 e0
1017 sage: e0*e1
1018 0
1019 sage: e0*e2
1020 0
1021 sage: e1*e1
1022 e1
1023 sage: e1*e2
1024 0
1025 sage: e2*e2
1026 e2
1027
1028 TESTS:
1029
1030 We can change the generator prefix::
1031
1032 sage: HadamardEJA(3, prefix='r').gens()
1033 (r0, r1, r2)
1034
1035 """
1036 def __init__(self, n, field=AA, **kwargs):
1037 V = VectorSpace(field, n)
1038 mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
1039 for i in range(n) ]
1040
1041 fdeja = super(HadamardEJA, self)
1042 return fdeja.__init__(field, mult_table, rank=n, **kwargs)
1043
1044 def inner_product(self, x, y):
1045 """
1046 Faster to reimplement than to use natural representations.
1047
1048 SETUP::
1049
1050 sage: from mjo.eja.eja_algebra import HadamardEJA
1051
1052 TESTS:
1053
1054 Ensure that this is the usual inner product for the algebras
1055 over `R^n`::
1056
1057 sage: set_random_seed()
1058 sage: J = HadamardEJA.random_instance()
1059 sage: x,y = J.random_elements(2)
1060 sage: X = x.natural_representation()
1061 sage: Y = y.natural_representation()
1062 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1063 True
1064
1065 """
1066 return x.to_vector().inner_product(y.to_vector())
1067
1068
1069 def random_eja(field=AA, nontrivial=False):
1070 """
1071 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1072
1073 SETUP::
1074
1075 sage: from mjo.eja.eja_algebra import random_eja
1076
1077 TESTS::
1078
1079 sage: random_eja()
1080 Euclidean Jordan algebra of dimension...
1081
1082 """
1083 eja_classes = KnownRankEJA.__subclasses__()
1084 if nontrivial:
1085 eja_classes.remove(TrivialEJA)
1086 classname = choice(eja_classes)
1087 return classname.random_instance(field=field)
1088
1089
1090
1091
1092
1093
1094 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
1095 @staticmethod
1096 def _max_test_case_size():
1097 # Play it safe, since this will be squared and the underlying
1098 # field can have dimension 4 (quaternions) too.
1099 return 2
1100
1101 def __init__(self, field, basis, rank, normalize_basis=True, **kwargs):
1102 """
1103 Compared to the superclass constructor, we take a basis instead of
1104 a multiplication table because the latter can be computed in terms
1105 of the former when the product is known (like it is here).
1106 """
1107 # Used in this class's fast _charpoly_coeff() override.
1108 self._basis_normalizers = None
1109
1110 # We're going to loop through this a few times, so now's a good
1111 # time to ensure that it isn't a generator expression.
1112 basis = tuple(basis)
1113
1114 if rank > 1 and normalize_basis:
1115 # We'll need sqrt(2) to normalize the basis, and this
1116 # winds up in the multiplication table, so the whole
1117 # algebra needs to be over the field extension.
1118 R = PolynomialRing(field, 'z')
1119 z = R.gen()
1120 p = z**2 - 2
1121 if p.is_irreducible():
1122 field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())
1123 basis = tuple( s.change_ring(field) for s in basis )
1124 self._basis_normalizers = tuple(
1125 ~(self.natural_inner_product(s,s).sqrt()) for s in basis )
1126 basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
1127
1128 Qs = self.multiplication_table_from_matrix_basis(basis)
1129
1130 fdeja = super(MatrixEuclideanJordanAlgebra, self)
1131 return fdeja.__init__(field,
1132 Qs,
1133 rank=rank,
1134 natural_basis=basis,
1135 **kwargs)
1136
1137
1138 def _rank_computation(self):
1139 r"""
1140 Override the parent method with something that tries to compute
1141 over a faster (non-extension) field.
1142 """
1143 if self._basis_normalizers is None:
1144 # We didn't normalize, so assume that the basis we started
1145 # with had entries in a nice field.
1146 return super(MatrixEuclideanJordanAlgebra, self)._rank_computation()
1147 else:
1148 basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
1149 self._basis_normalizers) )
1150
1151 # Do this over the rationals and convert back at the end.
1152 # Only works because we know the entries of the basis are
1153 # integers.
1154 J = MatrixEuclideanJordanAlgebra(QQ,
1155 basis,
1156 self.rank(),
1157 normalize_basis=False)
1158 return J._rank_computation()
1159
1160 @cached_method
1161 def _charpoly_coeff(self, i):
1162 """
1163 Override the parent method with something that tries to compute
1164 over a faster (non-extension) field.
1165 """
1166 if self._basis_normalizers is None:
1167 # We didn't normalize, so assume that the basis we started
1168 # with had entries in a nice field.
1169 return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i)
1170 else:
1171 basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
1172 self._basis_normalizers) )
1173
1174 # Do this over the rationals and convert back at the end.
1175 J = MatrixEuclideanJordanAlgebra(QQ,
1176 basis,
1177 self.rank(),
1178 normalize_basis=False)
1179 (_,x,_,_) = J._charpoly_matrix_system()
1180 p = J._charpoly_coeff(i)
1181 # p might be missing some vars, have to substitute "optionally"
1182 pairs = zip(x.base_ring().gens(), self._basis_normalizers)
1183 substitutions = { v: v*c for (v,c) in pairs }
1184 result = p.subs(substitutions)
1185
1186 # The result of "subs" can be either a coefficient-ring
1187 # element or a polynomial. Gotta handle both cases.
1188 if result in QQ:
1189 return self.base_ring()(result)
1190 else:
1191 return result.change_ring(self.base_ring())
1192
1193
1194 @staticmethod
1195 def multiplication_table_from_matrix_basis(basis):
1196 """
1197 At least three of the five simple Euclidean Jordan algebras have the
1198 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1199 multiplication on the right is matrix multiplication. Given a basis
1200 for the underlying matrix space, this function returns a
1201 multiplication table (obtained by looping through the basis
1202 elements) for an algebra of those matrices.
1203 """
1204 # In S^2, for example, we nominally have four coordinates even
1205 # though the space is of dimension three only. The vector space V
1206 # is supposed to hold the entire long vector, and the subspace W
1207 # of V will be spanned by the vectors that arise from symmetric
1208 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1209 field = basis[0].base_ring()
1210 dimension = basis[0].nrows()
1211
1212 V = VectorSpace(field, dimension**2)
1213 W = V.span_of_basis( _mat2vec(s) for s in basis )
1214 n = len(basis)
1215 mult_table = [[W.zero() for j in range(n)] for i in range(n)]
1216 for i in range(n):
1217 for j in range(n):
1218 mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
1219 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
1220
1221 return mult_table
1222
1223
1224 @staticmethod
1225 def real_embed(M):
1226 """
1227 Embed the matrix ``M`` into a space of real matrices.
1228
1229 The matrix ``M`` can have entries in any field at the moment:
1230 the real numbers, complex numbers, or quaternions. And although
1231 they are not a field, we can probably support octonions at some
1232 point, too. This function returns a real matrix that "acts like"
1233 the original with respect to matrix multiplication; i.e.
1234
1235 real_embed(M*N) = real_embed(M)*real_embed(N)
1236
1237 """
1238 raise NotImplementedError
1239
1240
1241 @staticmethod
1242 def real_unembed(M):
1243 """
1244 The inverse of :meth:`real_embed`.
1245 """
1246 raise NotImplementedError
1247
1248
1249 @classmethod
1250 def natural_inner_product(cls,X,Y):
1251 Xu = cls.real_unembed(X)
1252 Yu = cls.real_unembed(Y)
1253 tr = (Xu*Yu).trace()
1254
1255 if tr in RLF:
1256 # It's real already.
1257 return tr
1258
1259 # Otherwise, try the thing that works for complex numbers; and
1260 # if that doesn't work, the thing that works for quaternions.
1261 try:
1262 return tr.vector()[0] # real part, imag part is index 1
1263 except AttributeError:
1264 # A quaternions doesn't have a vector() method, but does
1265 # have coefficient_tuple() method that returns the
1266 # coefficients of 1, i, j, and k -- in that order.
1267 return tr.coefficient_tuple()[0]
1268
1269
1270 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1271 @staticmethod
1272 def real_embed(M):
1273 """
1274 The identity function, for embedding real matrices into real
1275 matrices.
1276 """
1277 return M
1278
1279 @staticmethod
1280 def real_unembed(M):
1281 """
1282 The identity function, for unembedding real matrices from real
1283 matrices.
1284 """
1285 return M
1286
1287
1288 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
1289 """
1290 The rank-n simple EJA consisting of real symmetric n-by-n
1291 matrices, the usual symmetric Jordan product, and the trace inner
1292 product. It has dimension `(n^2 + n)/2` over the reals.
1293
1294 SETUP::
1295
1296 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1297
1298 EXAMPLES::
1299
1300 sage: J = RealSymmetricEJA(2)
1301 sage: e0, e1, e2 = J.gens()
1302 sage: e0*e0
1303 e0
1304 sage: e1*e1
1305 1/2*e0 + 1/2*e2
1306 sage: e2*e2
1307 e2
1308
1309 In theory, our "field" can be any subfield of the reals::
1310
1311 sage: RealSymmetricEJA(2, RDF)
1312 Euclidean Jordan algebra of dimension 3 over Real Double Field
1313 sage: RealSymmetricEJA(2, RR)
1314 Euclidean Jordan algebra of dimension 3 over Real Field with
1315 53 bits of precision
1316
1317 TESTS:
1318
1319 The dimension of this algebra is `(n^2 + n) / 2`::
1320
1321 sage: set_random_seed()
1322 sage: n_max = RealSymmetricEJA._max_test_case_size()
1323 sage: n = ZZ.random_element(1, n_max)
1324 sage: J = RealSymmetricEJA(n)
1325 sage: J.dimension() == (n^2 + n)/2
1326 True
1327
1328 The Jordan multiplication is what we think it is::
1329
1330 sage: set_random_seed()
1331 sage: J = RealSymmetricEJA.random_instance()
1332 sage: x,y = J.random_elements(2)
1333 sage: actual = (x*y).natural_representation()
1334 sage: X = x.natural_representation()
1335 sage: Y = y.natural_representation()
1336 sage: expected = (X*Y + Y*X)/2
1337 sage: actual == expected
1338 True
1339 sage: J(expected) == x*y
1340 True
1341
1342 We can change the generator prefix::
1343
1344 sage: RealSymmetricEJA(3, prefix='q').gens()
1345 (q0, q1, q2, q3, q4, q5)
1346
1347 Our natural basis is normalized with respect to the natural inner
1348 product unless we specify otherwise::
1349
1350 sage: set_random_seed()
1351 sage: J = RealSymmetricEJA.random_instance()
1352 sage: all( b.norm() == 1 for b in J.gens() )
1353 True
1354
1355 Since our natural basis is normalized with respect to the natural
1356 inner product, and since we know that this algebra is an EJA, any
1357 left-multiplication operator's matrix will be symmetric because
1358 natural->EJA basis representation is an isometry and within the EJA
1359 the operator is self-adjoint by the Jordan axiom::
1360
1361 sage: set_random_seed()
1362 sage: x = RealSymmetricEJA.random_instance().random_element()
1363 sage: x.operator().matrix().is_symmetric()
1364 True
1365
1366 """
1367 @classmethod
1368 def _denormalized_basis(cls, n, field):
1369 """
1370 Return a basis for the space of real symmetric n-by-n matrices.
1371
1372 SETUP::
1373
1374 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1375
1376 TESTS::
1377
1378 sage: set_random_seed()
1379 sage: n = ZZ.random_element(1,5)
1380 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1381 sage: all( M.is_symmetric() for M in B)
1382 True
1383
1384 """
1385 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1386 # coordinates.
1387 S = []
1388 for i in range(n):
1389 for j in range(i+1):
1390 Eij = matrix(field, n, lambda k,l: k==i and l==j)
1391 if i == j:
1392 Sij = Eij
1393 else:
1394 Sij = Eij + Eij.transpose()
1395 S.append(Sij)
1396 return S
1397
1398
1399 @staticmethod
1400 def _max_test_case_size():
1401 return 4 # Dimension 10
1402
1403
1404 def __init__(self, n, field=AA, **kwargs):
1405 basis = self._denormalized_basis(n, field)
1406 super(RealSymmetricEJA, self).__init__(field, basis, n, **kwargs)
1407
1408
1409 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1410 @staticmethod
1411 def real_embed(M):
1412 """
1413 Embed the n-by-n complex matrix ``M`` into the space of real
1414 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1415 bi` to the block matrix ``[[a,b],[-b,a]]``.
1416
1417 SETUP::
1418
1419 sage: from mjo.eja.eja_algebra import \
1420 ....: ComplexMatrixEuclideanJordanAlgebra
1421
1422 EXAMPLES::
1423
1424 sage: F = QuadraticField(-1, 'I')
1425 sage: x1 = F(4 - 2*i)
1426 sage: x2 = F(1 + 2*i)
1427 sage: x3 = F(-i)
1428 sage: x4 = F(6)
1429 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1430 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1431 [ 4 -2| 1 2]
1432 [ 2 4|-2 1]
1433 [-----+-----]
1434 [ 0 -1| 6 0]
1435 [ 1 0| 0 6]
1436
1437 TESTS:
1438
1439 Embedding is a homomorphism (isomorphism, in fact)::
1440
1441 sage: set_random_seed()
1442 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1443 sage: n = ZZ.random_element(n_max)
1444 sage: F = QuadraticField(-1, 'I')
1445 sage: X = random_matrix(F, n)
1446 sage: Y = random_matrix(F, n)
1447 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1448 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1449 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1450 sage: Xe*Ye == XYe
1451 True
1452
1453 """
1454 n = M.nrows()
1455 if M.ncols() != n:
1456 raise ValueError("the matrix 'M' must be square")
1457
1458 # We don't need any adjoined elements...
1459 field = M.base_ring().base_ring()
1460
1461 blocks = []
1462 for z in M.list():
1463 a = z.list()[0] # real part, I guess
1464 b = z.list()[1] # imag part, I guess
1465 blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1466
1467 return matrix.block(field, n, blocks)
1468
1469
1470 @staticmethod
1471 def real_unembed(M):
1472 """
1473 The inverse of _embed_complex_matrix().
1474
1475 SETUP::
1476
1477 sage: from mjo.eja.eja_algebra import \
1478 ....: ComplexMatrixEuclideanJordanAlgebra
1479
1480 EXAMPLES::
1481
1482 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1483 ....: [-2, 1, -4, 3],
1484 ....: [ 9, 10, 11, 12],
1485 ....: [-10, 9, -12, 11] ])
1486 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1487 [ 2*I + 1 4*I + 3]
1488 [ 10*I + 9 12*I + 11]
1489
1490 TESTS:
1491
1492 Unembedding is the inverse of embedding::
1493
1494 sage: set_random_seed()
1495 sage: F = QuadraticField(-1, 'I')
1496 sage: M = random_matrix(F, 3)
1497 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1498 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1499 True
1500
1501 """
1502 n = ZZ(M.nrows())
1503 if M.ncols() != n:
1504 raise ValueError("the matrix 'M' must be square")
1505 if not n.mod(2).is_zero():
1506 raise ValueError("the matrix 'M' must be a complex embedding")
1507
1508 # If "M" was normalized, its base ring might have roots
1509 # adjoined and they can stick around after unembedding.
1510 field = M.base_ring()
1511 R = PolynomialRing(field, 'z')
1512 z = R.gen()
1513 if field is AA:
1514 # Sage doesn't know how to embed AA into QQbar, i.e. how
1515 # to adjoin sqrt(-1) to AA.
1516 F = QQbar
1517 else:
1518 F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
1519 i = F.gen()
1520
1521 # Go top-left to bottom-right (reading order), converting every
1522 # 2-by-2 block we see to a single complex element.
1523 elements = []
1524 for k in range(n/2):
1525 for j in range(n/2):
1526 submat = M[2*k:2*k+2,2*j:2*j+2]
1527 if submat[0,0] != submat[1,1]:
1528 raise ValueError('bad on-diagonal submatrix')
1529 if submat[0,1] != -submat[1,0]:
1530 raise ValueError('bad off-diagonal submatrix')
1531 z = submat[0,0] + submat[0,1]*i
1532 elements.append(z)
1533
1534 return matrix(F, n/2, elements)
1535
1536
1537 @classmethod
1538 def natural_inner_product(cls,X,Y):
1539 """
1540 Compute a natural inner product in this algebra directly from
1541 its real embedding.
1542
1543 SETUP::
1544
1545 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1546
1547 TESTS:
1548
1549 This gives the same answer as the slow, default method implemented
1550 in :class:`MatrixEuclideanJordanAlgebra`::
1551
1552 sage: set_random_seed()
1553 sage: J = ComplexHermitianEJA.random_instance()
1554 sage: x,y = J.random_elements(2)
1555 sage: Xe = x.natural_representation()
1556 sage: Ye = y.natural_representation()
1557 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1558 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1559 sage: expected = (X*Y).trace().real()
1560 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1561 sage: actual == expected
1562 True
1563
1564 """
1565 return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
1566
1567
1568 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
1569 """
1570 The rank-n simple EJA consisting of complex Hermitian n-by-n
1571 matrices over the real numbers, the usual symmetric Jordan product,
1572 and the real-part-of-trace inner product. It has dimension `n^2` over
1573 the reals.
1574
1575 SETUP::
1576
1577 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1578
1579 EXAMPLES:
1580
1581 In theory, our "field" can be any subfield of the reals::
1582
1583 sage: ComplexHermitianEJA(2, RDF)
1584 Euclidean Jordan algebra of dimension 4 over Real Double Field
1585 sage: ComplexHermitianEJA(2, RR)
1586 Euclidean Jordan algebra of dimension 4 over Real Field with
1587 53 bits of precision
1588
1589 TESTS:
1590
1591 The dimension of this algebra is `n^2`::
1592
1593 sage: set_random_seed()
1594 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1595 sage: n = ZZ.random_element(1, n_max)
1596 sage: J = ComplexHermitianEJA(n)
1597 sage: J.dimension() == n^2
1598 True
1599
1600 The Jordan multiplication is what we think it is::
1601
1602 sage: set_random_seed()
1603 sage: J = ComplexHermitianEJA.random_instance()
1604 sage: x,y = J.random_elements(2)
1605 sage: actual = (x*y).natural_representation()
1606 sage: X = x.natural_representation()
1607 sage: Y = y.natural_representation()
1608 sage: expected = (X*Y + Y*X)/2
1609 sage: actual == expected
1610 True
1611 sage: J(expected) == x*y
1612 True
1613
1614 We can change the generator prefix::
1615
1616 sage: ComplexHermitianEJA(2, prefix='z').gens()
1617 (z0, z1, z2, z3)
1618
1619 Our natural basis is normalized with respect to the natural inner
1620 product unless we specify otherwise::
1621
1622 sage: set_random_seed()
1623 sage: J = ComplexHermitianEJA.random_instance()
1624 sage: all( b.norm() == 1 for b in J.gens() )
1625 True
1626
1627 Since our natural basis is normalized with respect to the natural
1628 inner product, and since we know that this algebra is an EJA, any
1629 left-multiplication operator's matrix will be symmetric because
1630 natural->EJA basis representation is an isometry and within the EJA
1631 the operator is self-adjoint by the Jordan axiom::
1632
1633 sage: set_random_seed()
1634 sage: x = ComplexHermitianEJA.random_instance().random_element()
1635 sage: x.operator().matrix().is_symmetric()
1636 True
1637
1638 """
1639
1640 @classmethod
1641 def _denormalized_basis(cls, n, field):
1642 """
1643 Returns a basis for the space of complex Hermitian n-by-n matrices.
1644
1645 Why do we embed these? Basically, because all of numerical linear
1646 algebra assumes that you're working with vectors consisting of `n`
1647 entries from a field and scalars from the same field. There's no way
1648 to tell SageMath that (for example) the vectors contain complex
1649 numbers, while the scalar field is real.
1650
1651 SETUP::
1652
1653 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1654
1655 TESTS::
1656
1657 sage: set_random_seed()
1658 sage: n = ZZ.random_element(1,5)
1659 sage: field = QuadraticField(2, 'sqrt2')
1660 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1661 sage: all( M.is_symmetric() for M in B)
1662 True
1663
1664 """
1665 R = PolynomialRing(field, 'z')
1666 z = R.gen()
1667 F = field.extension(z**2 + 1, 'I')
1668 I = F.gen()
1669
1670 # This is like the symmetric case, but we need to be careful:
1671 #
1672 # * We want conjugate-symmetry, not just symmetry.
1673 # * The diagonal will (as a result) be real.
1674 #
1675 S = []
1676 for i in range(n):
1677 for j in range(i+1):
1678 Eij = matrix(F, n, lambda k,l: k==i and l==j)
1679 if i == j:
1680 Sij = cls.real_embed(Eij)
1681 S.append(Sij)
1682 else:
1683 # The second one has a minus because it's conjugated.
1684 Sij_real = cls.real_embed(Eij + Eij.transpose())
1685 S.append(Sij_real)
1686 Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
1687 S.append(Sij_imag)
1688
1689 # Since we embedded these, we can drop back to the "field" that we
1690 # started with instead of the complex extension "F".
1691 return ( s.change_ring(field) for s in S )
1692
1693
1694 def __init__(self, n, field=AA, **kwargs):
1695 basis = self._denormalized_basis(n,field)
1696 super(ComplexHermitianEJA,self).__init__(field, basis, n, **kwargs)
1697
1698
1699 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1700 @staticmethod
1701 def real_embed(M):
1702 """
1703 Embed the n-by-n quaternion matrix ``M`` into the space of real
1704 matrices of size 4n-by-4n by first sending each quaternion entry `z
1705 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1706 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1707 matrix.
1708
1709 SETUP::
1710
1711 sage: from mjo.eja.eja_algebra import \
1712 ....: QuaternionMatrixEuclideanJordanAlgebra
1713
1714 EXAMPLES::
1715
1716 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1717 sage: i,j,k = Q.gens()
1718 sage: x = 1 + 2*i + 3*j + 4*k
1719 sage: M = matrix(Q, 1, [[x]])
1720 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1721 [ 1 2 3 4]
1722 [-2 1 -4 3]
1723 [-3 4 1 -2]
1724 [-4 -3 2 1]
1725
1726 Embedding is a homomorphism (isomorphism, in fact)::
1727
1728 sage: set_random_seed()
1729 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1730 sage: n = ZZ.random_element(n_max)
1731 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1732 sage: X = random_matrix(Q, n)
1733 sage: Y = random_matrix(Q, n)
1734 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1735 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1736 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1737 sage: Xe*Ye == XYe
1738 True
1739
1740 """
1741 quaternions = M.base_ring()
1742 n = M.nrows()
1743 if M.ncols() != n:
1744 raise ValueError("the matrix 'M' must be square")
1745
1746 F = QuadraticField(-1, 'I')
1747 i = F.gen()
1748
1749 blocks = []
1750 for z in M.list():
1751 t = z.coefficient_tuple()
1752 a = t[0]
1753 b = t[1]
1754 c = t[2]
1755 d = t[3]
1756 cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
1757 [-c + d*i, a - b*i]])
1758 realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM)
1759 blocks.append(realM)
1760
1761 # We should have real entries by now, so use the realest field
1762 # we've got for the return value.
1763 return matrix.block(quaternions.base_ring(), n, blocks)
1764
1765
1766
1767 @staticmethod
1768 def real_unembed(M):
1769 """
1770 The inverse of _embed_quaternion_matrix().
1771
1772 SETUP::
1773
1774 sage: from mjo.eja.eja_algebra import \
1775 ....: QuaternionMatrixEuclideanJordanAlgebra
1776
1777 EXAMPLES::
1778
1779 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1780 ....: [-2, 1, -4, 3],
1781 ....: [-3, 4, 1, -2],
1782 ....: [-4, -3, 2, 1]])
1783 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1784 [1 + 2*i + 3*j + 4*k]
1785
1786 TESTS:
1787
1788 Unembedding is the inverse of embedding::
1789
1790 sage: set_random_seed()
1791 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1792 sage: M = random_matrix(Q, 3)
1793 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1794 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1795 True
1796
1797 """
1798 n = ZZ(M.nrows())
1799 if M.ncols() != n:
1800 raise ValueError("the matrix 'M' must be square")
1801 if not n.mod(4).is_zero():
1802 raise ValueError("the matrix 'M' must be a quaternion embedding")
1803
1804 # Use the base ring of the matrix to ensure that its entries can be
1805 # multiplied by elements of the quaternion algebra.
1806 field = M.base_ring()
1807 Q = QuaternionAlgebra(field,-1,-1)
1808 i,j,k = Q.gens()
1809
1810 # Go top-left to bottom-right (reading order), converting every
1811 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1812 # quaternion block.
1813 elements = []
1814 for l in range(n/4):
1815 for m in range(n/4):
1816 submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
1817 M[4*l:4*l+4,4*m:4*m+4] )
1818 if submat[0,0] != submat[1,1].conjugate():
1819 raise ValueError('bad on-diagonal submatrix')
1820 if submat[0,1] != -submat[1,0].conjugate():
1821 raise ValueError('bad off-diagonal submatrix')
1822 z = submat[0,0].real()
1823 z += submat[0,0].imag()*i
1824 z += submat[0,1].real()*j
1825 z += submat[0,1].imag()*k
1826 elements.append(z)
1827
1828 return matrix(Q, n/4, elements)
1829
1830
1831 @classmethod
1832 def natural_inner_product(cls,X,Y):
1833 """
1834 Compute a natural inner product in this algebra directly from
1835 its real embedding.
1836
1837 SETUP::
1838
1839 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1840
1841 TESTS:
1842
1843 This gives the same answer as the slow, default method implemented
1844 in :class:`MatrixEuclideanJordanAlgebra`::
1845
1846 sage: set_random_seed()
1847 sage: J = QuaternionHermitianEJA.random_instance()
1848 sage: x,y = J.random_elements(2)
1849 sage: Xe = x.natural_representation()
1850 sage: Ye = y.natural_representation()
1851 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1852 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1853 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1854 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1855 sage: actual == expected
1856 True
1857
1858 """
1859 return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4
1860
1861
1862 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
1863 KnownRankEJA):
1864 """
1865 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1866 matrices, the usual symmetric Jordan product, and the
1867 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1868 the reals.
1869
1870 SETUP::
1871
1872 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1873
1874 EXAMPLES:
1875
1876 In theory, our "field" can be any subfield of the reals::
1877
1878 sage: QuaternionHermitianEJA(2, RDF)
1879 Euclidean Jordan algebra of dimension 6 over Real Double Field
1880 sage: QuaternionHermitianEJA(2, RR)
1881 Euclidean Jordan algebra of dimension 6 over Real Field with
1882 53 bits of precision
1883
1884 TESTS:
1885
1886 The dimension of this algebra is `2*n^2 - n`::
1887
1888 sage: set_random_seed()
1889 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1890 sage: n = ZZ.random_element(1, n_max)
1891 sage: J = QuaternionHermitianEJA(n)
1892 sage: J.dimension() == 2*(n^2) - n
1893 True
1894
1895 The Jordan multiplication is what we think it is::
1896
1897 sage: set_random_seed()
1898 sage: J = QuaternionHermitianEJA.random_instance()
1899 sage: x,y = J.random_elements(2)
1900 sage: actual = (x*y).natural_representation()
1901 sage: X = x.natural_representation()
1902 sage: Y = y.natural_representation()
1903 sage: expected = (X*Y + Y*X)/2
1904 sage: actual == expected
1905 True
1906 sage: J(expected) == x*y
1907 True
1908
1909 We can change the generator prefix::
1910
1911 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1912 (a0, a1, a2, a3, a4, a5)
1913
1914 Our natural basis is normalized with respect to the natural inner
1915 product unless we specify otherwise::
1916
1917 sage: set_random_seed()
1918 sage: J = QuaternionHermitianEJA.random_instance()
1919 sage: all( b.norm() == 1 for b in J.gens() )
1920 True
1921
1922 Since our natural basis is normalized with respect to the natural
1923 inner product, and since we know that this algebra is an EJA, any
1924 left-multiplication operator's matrix will be symmetric because
1925 natural->EJA basis representation is an isometry and within the EJA
1926 the operator is self-adjoint by the Jordan axiom::
1927
1928 sage: set_random_seed()
1929 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1930 sage: x.operator().matrix().is_symmetric()
1931 True
1932
1933 """
1934 @classmethod
1935 def _denormalized_basis(cls, n, field):
1936 """
1937 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1938
1939 Why do we embed these? Basically, because all of numerical
1940 linear algebra assumes that you're working with vectors consisting
1941 of `n` entries from a field and scalars from the same field. There's
1942 no way to tell SageMath that (for example) the vectors contain
1943 complex numbers, while the scalar field is real.
1944
1945 SETUP::
1946
1947 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1948
1949 TESTS::
1950
1951 sage: set_random_seed()
1952 sage: n = ZZ.random_element(1,5)
1953 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1954 sage: all( M.is_symmetric() for M in B )
1955 True
1956
1957 """
1958 Q = QuaternionAlgebra(QQ,-1,-1)
1959 I,J,K = Q.gens()
1960
1961 # This is like the symmetric case, but we need to be careful:
1962 #
1963 # * We want conjugate-symmetry, not just symmetry.
1964 # * The diagonal will (as a result) be real.
1965 #
1966 S = []
1967 for i in range(n):
1968 for j in range(i+1):
1969 Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1970 if i == j:
1971 Sij = cls.real_embed(Eij)
1972 S.append(Sij)
1973 else:
1974 # The second, third, and fourth ones have a minus
1975 # because they're conjugated.
1976 Sij_real = cls.real_embed(Eij + Eij.transpose())
1977 S.append(Sij_real)
1978 Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
1979 S.append(Sij_I)
1980 Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
1981 S.append(Sij_J)
1982 Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
1983 S.append(Sij_K)
1984
1985 # Since we embedded these, we can drop back to the "field" that we
1986 # started with instead of the quaternion algebra "Q".
1987 return ( s.change_ring(field) for s in S )
1988
1989
1990 def __init__(self, n, field=AA, **kwargs):
1991 basis = self._denormalized_basis(n,field)
1992 super(QuaternionHermitianEJA,self).__init__(field, basis, n, **kwargs)
1993
1994
1995 class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
1996 r"""
1997 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1998 with the half-trace inner product and jordan product ``x*y =
1999 (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
2000 symmetric positive-definite "bilinear form" matrix. It has
2001 dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
2002 when ``B`` is the identity matrix of order ``n-1``.
2003
2004 SETUP::
2005
2006 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2007 ....: JordanSpinEJA)
2008
2009 EXAMPLES:
2010
2011 When no bilinear form is specified, the identity matrix is used,
2012 and the resulting algebra is the Jordan spin algebra::
2013
2014 sage: J0 = BilinearFormEJA(3)
2015 sage: J1 = JordanSpinEJA(3)
2016 sage: J0.multiplication_table() == J0.multiplication_table()
2017 True
2018
2019 TESTS:
2020
2021 We can create a zero-dimensional algebra::
2022
2023 sage: J = BilinearFormEJA(0)
2024 sage: J.basis()
2025 Finite family {}
2026
2027 We can check the multiplication condition given in the Jordan, von
2028 Neumann, and Wigner paper (and also discussed on my "On the
2029 symmetry..." paper). Note that this relies heavily on the standard
2030 choice of basis, as does anything utilizing the bilinear form matrix::
2031
2032 sage: set_random_seed()
2033 sage: n = ZZ.random_element(5)
2034 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2035 sage: B = M.transpose()*M
2036 sage: J = BilinearFormEJA(n, B=B)
2037 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2038 sage: V = J.vector_space()
2039 sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
2040 ....: for ei in eis ]
2041 sage: actual = [ sis[i]*sis[j]
2042 ....: for i in range(n-1)
2043 ....: for j in range(n-1) ]
2044 sage: expected = [ J.one() if i == j else J.zero()
2045 ....: for i in range(n-1)
2046 ....: for j in range(n-1) ]
2047 sage: actual == expected
2048 True
2049 """
2050 def __init__(self, n, field=AA, B=None, **kwargs):
2051 if B is None:
2052 self._B = matrix.identity(field, max(0,n-1))
2053 else:
2054 self._B = B
2055
2056 V = VectorSpace(field, n)
2057 mult_table = [[V.zero() for j in range(n)] for i in range(n)]
2058 for i in range(n):
2059 for j in range(n):
2060 x = V.gen(i)
2061 y = V.gen(j)
2062 x0 = x[0]
2063 xbar = x[1:]
2064 y0 = y[0]
2065 ybar = y[1:]
2066 z0 = x0*y0 + (self._B*xbar).inner_product(ybar)
2067 zbar = y0*xbar + x0*ybar
2068 z = V([z0] + zbar.list())
2069 mult_table[i][j] = z
2070
2071 # The rank of this algebra is two, unless we're in a
2072 # one-dimensional ambient space (because the rank is bounded
2073 # by the ambient dimension).
2074 fdeja = super(BilinearFormEJA, self)
2075 return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
2076
2077 def inner_product(self, x, y):
2078 r"""
2079 Half of the trace inner product.
2080
2081 This is defined so that the special case of the Jordan spin
2082 algebra gets the usual inner product.
2083
2084 SETUP::
2085
2086 sage: from mjo.eja.eja_algebra import BilinearFormEJA
2087
2088 TESTS:
2089
2090 Ensure that this is one-half of the trace inner-product when
2091 the algebra isn't just the reals (when ``n`` isn't one). This
2092 is in Faraut and Koranyi, and also my "On the symmetry..."
2093 paper::
2094
2095 sage: set_random_seed()
2096 sage: n = ZZ.random_element(2,5)
2097 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2098 sage: B = M.transpose()*M
2099 sage: J = BilinearFormEJA(n, B=B)
2100 sage: x = J.random_element()
2101 sage: y = J.random_element()
2102 sage: x.inner_product(y) == (x*y).trace()/2
2103 True
2104
2105 """
2106 xvec = x.to_vector()
2107 xbar = xvec[1:]
2108 yvec = y.to_vector()
2109 ybar = yvec[1:]
2110 return x[0]*y[0] + (self._B*xbar).inner_product(ybar)
2111
2112
2113 class JordanSpinEJA(BilinearFormEJA):
2114 """
2115 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2116 with the usual inner product and jordan product ``x*y =
2117 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2118 the reals.
2119
2120 SETUP::
2121
2122 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2123
2124 EXAMPLES:
2125
2126 This multiplication table can be verified by hand::
2127
2128 sage: J = JordanSpinEJA(4)
2129 sage: e0,e1,e2,e3 = J.gens()
2130 sage: e0*e0
2131 e0
2132 sage: e0*e1
2133 e1
2134 sage: e0*e2
2135 e2
2136 sage: e0*e3
2137 e3
2138 sage: e1*e2
2139 0
2140 sage: e1*e3
2141 0
2142 sage: e2*e3
2143 0
2144
2145 We can change the generator prefix::
2146
2147 sage: JordanSpinEJA(2, prefix='B').gens()
2148 (B0, B1)
2149
2150 TESTS:
2151
2152 Ensure that we have the usual inner product on `R^n`::
2153
2154 sage: set_random_seed()
2155 sage: J = JordanSpinEJA.random_instance()
2156 sage: x,y = J.random_elements(2)
2157 sage: X = x.natural_representation()
2158 sage: Y = y.natural_representation()
2159 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
2160 True
2161
2162 """
2163 def __init__(self, n, field=AA, **kwargs):
2164 # This is a special case of the BilinearFormEJA with the identity
2165 # matrix as its bilinear form.
2166 return super(JordanSpinEJA, self).__init__(n, field, **kwargs)
2167
2168
2169 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
2170 """
2171 The trivial Euclidean Jordan algebra consisting of only a zero element.
2172
2173 SETUP::
2174
2175 sage: from mjo.eja.eja_algebra import TrivialEJA
2176
2177 EXAMPLES::
2178
2179 sage: J = TrivialEJA()
2180 sage: J.dimension()
2181 0
2182 sage: J.zero()
2183 0
2184 sage: J.one()
2185 0
2186 sage: 7*J.one()*12*J.one()
2187 0
2188 sage: J.one().inner_product(J.one())
2189 0
2190 sage: J.one().norm()
2191 0
2192 sage: J.one().subalgebra_generated_by()
2193 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2194 sage: J.rank()
2195 0
2196
2197 """
2198 def __init__(self, field=AA, **kwargs):
2199 mult_table = []
2200 fdeja = super(TrivialEJA, self)
2201 # The rank is zero using my definition, namely the dimension of the
2202 # largest subalgebra generated by any element.
2203 return fdeja.__init__(field, mult_table, rank=0, **kwargs)