]> gitweb.michael.orlitzky.com - sage.d.git/blob - eja_algebra.py
ec3dd1141a9f497acca26f6d09c830420dfb5d5d
[sage.d.git] / 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 SETUP::
9
10 sage: from mjo.eja.eja_algebra import random_eja
11
12 EXAMPLES::
13
14 sage: random_eja()
15 Euclidean Jordan algebra of dimension...
16
17 """
18
19 from itertools import repeat
20
21 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
22 from sage.categories.magmatic_algebras import MagmaticAlgebras
23 from sage.combinat.free_module import CombinatorialFreeModule
24 from sage.matrix.constructor import matrix
25 from sage.matrix.matrix_space import MatrixSpace
26 from sage.misc.cachefunc import cached_method
27 from sage.misc.table import table
28 from sage.modules.free_module import FreeModule, VectorSpace
29 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
30 PolynomialRing,
31 QuadraticField)
32 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
33 from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
34 from mjo.eja.eja_utils import _mat2vec
35
36 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
37 r"""
38 The lowest-level class for representing a Euclidean Jordan algebra.
39 """
40 def _coerce_map_from_base_ring(self):
41 """
42 Disable the map from the base ring into the algebra.
43
44 Performing a nonsense conversion like this automatically
45 is counterpedagogical. The fallback is to try the usual
46 element constructor, which should also fail.
47
48 SETUP::
49
50 sage: from mjo.eja.eja_algebra import random_eja
51
52 TESTS::
53
54 sage: set_random_seed()
55 sage: J = random_eja()
56 sage: J(1)
57 Traceback (most recent call last):
58 ...
59 ValueError: not an element of this algebra
60
61 """
62 return None
63
64 def __init__(self,
65 field,
66 multiplication_table,
67 inner_product_table,
68 prefix='e',
69 category=None,
70 matrix_basis=None,
71 check_field=True,
72 check_axioms=True):
73 """
74 INPUT:
75
76 * field -- the scalar field for this algebra (must be real)
77
78 * multiplication_table -- the multiplication table for this
79 algebra's implicit basis. Only the lower-triangular portion
80 of the table is used, since the multiplication is assumed
81 to be commutative.
82
83 SETUP::
84
85 sage: from mjo.eja.eja_algebra import (
86 ....: FiniteDimensionalEuclideanJordanAlgebra,
87 ....: JordanSpinEJA,
88 ....: random_eja)
89
90 EXAMPLES:
91
92 By definition, Jordan multiplication commutes::
93
94 sage: set_random_seed()
95 sage: J = random_eja()
96 sage: x,y = J.random_elements(2)
97 sage: x*y == y*x
98 True
99
100 An error is raised if the Jordan product is not commutative::
101
102 sage: JP = ((1,2),(0,0))
103 sage: IP = ((1,0),(0,1))
104 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
105 Traceback (most recent call last):
106 ...
107 ValueError: Jordan product is not commutative
108
109 An error is raised if the inner-product is not commutative::
110
111 sage: JP = ((1,0),(0,1))
112 sage: IP = ((1,2),(0,0))
113 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP)
114 Traceback (most recent call last):
115 ...
116 ValueError: inner-product is not commutative
117
118 TESTS:
119
120 The ``field`` we're given must be real with ``check_field=True``::
121
122 sage: JordanSpinEJA(2,QQbar)
123 Traceback (most recent call last):
124 ...
125 ValueError: scalar field is not real
126
127 The multiplication table must be square with ``check_axioms=True``::
128
129 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),))
130 Traceback (most recent call last):
131 ...
132 ValueError: multiplication table is not square
133
134 The multiplication and inner-product tables must be the same
135 size (and in particular, the inner-product table must also be
136 square) with ``check_axioms=True``::
137
138 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(()))
139 Traceback (most recent call last):
140 ...
141 ValueError: multiplication and inner-product tables are
142 different sizes
143 sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),))
144 Traceback (most recent call last):
145 ...
146 ValueError: multiplication and inner-product tables are
147 different sizes
148
149 """
150 if check_field:
151 if not field.is_subring(RR):
152 # Note: this does return true for the real algebraic
153 # field, the rationals, and any quadratic field where
154 # we've specified a real embedding.
155 raise ValueError("scalar field is not real")
156
157
158 # The multiplication and inner-product tables should be square
159 # if the user wants us to verify them. And we verify them as
160 # soon as possible, because we want to exploit their symmetry.
161 n = len(multiplication_table)
162 if check_axioms:
163 if not all( len(l) == n for l in multiplication_table ):
164 raise ValueError("multiplication table is not square")
165
166 # If the multiplication table is square, we can check if
167 # the inner-product table is square by comparing it to the
168 # multiplication table's dimensions.
169 msg = "multiplication and inner-product tables are different sizes"
170 if not len(inner_product_table) == n:
171 raise ValueError(msg)
172
173 if not all( len(l) == n for l in inner_product_table ):
174 raise ValueError(msg)
175
176 # Check commutativity of the Jordan product (symmetry of
177 # the multiplication table) and the commutativity of the
178 # inner-product (symmetry of the inner-product table)
179 # first if we're going to check them at all.. This has to
180 # be done before we define product_on_basis(), because
181 # that method assumes that self._multiplication_table is
182 # symmetric. And it has to be done before we build
183 # self._inner_product_matrix, because the process used to
184 # construct it assumes symmetry as well.
185 if not all( multiplication_table[j][i]
186 == multiplication_table[i][j]
187 for i in range(n)
188 for j in range(i+1) ):
189 raise ValueError("Jordan product is not commutative")
190
191 if not all( inner_product_table[j][i]
192 == inner_product_table[i][j]
193 for i in range(n)
194 for j in range(i+1) ):
195 raise ValueError("inner-product is not commutative")
196
197 self._matrix_basis = matrix_basis
198
199 if category is None:
200 category = MagmaticAlgebras(field).FiniteDimensional()
201 category = category.WithBasis().Unital()
202
203 fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
204 fda.__init__(field,
205 range(n),
206 prefix=prefix,
207 category=category)
208 self.print_options(bracket='')
209
210 # The multiplication table we're given is necessarily in terms
211 # of vectors, because we don't have an algebra yet for
212 # anything to be an element of. However, it's faster in the
213 # long run to have the multiplication table be in terms of
214 # algebra elements. We do this after calling the superclass
215 # constructor so that from_vector() knows what to do.
216 #
217 # Note: we take advantage of symmetry here, and only store
218 # the lower-triangular portion of the table.
219 self._multiplication_table = [ [ self.vector_space().zero()
220 for j in range(i+1) ]
221 for i in range(n) ]
222
223 for i in range(n):
224 for j in range(i+1):
225 elt = self.from_vector(multiplication_table[i][j])
226 self._multiplication_table[i][j] = elt
227
228 self._multiplication_table = tuple(map(tuple, self._multiplication_table))
229
230 # Save our inner product as a matrix, since the efficiency of
231 # matrix multiplication will usually outweigh the fact that we
232 # have to store a redundant upper- or lower-triangular part.
233 # Pre-cache the fact that these are Hermitian (real symmetric,
234 # in fact) in case some e.g. matrix multiplication routine can
235 # take advantage of it.
236 ip_matrix_constructor = lambda i,j: inner_product_table[i][j] if j <= i else inner_product_table[j][i]
237 self._inner_product_matrix = matrix(field, n, ip_matrix_constructor)
238 self._inner_product_matrix._cache = {'hermitian': True}
239 self._inner_product_matrix.set_immutable()
240
241 if check_axioms:
242 if not self._is_jordanian():
243 raise ValueError("Jordan identity does not hold")
244 if not self._inner_product_is_associative():
245 raise ValueError("inner product is not associative")
246
247 def _element_constructor_(self, elt):
248 """
249 Construct an element of this algebra from its vector or matrix
250 representation.
251
252 This gets called only after the parent element _call_ method
253 fails to find a coercion for the argument.
254
255 SETUP::
256
257 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
258 ....: HadamardEJA,
259 ....: RealSymmetricEJA)
260
261 EXAMPLES:
262
263 The identity in `S^n` is converted to the identity in the EJA::
264
265 sage: J = RealSymmetricEJA(3)
266 sage: I = matrix.identity(QQ,3)
267 sage: J(I) == J.one()
268 True
269
270 This skew-symmetric matrix can't be represented in the EJA::
271
272 sage: J = RealSymmetricEJA(3)
273 sage: A = matrix(QQ,3, lambda i,j: i-j)
274 sage: J(A)
275 Traceback (most recent call last):
276 ...
277 ValueError: not an element of this algebra
278
279 TESTS:
280
281 Ensure that we can convert any element of the two non-matrix
282 simple algebras (whose matrix representations are columns)
283 back and forth faithfully::
284
285 sage: set_random_seed()
286 sage: J = HadamardEJA.random_instance()
287 sage: x = J.random_element()
288 sage: J(x.to_vector().column()) == x
289 True
290 sage: J = JordanSpinEJA.random_instance()
291 sage: x = J.random_element()
292 sage: J(x.to_vector().column()) == x
293 True
294 """
295 msg = "not an element of this algebra"
296 if elt == 0:
297 # The superclass implementation of random_element()
298 # needs to be able to coerce "0" into the algebra.
299 return self.zero()
300 elif elt in self.base_ring():
301 # Ensure that no base ring -> algebra coercion is performed
302 # by this method. There's some stupidity in sage that would
303 # otherwise propagate to this method; for example, sage thinks
304 # that the integer 3 belongs to the space of 2-by-2 matrices.
305 raise ValueError(msg)
306
307 if elt not in self.matrix_space():
308 raise ValueError(msg)
309
310 # Thanks for nothing! Matrix spaces aren't vector spaces in
311 # Sage, so we have to figure out its matrix-basis coordinates
312 # ourselves. We use the basis space's ring instead of the
313 # element's ring because the basis space might be an algebraic
314 # closure whereas the base ring of the 3-by-3 identity matrix
315 # could be QQ instead of QQbar.
316 V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
317 W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() )
318
319 try:
320 coords = W.coordinate_vector(_mat2vec(elt))
321 except ArithmeticError: # vector is not in free module
322 raise ValueError(msg)
323
324 return self.from_vector(coords)
325
326 def _repr_(self):
327 """
328 Return a string representation of ``self``.
329
330 SETUP::
331
332 sage: from mjo.eja.eja_algebra import JordanSpinEJA
333
334 TESTS:
335
336 Ensure that it says what we think it says::
337
338 sage: JordanSpinEJA(2, field=AA)
339 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
340 sage: JordanSpinEJA(3, field=RDF)
341 Euclidean Jordan algebra of dimension 3 over Real Double Field
342
343 """
344 fmt = "Euclidean Jordan algebra of dimension {} over {}"
345 return fmt.format(self.dimension(), self.base_ring())
346
347 def product_on_basis(self, i, j):
348 # We only stored the lower-triangular portion of the
349 # multiplication table.
350 if j <= i:
351 return self._multiplication_table[i][j]
352 else:
353 return self._multiplication_table[j][i]
354
355 def _is_commutative(self):
356 r"""
357 Whether or not this algebra's multiplication table is commutative.
358
359 This method should of course always return ``True``, unless
360 this algebra was constructed with ``check_axioms=False`` and
361 passed an invalid multiplication table.
362 """
363 return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
364 for i in range(self.dimension())
365 for j in range(self.dimension()) )
366
367 def _is_jordanian(self):
368 r"""
369 Whether or not this algebra's multiplication table respects the
370 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
371
372 We only check one arrangement of `x` and `y`, so for a
373 ``True`` result to be truly true, you should also check
374 :meth:`_is_commutative`. This method should of course always
375 return ``True``, unless this algebra was constructed with
376 ``check_axioms=False`` and passed an invalid multiplication table.
377 """
378 return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
379 ==
380 (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
381 for i in range(self.dimension())
382 for j in range(self.dimension()) )
383
384 def _inner_product_is_associative(self):
385 r"""
386 Return whether or not this algebra's inner product `B` is
387 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
388
389 This method should of course always return ``True``, unless
390 this algebra was constructed with ``check_axioms=False`` and
391 passed an invalid multiplication table.
392 """
393
394 # Used to check whether or not something is zero in an inexact
395 # ring. This number is sufficient to allow the construction of
396 # QuaternionHermitianEJA(2, RDF) with check_axioms=True.
397 epsilon = 1e-16
398
399 for i in range(self.dimension()):
400 for j in range(self.dimension()):
401 for k in range(self.dimension()):
402 x = self.monomial(i)
403 y = self.monomial(j)
404 z = self.monomial(k)
405 diff = (x*y).inner_product(z) - x.inner_product(y*z)
406
407 if self.base_ring().is_exact():
408 if diff != 0:
409 return False
410 else:
411 if diff.abs() > epsilon:
412 return False
413
414 return True
415
416 @cached_method
417 def characteristic_polynomial_of(self):
418 """
419 Return the algebra's "characteristic polynomial of" function,
420 which is itself a multivariate polynomial that, when evaluated
421 at the coordinates of some algebra element, returns that
422 element's characteristic polynomial.
423
424 The resulting polynomial has `n+1` variables, where `n` is the
425 dimension of this algebra. The first `n` variables correspond to
426 the coordinates of an algebra element: when evaluated at the
427 coordinates of an algebra element with respect to a certain
428 basis, the result is a univariate polynomial (in the one
429 remaining variable ``t``), namely the characteristic polynomial
430 of that element.
431
432 SETUP::
433
434 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
435
436 EXAMPLES:
437
438 The characteristic polynomial in the spin algebra is given in
439 Alizadeh, Example 11.11::
440
441 sage: J = JordanSpinEJA(3)
442 sage: p = J.characteristic_polynomial_of(); p
443 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
444 sage: xvec = J.one().to_vector()
445 sage: p(*xvec)
446 t^2 - 2*t + 1
447
448 By definition, the characteristic polynomial is a monic
449 degree-zero polynomial in a rank-zero algebra. Note that
450 Cayley-Hamilton is indeed satisfied since the polynomial
451 ``1`` evaluates to the identity element of the algebra on
452 any argument::
453
454 sage: J = TrivialEJA()
455 sage: J.characteristic_polynomial_of()
456 1
457
458 """
459 r = self.rank()
460 n = self.dimension()
461
462 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
463 a = self._charpoly_coefficients()
464
465 # We go to a bit of trouble here to reorder the
466 # indeterminates, so that it's easier to evaluate the
467 # characteristic polynomial at x's coordinates and get back
468 # something in terms of t, which is what we want.
469 S = PolynomialRing(self.base_ring(),'t')
470 t = S.gen(0)
471 if r > 0:
472 R = a[0].parent()
473 S = PolynomialRing(S, R.variable_names())
474 t = S(t)
475
476 return (t**r + sum( a[k]*(t**k) for k in range(r) ))
477
478 def coordinate_polynomial_ring(self):
479 r"""
480 The multivariate polynomial ring in which this algebra's
481 :meth:`characteristic_polynomial_of` lives.
482
483 SETUP::
484
485 sage: from mjo.eja.eja_algebra import (HadamardEJA,
486 ....: RealSymmetricEJA)
487
488 EXAMPLES::
489
490 sage: J = HadamardEJA(2)
491 sage: J.coordinate_polynomial_ring()
492 Multivariate Polynomial Ring in X1, X2...
493 sage: J = RealSymmetricEJA(3,QQ,orthonormalize=False)
494 sage: J.coordinate_polynomial_ring()
495 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
496
497 """
498 var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) )
499 return PolynomialRing(self.base_ring(), var_names)
500
501 def inner_product(self, x, y):
502 """
503 The inner product associated with this Euclidean Jordan algebra.
504
505 Defaults to the trace inner product, but can be overridden by
506 subclasses if they are sure that the necessary properties are
507 satisfied.
508
509 SETUP::
510
511 sage: from mjo.eja.eja_algebra import (random_eja,
512 ....: HadamardEJA,
513 ....: BilinearFormEJA)
514
515 EXAMPLES:
516
517 Our inner product is "associative," which means the following for
518 a symmetric bilinear form::
519
520 sage: set_random_seed()
521 sage: J = random_eja()
522 sage: x,y,z = J.random_elements(3)
523 sage: (x*y).inner_product(z) == y.inner_product(x*z)
524 True
525
526 TESTS:
527
528 Ensure that this is the usual inner product for the algebras
529 over `R^n`::
530
531 sage: set_random_seed()
532 sage: J = HadamardEJA.random_instance()
533 sage: x,y = J.random_elements(2)
534 sage: actual = x.inner_product(y)
535 sage: expected = x.to_vector().inner_product(y.to_vector())
536 sage: actual == expected
537 True
538
539 Ensure that this is one-half of the trace inner-product in a
540 BilinearFormEJA that isn't just the reals (when ``n`` isn't
541 one). This is in Faraut and Koranyi, and also my "On the
542 symmetry..." paper::
543
544 sage: set_random_seed()
545 sage: J = BilinearFormEJA.random_instance()
546 sage: n = J.dimension()
547 sage: x = J.random_element()
548 sage: y = J.random_element()
549 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
550 True
551 """
552 B = self._inner_product_matrix
553 return (B*x.to_vector()).inner_product(y.to_vector())
554
555
556 def is_trivial(self):
557 """
558 Return whether or not this algebra is trivial.
559
560 A trivial algebra contains only the zero element.
561
562 SETUP::
563
564 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
565 ....: TrivialEJA)
566
567 EXAMPLES::
568
569 sage: J = ComplexHermitianEJA(3)
570 sage: J.is_trivial()
571 False
572
573 ::
574
575 sage: J = TrivialEJA()
576 sage: J.is_trivial()
577 True
578
579 """
580 return self.dimension() == 0
581
582
583 def multiplication_table(self):
584 """
585 Return a visual representation of this algebra's multiplication
586 table (on basis elements).
587
588 SETUP::
589
590 sage: from mjo.eja.eja_algebra import JordanSpinEJA
591
592 EXAMPLES::
593
594 sage: J = JordanSpinEJA(4)
595 sage: J.multiplication_table()
596 +----++----+----+----+----+
597 | * || e0 | e1 | e2 | e3 |
598 +====++====+====+====+====+
599 | e0 || e0 | e1 | e2 | e3 |
600 +----++----+----+----+----+
601 | e1 || e1 | e0 | 0 | 0 |
602 +----++----+----+----+----+
603 | e2 || e2 | 0 | e0 | 0 |
604 +----++----+----+----+----+
605 | e3 || e3 | 0 | 0 | e0 |
606 +----++----+----+----+----+
607
608 """
609 n = self.dimension()
610 M = [ [ self.zero() for j in range(n) ]
611 for i in range(n) ]
612 for i in range(n):
613 for j in range(i+1):
614 M[i][j] = self._multiplication_table[i][j]
615 M[j][i] = M[i][j]
616
617 for i in range(n):
618 # Prepend the left "header" column entry Can't do this in
619 # the loop because it messes up the symmetry.
620 M[i] = [self.monomial(i)] + M[i]
621
622 # Prepend the header row.
623 M = [["*"] + list(self.gens())] + M
624 return table(M, header_row=True, header_column=True, frame=True)
625
626
627 def matrix_basis(self):
628 """
629 Return an (often more natural) representation of this algebras
630 basis as an ordered tuple of matrices.
631
632 Every finite-dimensional Euclidean Jordan Algebra is a, up to
633 Jordan isomorphism, a direct sum of five simple
634 algebras---four of which comprise Hermitian matrices. And the
635 last type of algebra can of course be thought of as `n`-by-`1`
636 column matrices (ambiguusly called column vectors) to avoid
637 special cases. As a result, matrices (and column vectors) are
638 a natural representation format for Euclidean Jordan algebra
639 elements.
640
641 But, when we construct an algebra from a basis of matrices,
642 those matrix representations are lost in favor of coordinate
643 vectors *with respect to* that basis. We could eventually
644 convert back if we tried hard enough, but having the original
645 representations handy is valuable enough that we simply store
646 them and return them from this method.
647
648 Why implement this for non-matrix algebras? Avoiding special
649 cases for the :class:`BilinearFormEJA` pays with simplicity in
650 its own right. But mainly, we would like to be able to assume
651 that elements of a :class:`DirectSumEJA` can be displayed
652 nicely, without having to have special classes for direct sums
653 one of whose components was a matrix algebra.
654
655 SETUP::
656
657 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
658 ....: RealSymmetricEJA)
659
660 EXAMPLES::
661
662 sage: J = RealSymmetricEJA(2)
663 sage: J.basis()
664 Finite family {0: e0, 1: e1, 2: e2}
665 sage: J.matrix_basis()
666 (
667 [1 0] [ 0 0.7071067811865475?] [0 0]
668 [0 0], [0.7071067811865475? 0], [0 1]
669 )
670
671 ::
672
673 sage: J = JordanSpinEJA(2)
674 sage: J.basis()
675 Finite family {0: e0, 1: e1}
676 sage: J.matrix_basis()
677 (
678 [1] [0]
679 [0], [1]
680 )
681 """
682 if self._matrix_basis is None:
683 M = self.matrix_space()
684 return tuple( M(b.to_vector()) for b in self.basis() )
685 else:
686 return self._matrix_basis
687
688
689 def matrix_space(self):
690 """
691 Return the matrix space in which this algebra's elements live, if
692 we think of them as matrices (including column vectors of the
693 appropriate size).
694
695 Generally this will be an `n`-by-`1` column-vector space,
696 except when the algebra is trivial. There it's `n`-by-`n`
697 (where `n` is zero), to ensure that two elements of the matrix
698 space (empty matrices) can be multiplied.
699
700 Matrix algebras override this with something more useful.
701 """
702 if self.is_trivial():
703 return MatrixSpace(self.base_ring(), 0)
704 elif self._matrix_basis is None or len(self._matrix_basis) == 0:
705 return MatrixSpace(self.base_ring(), self.dimension(), 1)
706 else:
707 return self._matrix_basis[0].matrix_space()
708
709
710 @cached_method
711 def one(self):
712 """
713 Return the unit element of this algebra.
714
715 SETUP::
716
717 sage: from mjo.eja.eja_algebra import (HadamardEJA,
718 ....: random_eja)
719
720 EXAMPLES::
721
722 sage: J = HadamardEJA(5)
723 sage: J.one()
724 e0 + e1 + e2 + e3 + e4
725
726 TESTS:
727
728 The identity element acts like the identity::
729
730 sage: set_random_seed()
731 sage: J = random_eja()
732 sage: x = J.random_element()
733 sage: J.one()*x == x and x*J.one() == x
734 True
735
736 The matrix of the unit element's operator is the identity::
737
738 sage: set_random_seed()
739 sage: J = random_eja()
740 sage: actual = J.one().operator().matrix()
741 sage: expected = matrix.identity(J.base_ring(), J.dimension())
742 sage: actual == expected
743 True
744
745 Ensure that the cached unit element (often precomputed by
746 hand) agrees with the computed one::
747
748 sage: set_random_seed()
749 sage: J = random_eja()
750 sage: cached = J.one()
751 sage: J.one.clear_cache()
752 sage: J.one() == cached
753 True
754
755 """
756 # We can brute-force compute the matrices of the operators
757 # that correspond to the basis elements of this algebra.
758 # If some linear combination of those basis elements is the
759 # algebra identity, then the same linear combination of
760 # their matrices has to be the identity matrix.
761 #
762 # Of course, matrices aren't vectors in sage, so we have to
763 # appeal to the "long vectors" isometry.
764 oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
765
766 # Now we use basic linear algebra to find the coefficients,
767 # of the matrices-as-vectors-linear-combination, which should
768 # work for the original algebra basis too.
769 A = matrix(self.base_ring(), oper_vecs)
770
771 # We used the isometry on the left-hand side already, but we
772 # still need to do it for the right-hand side. Recall that we
773 # wanted something that summed to the identity matrix.
774 b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
775
776 # Now if there's an identity element in the algebra, this
777 # should work. We solve on the left to avoid having to
778 # transpose the matrix "A".
779 return self.from_vector(A.solve_left(b))
780
781
782 def peirce_decomposition(self, c):
783 """
784 The Peirce decomposition of this algebra relative to the
785 idempotent ``c``.
786
787 In the future, this can be extended to a complete system of
788 orthogonal idempotents.
789
790 INPUT:
791
792 - ``c`` -- an idempotent of this algebra.
793
794 OUTPUT:
795
796 A triple (J0, J5, J1) containing two subalgebras and one subspace
797 of this algebra,
798
799 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
800 corresponding to the eigenvalue zero.
801
802 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
803 corresponding to the eigenvalue one-half.
804
805 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
806 corresponding to the eigenvalue one.
807
808 These are the only possible eigenspaces for that operator, and this
809 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
810 orthogonal, and are subalgebras of this algebra with the appropriate
811 restrictions.
812
813 SETUP::
814
815 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
816
817 EXAMPLES:
818
819 The canonical example comes from the symmetric matrices, which
820 decompose into diagonal and off-diagonal parts::
821
822 sage: J = RealSymmetricEJA(3)
823 sage: C = matrix(QQ, [ [1,0,0],
824 ....: [0,1,0],
825 ....: [0,0,0] ])
826 sage: c = J(C)
827 sage: J0,J5,J1 = J.peirce_decomposition(c)
828 sage: J0
829 Euclidean Jordan algebra of dimension 1...
830 sage: J5
831 Vector space of degree 6 and dimension 2...
832 sage: J1
833 Euclidean Jordan algebra of dimension 3...
834 sage: J0.one().to_matrix()
835 [0 0 0]
836 [0 0 0]
837 [0 0 1]
838 sage: orig_df = AA.options.display_format
839 sage: AA.options.display_format = 'radical'
840 sage: J.from_vector(J5.basis()[0]).to_matrix()
841 [ 0 0 1/2*sqrt(2)]
842 [ 0 0 0]
843 [1/2*sqrt(2) 0 0]
844 sage: J.from_vector(J5.basis()[1]).to_matrix()
845 [ 0 0 0]
846 [ 0 0 1/2*sqrt(2)]
847 [ 0 1/2*sqrt(2) 0]
848 sage: AA.options.display_format = orig_df
849 sage: J1.one().to_matrix()
850 [1 0 0]
851 [0 1 0]
852 [0 0 0]
853
854 TESTS:
855
856 Every algebra decomposes trivially with respect to its identity
857 element::
858
859 sage: set_random_seed()
860 sage: J = random_eja()
861 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
862 sage: J0.dimension() == 0 and J5.dimension() == 0
863 True
864 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
865 True
866
867 The decomposition is into eigenspaces, and its components are
868 therefore necessarily orthogonal. Moreover, the identity
869 elements in the two subalgebras are the projections onto their
870 respective subspaces of the superalgebra's identity element::
871
872 sage: set_random_seed()
873 sage: J = random_eja()
874 sage: x = J.random_element()
875 sage: if not J.is_trivial():
876 ....: while x.is_nilpotent():
877 ....: x = J.random_element()
878 sage: c = x.subalgebra_idempotent()
879 sage: J0,J5,J1 = J.peirce_decomposition(c)
880 sage: ipsum = 0
881 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
882 ....: w = w.superalgebra_element()
883 ....: y = J.from_vector(y)
884 ....: z = z.superalgebra_element()
885 ....: ipsum += w.inner_product(y).abs()
886 ....: ipsum += w.inner_product(z).abs()
887 ....: ipsum += y.inner_product(z).abs()
888 sage: ipsum
889 0
890 sage: J1(c) == J1.one()
891 True
892 sage: J0(J.one() - c) == J0.one()
893 True
894
895 """
896 if not c.is_idempotent():
897 raise ValueError("element is not idempotent: %s" % c)
898
899 from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
900
901 # Default these to what they should be if they turn out to be
902 # trivial, because eigenspaces_left() won't return eigenvalues
903 # corresponding to trivial spaces (e.g. it returns only the
904 # eigenspace corresponding to lambda=1 if you take the
905 # decomposition relative to the identity element).
906 trivial = FiniteDimensionalEuclideanJordanSubalgebra(self, ())
907 J0 = trivial # eigenvalue zero
908 J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
909 J1 = trivial # eigenvalue one
910
911 for (eigval, eigspace) in c.operator().matrix().right_eigenspaces():
912 if eigval == ~(self.base_ring()(2)):
913 J5 = eigspace
914 else:
915 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
916 subalg = FiniteDimensionalEuclideanJordanSubalgebra(self,
917 gens,
918 check_axioms=False)
919 if eigval == 0:
920 J0 = subalg
921 elif eigval == 1:
922 J1 = subalg
923 else:
924 raise ValueError("unexpected eigenvalue: %s" % eigval)
925
926 return (J0, J5, J1)
927
928
929 def random_element(self, thorough=False):
930 r"""
931 Return a random element of this algebra.
932
933 Our algebra superclass method only returns a linear
934 combination of at most two basis elements. We instead
935 want the vector space "random element" method that
936 returns a more diverse selection.
937
938 INPUT:
939
940 - ``thorough`` -- (boolean; default False) whether or not we
941 should generate irrational coefficients for the random
942 element when our base ring is irrational; this slows the
943 algebra operations to a crawl, but any truly random method
944 should include them
945
946 """
947 # For a general base ring... maybe we can trust this to do the
948 # right thing? Unlikely, but.
949 V = self.vector_space()
950 v = V.random_element()
951
952 if self.base_ring() is AA:
953 # The "random element" method of the algebraic reals is
954 # stupid at the moment, and only returns integers between
955 # -2 and 2, inclusive:
956 #
957 # https://trac.sagemath.org/ticket/30875
958 #
959 # Instead, we implement our own "random vector" method,
960 # and then coerce that into the algebra. We use the vector
961 # space degree here instead of the dimension because a
962 # subalgebra could (for example) be spanned by only two
963 # vectors, each with five coordinates. We need to
964 # generate all five coordinates.
965 if thorough:
966 v *= QQbar.random_element().real()
967 else:
968 v *= QQ.random_element()
969
970 return self.from_vector(V.coordinate_vector(v))
971
972 def random_elements(self, count, thorough=False):
973 """
974 Return ``count`` random elements as a tuple.
975
976 INPUT:
977
978 - ``thorough`` -- (boolean; default False) whether or not we
979 should generate irrational coefficients for the random
980 elements when our base ring is irrational; this slows the
981 algebra operations to a crawl, but any truly random method
982 should include them
983
984 SETUP::
985
986 sage: from mjo.eja.eja_algebra import JordanSpinEJA
987
988 EXAMPLES::
989
990 sage: J = JordanSpinEJA(3)
991 sage: x,y,z = J.random_elements(3)
992 sage: all( [ x in J, y in J, z in J ])
993 True
994 sage: len( J.random_elements(10) ) == 10
995 True
996
997 """
998 return tuple( self.random_element(thorough)
999 for idx in range(count) )
1000
1001
1002 @cached_method
1003 def _charpoly_coefficients(self):
1004 r"""
1005 The `r` polynomial coefficients of the "characteristic polynomial
1006 of" function.
1007 """
1008 n = self.dimension()
1009 R = self.coordinate_polynomial_ring()
1010 vars = R.gens()
1011 F = R.fraction_field()
1012
1013 def L_x_i_j(i,j):
1014 # From a result in my book, these are the entries of the
1015 # basis representation of L_x.
1016 return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
1017 for k in range(n) )
1018
1019 L_x = matrix(F, n, n, L_x_i_j)
1020
1021 r = None
1022 if self.rank.is_in_cache():
1023 r = self.rank()
1024 # There's no need to pad the system with redundant
1025 # columns if we *know* they'll be redundant.
1026 n = r
1027
1028 # Compute an extra power in case the rank is equal to
1029 # the dimension (otherwise, we would stop at x^(r-1)).
1030 x_powers = [ (L_x**k)*self.one().to_vector()
1031 for k in range(n+1) ]
1032 A = matrix.column(F, x_powers[:n])
1033 AE = A.extended_echelon_form()
1034 E = AE[:,n:]
1035 A_rref = AE[:,:n]
1036 if r is None:
1037 r = A_rref.rank()
1038 b = x_powers[r]
1039
1040 # The theory says that only the first "r" coefficients are
1041 # nonzero, and they actually live in the original polynomial
1042 # ring and not the fraction field. We negate them because
1043 # in the actual characteristic polynomial, they get moved
1044 # to the other side where x^r lives.
1045 return -A_rref.solve_right(E*b).change_ring(R)[:r]
1046
1047 @cached_method
1048 def rank(self):
1049 r"""
1050 Return the rank of this EJA.
1051
1052 This is a cached method because we know the rank a priori for
1053 all of the algebras we can construct. Thus we can avoid the
1054 expensive ``_charpoly_coefficients()`` call unless we truly
1055 need to compute the whole characteristic polynomial.
1056
1057 SETUP::
1058
1059 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1060 ....: JordanSpinEJA,
1061 ....: RealSymmetricEJA,
1062 ....: ComplexHermitianEJA,
1063 ....: QuaternionHermitianEJA,
1064 ....: random_eja)
1065
1066 EXAMPLES:
1067
1068 The rank of the Jordan spin algebra is always two::
1069
1070 sage: JordanSpinEJA(2).rank()
1071 2
1072 sage: JordanSpinEJA(3).rank()
1073 2
1074 sage: JordanSpinEJA(4).rank()
1075 2
1076
1077 The rank of the `n`-by-`n` Hermitian real, complex, or
1078 quaternion matrices is `n`::
1079
1080 sage: RealSymmetricEJA(4).rank()
1081 4
1082 sage: ComplexHermitianEJA(3).rank()
1083 3
1084 sage: QuaternionHermitianEJA(2).rank()
1085 2
1086
1087 TESTS:
1088
1089 Ensure that every EJA that we know how to construct has a
1090 positive integer rank, unless the algebra is trivial in
1091 which case its rank will be zero::
1092
1093 sage: set_random_seed()
1094 sage: J = random_eja()
1095 sage: r = J.rank()
1096 sage: r in ZZ
1097 True
1098 sage: r > 0 or (r == 0 and J.is_trivial())
1099 True
1100
1101 Ensure that computing the rank actually works, since the ranks
1102 of all simple algebras are known and will be cached by default::
1103
1104 sage: set_random_seed() # long time
1105 sage: J = random_eja() # long time
1106 sage: caches = J.rank() # long time
1107 sage: J.rank.clear_cache() # long time
1108 sage: J.rank() == cached # long time
1109 True
1110
1111 """
1112 return len(self._charpoly_coefficients())
1113
1114
1115 def vector_space(self):
1116 """
1117 Return the vector space that underlies this algebra.
1118
1119 SETUP::
1120
1121 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1122
1123 EXAMPLES::
1124
1125 sage: J = RealSymmetricEJA(2)
1126 sage: J.vector_space()
1127 Vector space of dimension 3 over...
1128
1129 """
1130 return self.zero().to_vector().parent().ambient_vector_space()
1131
1132
1133 Element = FiniteDimensionalEuclideanJordanAlgebraElement
1134
1135 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
1136 r"""
1137 New class for algebras whose supplied basis elements have all rational entries.
1138
1139 SETUP::
1140
1141 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1142
1143 EXAMPLES:
1144
1145 The supplied basis is orthonormalized by default::
1146
1147 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1148 sage: J = BilinearFormEJA(B)
1149 sage: J.matrix_basis()
1150 (
1151 [1] [ 0] [ 0]
1152 [0] [1/5] [32/5]
1153 [0], [ 0], [ 5]
1154 )
1155
1156 """
1157 def __init__(self,
1158 field,
1159 basis,
1160 jordan_product,
1161 inner_product,
1162 orthonormalize=True,
1163 prefix='e',
1164 category=None,
1165 check_field=True,
1166 check_axioms=True):
1167
1168 n = len(basis)
1169 vector_basis = basis
1170
1171 from sage.structure.element import is_Matrix
1172 basis_is_matrices = False
1173
1174 degree = 0
1175 if n > 0:
1176 if is_Matrix(basis[0]):
1177 basis_is_matrices = True
1178 from mjo.eja.eja_utils import _vec2mat
1179 vector_basis = tuple( map(_mat2vec,basis) )
1180 degree = basis[0].nrows()**2
1181 else:
1182 degree = basis[0].degree()
1183
1184 V = VectorSpace(field, degree)
1185
1186 # If we were asked to orthonormalize, and if the orthonormal
1187 # basis is different from the given one, then we also want to
1188 # compute multiplication and inner-product tables for the
1189 # deorthonormalized basis. These can be used later to
1190 # construct a deorthonormalized copy of this algebra over QQ
1191 # in which several operations are much faster.
1192 self._rational_algebra = None
1193
1194 if orthonormalize:
1195 if self.base_ring() is not QQ:
1196 # There's no point in constructing the extra algebra if this
1197 # one is already rational. If the original basis is rational
1198 # but normalization would make it irrational, then this whole
1199 # constructor will just fail anyway as it tries to stick an
1200 # irrational number into a rational algebra.
1201 #
1202 # Note: the same Jordan and inner-products work here,
1203 # because they are necessarily defined with respect to
1204 # ambient coordinates and not any particular basis.
1205 self._rational_algebra = RationalBasisEuclideanJordanAlgebra(
1206 QQ,
1207 basis,
1208 jordan_product,
1209 inner_product,
1210 orthonormalize=False,
1211 prefix=prefix,
1212 category=category,
1213 check_field=False,
1214 check_axioms=False)
1215
1216 # Compute the deorthonormalized tables before we orthonormalize
1217 # the given basis.
1218 W = V.span_of_basis( vector_basis )
1219
1220 # Note: the Jordan and inner-products are defined in terms
1221 # of the ambient basis. It's important that their arguments
1222 # are in ambient coordinates as well.
1223 for i in range(n):
1224 for j in range(i+1):
1225 # given basis w.r.t. ambient coords
1226 q_i = vector_basis[i]
1227 q_j = vector_basis[j]
1228
1229 if basis_is_matrices:
1230 q_i = _vec2mat(q_i)
1231 q_j = _vec2mat(q_j)
1232
1233 elt = jordan_product(q_i, q_j)
1234 ip = inner_product(q_i, q_j)
1235
1236 if basis_is_matrices:
1237 # do another mat2vec because the multiplication
1238 # table is in terms of vectors
1239 elt = _mat2vec(elt)
1240
1241 # We overwrite the name "vector_basis" in a second, but never modify it
1242 # in place, to this effectively makes a copy of it.
1243 deortho_vector_basis = vector_basis
1244 self._deortho_matrix = None
1245
1246 if orthonormalize:
1247 from mjo.eja.eja_utils import gram_schmidt
1248 if basis_is_matrices:
1249 vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y))
1250 vector_basis = gram_schmidt(vector_basis, vector_ip)
1251 else:
1252 vector_basis = gram_schmidt(vector_basis, inner_product)
1253
1254 W = V.span_of_basis( vector_basis )
1255
1256 # Normalize the "matrix" basis, too!
1257 basis = vector_basis
1258
1259 if basis_is_matrices:
1260 basis = tuple( map(_vec2mat,basis) )
1261
1262 W = V.span_of_basis( vector_basis )
1263
1264 # Now "W" is the vector space of our algebra coordinates. The
1265 # variables "X1", "X2",... refer to the entries of vectors in
1266 # W. Thus to convert back and forth between the orthonormal
1267 # coordinates and the given ones, we need to stick the original
1268 # basis in W.
1269 U = V.span_of_basis( deortho_vector_basis )
1270 self._deortho_matrix = matrix( U.coordinate_vector(q)
1271 for q in vector_basis )
1272
1273 # If the superclass constructor is going to verify the
1274 # symmetry of this table, it has better at least be
1275 # square...
1276 if check_axioms:
1277 mult_table = [ [0 for j in range(n)] for i in range(n) ]
1278 ip_table = [ [0 for j in range(n)] for i in range(n) ]
1279 else:
1280 mult_table = [ [0 for j in range(i+1)] for i in range(n) ]
1281 ip_table = [ [0 for j in range(i+1)] for i in range(n) ]
1282
1283 # Note: the Jordan and inner-products are defined in terms
1284 # of the ambient basis. It's important that their arguments
1285 # are in ambient coordinates as well.
1286 for i in range(n):
1287 for j in range(i+1):
1288 # ortho basis w.r.t. ambient coords
1289 q_i = vector_basis[i]
1290 q_j = vector_basis[j]
1291
1292 if basis_is_matrices:
1293 q_i = _vec2mat(q_i)
1294 q_j = _vec2mat(q_j)
1295
1296 elt = jordan_product(q_i, q_j)
1297 ip = inner_product(q_i, q_j)
1298
1299 if basis_is_matrices:
1300 # do another mat2vec because the multiplication
1301 # table is in terms of vectors
1302 elt = _mat2vec(elt)
1303
1304 elt = W.coordinate_vector(elt)
1305 mult_table[i][j] = elt
1306 ip_table[i][j] = ip
1307 if check_axioms:
1308 # The tables are square if we're verifying that they
1309 # are commutative.
1310 mult_table[j][i] = elt
1311 ip_table[j][i] = ip
1312
1313 if basis_is_matrices:
1314 for m in basis:
1315 m.set_immutable()
1316 else:
1317 basis = tuple( x.column() for x in basis )
1318
1319 super().__init__(field,
1320 mult_table,
1321 ip_table,
1322 prefix,
1323 category,
1324 basis, # matrix basis
1325 check_field,
1326 check_axioms)
1327
1328 @cached_method
1329 def _charpoly_coefficients(self):
1330 r"""
1331 SETUP::
1332
1333 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1334 ....: JordanSpinEJA)
1335
1336 EXAMPLES:
1337
1338 The base ring of the resulting polynomial coefficients is what
1339 it should be, and not the rationals (unless the algebra was
1340 already over the rationals)::
1341
1342 sage: J = JordanSpinEJA(3)
1343 sage: J._charpoly_coefficients()
1344 (X1^2 - X2^2 - X3^2, -2*X1)
1345 sage: a0 = J._charpoly_coefficients()[0]
1346 sage: J.base_ring()
1347 Algebraic Real Field
1348 sage: a0.base_ring()
1349 Algebraic Real Field
1350
1351 """
1352 if self.base_ring() is QQ:
1353 # There's no need to construct *another* algebra over the
1354 # rationals if this one is already over the rationals.
1355 superclass = super(RationalBasisEuclideanJordanAlgebra, self)
1356 return superclass._charpoly_coefficients()
1357
1358 # Do the computation over the rationals. The answer will be
1359 # the same, because all we've done is a change of basis.
1360 # Then, change back from QQ to our real base ring
1361 a = ( a_i.change_ring(self.base_ring())
1362 for a_i in self._rational_algebra._charpoly_coefficients() )
1363
1364 # Now convert the coordinate variables back to the
1365 # deorthonormalized ones.
1366 R = self.coordinate_polynomial_ring()
1367 from sage.modules.free_module_element import vector
1368 X = vector(R, R.gens())
1369 BX = self._deortho_matrix*X
1370
1371 subs_dict = { X[i]: BX[i] for i in range(len(X)) }
1372 return tuple( a_i.subs(subs_dict) for a_i in a )
1373
1374 class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
1375 r"""
1376 A class for the Euclidean Jordan algebras that we know by name.
1377
1378 These are the Jordan algebras whose basis, multiplication table,
1379 rank, and so on are known a priori. More to the point, they are
1380 the Euclidean Jordan algebras for which we are able to conjure up
1381 a "random instance."
1382
1383 SETUP::
1384
1385 sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
1386
1387 TESTS:
1388
1389 Our basis is normalized with respect to the algebra's inner
1390 product, unless we specify otherwise::
1391
1392 sage: set_random_seed()
1393 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1394 sage: all( b.norm() == 1 for b in J.gens() )
1395 True
1396
1397 Since our basis is orthonormal with respect to the algebra's inner
1398 product, and since we know that this algebra is an EJA, any
1399 left-multiplication operator's matrix will be symmetric because
1400 natural->EJA basis representation is an isometry and within the
1401 EJA the operator is self-adjoint by the Jordan axiom::
1402
1403 sage: set_random_seed()
1404 sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
1405 sage: x = J.random_element()
1406 sage: x.operator().is_self_adjoint()
1407 True
1408 """
1409
1410 @staticmethod
1411 def _max_random_instance_size():
1412 """
1413 Return an integer "size" that is an upper bound on the size of
1414 this algebra when it is used in a random test
1415 case. Unfortunately, the term "size" is ambiguous -- when
1416 dealing with `R^n` under either the Hadamard or Jordan spin
1417 product, the "size" refers to the dimension `n`. When dealing
1418 with a matrix algebra (real symmetric or complex/quaternion
1419 Hermitian), it refers to the size of the matrix, which is far
1420 less than the dimension of the underlying vector space.
1421
1422 This method must be implemented in each subclass.
1423 """
1424 raise NotImplementedError
1425
1426 @classmethod
1427 def random_instance(cls, *args, **kwargs):
1428 """
1429 Return a random instance of this type of algebra.
1430
1431 This method should be implemented in each subclass.
1432 """
1433 from sage.misc.prandom import choice
1434 eja_class = choice(cls.__subclasses__())
1435
1436 # These all bubble up to the RationalBasisEuclideanJordanAlgebra
1437 # superclass constructor, so any (kw)args valid there are also
1438 # valid here.
1439 return eja_class.random_instance(*args, **kwargs)
1440
1441
1442 class MatrixEuclideanJordanAlgebra:
1443 @staticmethod
1444 def real_embed(M):
1445 """
1446 Embed the matrix ``M`` into a space of real matrices.
1447
1448 The matrix ``M`` can have entries in any field at the moment:
1449 the real numbers, complex numbers, or quaternions. And although
1450 they are not a field, we can probably support octonions at some
1451 point, too. This function returns a real matrix that "acts like"
1452 the original with respect to matrix multiplication; i.e.
1453
1454 real_embed(M*N) = real_embed(M)*real_embed(N)
1455
1456 """
1457 raise NotImplementedError
1458
1459
1460 @staticmethod
1461 def real_unembed(M):
1462 """
1463 The inverse of :meth:`real_embed`.
1464 """
1465 raise NotImplementedError
1466
1467 @staticmethod
1468 def jordan_product(X,Y):
1469 return (X*Y + Y*X)/2
1470
1471 @classmethod
1472 def trace_inner_product(cls,X,Y):
1473 Xu = cls.real_unembed(X)
1474 Yu = cls.real_unembed(Y)
1475 tr = (Xu*Yu).trace()
1476
1477 try:
1478 # Works in QQ, AA, RDF, et cetera.
1479 return tr.real()
1480 except AttributeError:
1481 # A quaternion doesn't have a real() method, but does
1482 # have coefficient_tuple() method that returns the
1483 # coefficients of 1, i, j, and k -- in that order.
1484 return tr.coefficient_tuple()[0]
1485
1486
1487 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1488 @staticmethod
1489 def real_embed(M):
1490 """
1491 The identity function, for embedding real matrices into real
1492 matrices.
1493 """
1494 return M
1495
1496 @staticmethod
1497 def real_unembed(M):
1498 """
1499 The identity function, for unembedding real matrices from real
1500 matrices.
1501 """
1502 return M
1503
1504
1505 class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
1506 RealMatrixEuclideanJordanAlgebra):
1507 """
1508 The rank-n simple EJA consisting of real symmetric n-by-n
1509 matrices, the usual symmetric Jordan product, and the trace inner
1510 product. It has dimension `(n^2 + n)/2` over the reals.
1511
1512 SETUP::
1513
1514 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1515
1516 EXAMPLES::
1517
1518 sage: J = RealSymmetricEJA(2)
1519 sage: e0, e1, e2 = J.gens()
1520 sage: e0*e0
1521 e0
1522 sage: e1*e1
1523 1/2*e0 + 1/2*e2
1524 sage: e2*e2
1525 e2
1526
1527 In theory, our "field" can be any subfield of the reals::
1528
1529 sage: RealSymmetricEJA(2, RDF)
1530 Euclidean Jordan algebra of dimension 3 over Real Double Field
1531 sage: RealSymmetricEJA(2, RR)
1532 Euclidean Jordan algebra of dimension 3 over Real Field with
1533 53 bits of precision
1534
1535 TESTS:
1536
1537 The dimension of this algebra is `(n^2 + n) / 2`::
1538
1539 sage: set_random_seed()
1540 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1541 sage: n = ZZ.random_element(1, n_max)
1542 sage: J = RealSymmetricEJA(n)
1543 sage: J.dimension() == (n^2 + n)/2
1544 True
1545
1546 The Jordan multiplication is what we think it is::
1547
1548 sage: set_random_seed()
1549 sage: J = RealSymmetricEJA.random_instance()
1550 sage: x,y = J.random_elements(2)
1551 sage: actual = (x*y).to_matrix()
1552 sage: X = x.to_matrix()
1553 sage: Y = y.to_matrix()
1554 sage: expected = (X*Y + Y*X)/2
1555 sage: actual == expected
1556 True
1557 sage: J(expected) == x*y
1558 True
1559
1560 We can change the generator prefix::
1561
1562 sage: RealSymmetricEJA(3, prefix='q').gens()
1563 (q0, q1, q2, q3, q4, q5)
1564
1565 We can construct the (trivial) algebra of rank zero::
1566
1567 sage: RealSymmetricEJA(0)
1568 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1569
1570 """
1571 @classmethod
1572 def _denormalized_basis(cls, n, field):
1573 """
1574 Return a basis for the space of real symmetric n-by-n matrices.
1575
1576 SETUP::
1577
1578 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1579
1580 TESTS::
1581
1582 sage: set_random_seed()
1583 sage: n = ZZ.random_element(1,5)
1584 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1585 sage: all( M.is_symmetric() for M in B)
1586 True
1587
1588 """
1589 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1590 # coordinates.
1591 S = []
1592 for i in range(n):
1593 for j in range(i+1):
1594 Eij = matrix(field, n, lambda k,l: k==i and l==j)
1595 if i == j:
1596 Sij = Eij
1597 else:
1598 Sij = Eij + Eij.transpose()
1599 S.append(Sij)
1600 return tuple(S)
1601
1602
1603 @staticmethod
1604 def _max_random_instance_size():
1605 return 4 # Dimension 10
1606
1607 @classmethod
1608 def random_instance(cls, field=AA, **kwargs):
1609 """
1610 Return a random instance of this type of algebra.
1611 """
1612 n = ZZ.random_element(cls._max_random_instance_size() + 1)
1613 return cls(n, field, **kwargs)
1614
1615 def __init__(self, n, field=AA, **kwargs):
1616 basis = self._denormalized_basis(n, field)
1617 super(RealSymmetricEJA, self).__init__(field,
1618 basis,
1619 self.jordan_product,
1620 self.trace_inner_product,
1621 **kwargs)
1622 self.rank.set_cache(n)
1623 self.one.set_cache(self(matrix.identity(field,n)))
1624
1625
1626 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1627 @staticmethod
1628 def real_embed(M):
1629 """
1630 Embed the n-by-n complex matrix ``M`` into the space of real
1631 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1632 bi` to the block matrix ``[[a,b],[-b,a]]``.
1633
1634 SETUP::
1635
1636 sage: from mjo.eja.eja_algebra import \
1637 ....: ComplexMatrixEuclideanJordanAlgebra
1638
1639 EXAMPLES::
1640
1641 sage: F = QuadraticField(-1, 'I')
1642 sage: x1 = F(4 - 2*i)
1643 sage: x2 = F(1 + 2*i)
1644 sage: x3 = F(-i)
1645 sage: x4 = F(6)
1646 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1647 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1648 [ 4 -2| 1 2]
1649 [ 2 4|-2 1]
1650 [-----+-----]
1651 [ 0 -1| 6 0]
1652 [ 1 0| 0 6]
1653
1654 TESTS:
1655
1656 Embedding is a homomorphism (isomorphism, in fact)::
1657
1658 sage: set_random_seed()
1659 sage: n = ZZ.random_element(3)
1660 sage: F = QuadraticField(-1, 'I')
1661 sage: X = random_matrix(F, n)
1662 sage: Y = random_matrix(F, n)
1663 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1664 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1665 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1666 sage: Xe*Ye == XYe
1667 True
1668
1669 """
1670 n = M.nrows()
1671 if M.ncols() != n:
1672 raise ValueError("the matrix 'M' must be square")
1673
1674 # We don't need any adjoined elements...
1675 field = M.base_ring().base_ring()
1676
1677 blocks = []
1678 for z in M.list():
1679 a = z.list()[0] # real part, I guess
1680 b = z.list()[1] # imag part, I guess
1681 blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1682
1683 return matrix.block(field, n, blocks)
1684
1685
1686 @staticmethod
1687 def real_unembed(M):
1688 """
1689 The inverse of _embed_complex_matrix().
1690
1691 SETUP::
1692
1693 sage: from mjo.eja.eja_algebra import \
1694 ....: ComplexMatrixEuclideanJordanAlgebra
1695
1696 EXAMPLES::
1697
1698 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1699 ....: [-2, 1, -4, 3],
1700 ....: [ 9, 10, 11, 12],
1701 ....: [-10, 9, -12, 11] ])
1702 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1703 [ 2*I + 1 4*I + 3]
1704 [ 10*I + 9 12*I + 11]
1705
1706 TESTS:
1707
1708 Unembedding is the inverse of embedding::
1709
1710 sage: set_random_seed()
1711 sage: F = QuadraticField(-1, 'I')
1712 sage: M = random_matrix(F, 3)
1713 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1714 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1715 True
1716
1717 """
1718 n = ZZ(M.nrows())
1719 if M.ncols() != n:
1720 raise ValueError("the matrix 'M' must be square")
1721 if not n.mod(2).is_zero():
1722 raise ValueError("the matrix 'M' must be a complex embedding")
1723
1724 # If "M" was normalized, its base ring might have roots
1725 # adjoined and they can stick around after unembedding.
1726 field = M.base_ring()
1727 R = PolynomialRing(field, 'z')
1728 z = R.gen()
1729 if field is AA:
1730 # Sage doesn't know how to embed AA into QQbar, i.e. how
1731 # to adjoin sqrt(-1) to AA.
1732 F = QQbar
1733 else:
1734 F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
1735 i = F.gen()
1736
1737 # Go top-left to bottom-right (reading order), converting every
1738 # 2-by-2 block we see to a single complex element.
1739 elements = []
1740 for k in range(n/2):
1741 for j in range(n/2):
1742 submat = M[2*k:2*k+2,2*j:2*j+2]
1743 if submat[0,0] != submat[1,1]:
1744 raise ValueError('bad on-diagonal submatrix')
1745 if submat[0,1] != -submat[1,0]:
1746 raise ValueError('bad off-diagonal submatrix')
1747 z = submat[0,0] + submat[0,1]*i
1748 elements.append(z)
1749
1750 return matrix(F, n/2, elements)
1751
1752
1753 @classmethod
1754 def trace_inner_product(cls,X,Y):
1755 """
1756 Compute a matrix inner product in this algebra directly from
1757 its real embedding.
1758
1759 SETUP::
1760
1761 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1762
1763 TESTS:
1764
1765 This gives the same answer as the slow, default method implemented
1766 in :class:`MatrixEuclideanJordanAlgebra`::
1767
1768 sage: set_random_seed()
1769 sage: J = ComplexHermitianEJA.random_instance()
1770 sage: x,y = J.random_elements(2)
1771 sage: Xe = x.to_matrix()
1772 sage: Ye = y.to_matrix()
1773 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1774 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1775 sage: expected = (X*Y).trace().real()
1776 sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye)
1777 sage: actual == expected
1778 True
1779
1780 """
1781 return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2
1782
1783
1784 class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
1785 ComplexMatrixEuclideanJordanAlgebra):
1786 """
1787 The rank-n simple EJA consisting of complex Hermitian n-by-n
1788 matrices over the real numbers, the usual symmetric Jordan product,
1789 and the real-part-of-trace inner product. It has dimension `n^2` over
1790 the reals.
1791
1792 SETUP::
1793
1794 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1795
1796 EXAMPLES:
1797
1798 In theory, our "field" can be any subfield of the reals::
1799
1800 sage: ComplexHermitianEJA(2, RDF)
1801 Euclidean Jordan algebra of dimension 4 over Real Double Field
1802 sage: ComplexHermitianEJA(2, RR)
1803 Euclidean Jordan algebra of dimension 4 over Real Field with
1804 53 bits of precision
1805
1806 TESTS:
1807
1808 The dimension of this algebra is `n^2`::
1809
1810 sage: set_random_seed()
1811 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
1812 sage: n = ZZ.random_element(1, n_max)
1813 sage: J = ComplexHermitianEJA(n)
1814 sage: J.dimension() == n^2
1815 True
1816
1817 The Jordan multiplication is what we think it is::
1818
1819 sage: set_random_seed()
1820 sage: J = ComplexHermitianEJA.random_instance()
1821 sage: x,y = J.random_elements(2)
1822 sage: actual = (x*y).to_matrix()
1823 sage: X = x.to_matrix()
1824 sage: Y = y.to_matrix()
1825 sage: expected = (X*Y + Y*X)/2
1826 sage: actual == expected
1827 True
1828 sage: J(expected) == x*y
1829 True
1830
1831 We can change the generator prefix::
1832
1833 sage: ComplexHermitianEJA(2, prefix='z').gens()
1834 (z0, z1, z2, z3)
1835
1836 We can construct the (trivial) algebra of rank zero::
1837
1838 sage: ComplexHermitianEJA(0)
1839 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1840
1841 """
1842
1843 @classmethod
1844 def _denormalized_basis(cls, n, field):
1845 """
1846 Returns a basis for the space of complex Hermitian n-by-n matrices.
1847
1848 Why do we embed these? Basically, because all of numerical linear
1849 algebra assumes that you're working with vectors consisting of `n`
1850 entries from a field and scalars from the same field. There's no way
1851 to tell SageMath that (for example) the vectors contain complex
1852 numbers, while the scalar field is real.
1853
1854 SETUP::
1855
1856 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1857
1858 TESTS::
1859
1860 sage: set_random_seed()
1861 sage: n = ZZ.random_element(1,5)
1862 sage: field = QuadraticField(2, 'sqrt2')
1863 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1864 sage: all( M.is_symmetric() for M in B)
1865 True
1866
1867 """
1868 R = PolynomialRing(field, 'z')
1869 z = R.gen()
1870 F = field.extension(z**2 + 1, 'I')
1871 I = F.gen()
1872
1873 # This is like the symmetric case, but we need to be careful:
1874 #
1875 # * We want conjugate-symmetry, not just symmetry.
1876 # * The diagonal will (as a result) be real.
1877 #
1878 S = []
1879 for i in range(n):
1880 for j in range(i+1):
1881 Eij = matrix(F, n, lambda k,l: k==i and l==j)
1882 if i == j:
1883 Sij = cls.real_embed(Eij)
1884 S.append(Sij)
1885 else:
1886 # The second one has a minus because it's conjugated.
1887 Sij_real = cls.real_embed(Eij + Eij.transpose())
1888 S.append(Sij_real)
1889 Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
1890 S.append(Sij_imag)
1891
1892 # Since we embedded these, we can drop back to the "field" that we
1893 # started with instead of the complex extension "F".
1894 return tuple( s.change_ring(field) for s in S )
1895
1896
1897 def __init__(self, n, field=AA, **kwargs):
1898 basis = self._denormalized_basis(n,field)
1899 super(ComplexHermitianEJA, self).__init__(field,
1900 basis,
1901 self.jordan_product,
1902 self.trace_inner_product,
1903 **kwargs)
1904 self.rank.set_cache(n)
1905 # TODO: pre-cache the identity!
1906
1907 @staticmethod
1908 def _max_random_instance_size():
1909 return 3 # Dimension 9
1910
1911 @classmethod
1912 def random_instance(cls, field=AA, **kwargs):
1913 """
1914 Return a random instance of this type of algebra.
1915 """
1916 n = ZZ.random_element(cls._max_random_instance_size() + 1)
1917 return cls(n, field, **kwargs)
1918
1919 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
1920 @staticmethod
1921 def real_embed(M):
1922 """
1923 Embed the n-by-n quaternion matrix ``M`` into the space of real
1924 matrices of size 4n-by-4n by first sending each quaternion entry `z
1925 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1926 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1927 matrix.
1928
1929 SETUP::
1930
1931 sage: from mjo.eja.eja_algebra import \
1932 ....: QuaternionMatrixEuclideanJordanAlgebra
1933
1934 EXAMPLES::
1935
1936 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1937 sage: i,j,k = Q.gens()
1938 sage: x = 1 + 2*i + 3*j + 4*k
1939 sage: M = matrix(Q, 1, [[x]])
1940 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1941 [ 1 2 3 4]
1942 [-2 1 -4 3]
1943 [-3 4 1 -2]
1944 [-4 -3 2 1]
1945
1946 Embedding is a homomorphism (isomorphism, in fact)::
1947
1948 sage: set_random_seed()
1949 sage: n = ZZ.random_element(2)
1950 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1951 sage: X = random_matrix(Q, n)
1952 sage: Y = random_matrix(Q, n)
1953 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1954 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1955 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1956 sage: Xe*Ye == XYe
1957 True
1958
1959 """
1960 quaternions = M.base_ring()
1961 n = M.nrows()
1962 if M.ncols() != n:
1963 raise ValueError("the matrix 'M' must be square")
1964
1965 F = QuadraticField(-1, 'I')
1966 i = F.gen()
1967
1968 blocks = []
1969 for z in M.list():
1970 t = z.coefficient_tuple()
1971 a = t[0]
1972 b = t[1]
1973 c = t[2]
1974 d = t[3]
1975 cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
1976 [-c + d*i, a - b*i]])
1977 realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM)
1978 blocks.append(realM)
1979
1980 # We should have real entries by now, so use the realest field
1981 # we've got for the return value.
1982 return matrix.block(quaternions.base_ring(), n, blocks)
1983
1984
1985
1986 @staticmethod
1987 def real_unembed(M):
1988 """
1989 The inverse of _embed_quaternion_matrix().
1990
1991 SETUP::
1992
1993 sage: from mjo.eja.eja_algebra import \
1994 ....: QuaternionMatrixEuclideanJordanAlgebra
1995
1996 EXAMPLES::
1997
1998 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1999 ....: [-2, 1, -4, 3],
2000 ....: [-3, 4, 1, -2],
2001 ....: [-4, -3, 2, 1]])
2002 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
2003 [1 + 2*i + 3*j + 4*k]
2004
2005 TESTS:
2006
2007 Unembedding is the inverse of embedding::
2008
2009 sage: set_random_seed()
2010 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2011 sage: M = random_matrix(Q, 3)
2012 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
2013 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
2014 True
2015
2016 """
2017 n = ZZ(M.nrows())
2018 if M.ncols() != n:
2019 raise ValueError("the matrix 'M' must be square")
2020 if not n.mod(4).is_zero():
2021 raise ValueError("the matrix 'M' must be a quaternion embedding")
2022
2023 # Use the base ring of the matrix to ensure that its entries can be
2024 # multiplied by elements of the quaternion algebra.
2025 field = M.base_ring()
2026 Q = QuaternionAlgebra(field,-1,-1)
2027 i,j,k = Q.gens()
2028
2029 # Go top-left to bottom-right (reading order), converting every
2030 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2031 # quaternion block.
2032 elements = []
2033 for l in range(n/4):
2034 for m in range(n/4):
2035 submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
2036 M[4*l:4*l+4,4*m:4*m+4] )
2037 if submat[0,0] != submat[1,1].conjugate():
2038 raise ValueError('bad on-diagonal submatrix')
2039 if submat[0,1] != -submat[1,0].conjugate():
2040 raise ValueError('bad off-diagonal submatrix')
2041 z = submat[0,0].real()
2042 z += submat[0,0].imag()*i
2043 z += submat[0,1].real()*j
2044 z += submat[0,1].imag()*k
2045 elements.append(z)
2046
2047 return matrix(Q, n/4, elements)
2048
2049
2050 @classmethod
2051 def trace_inner_product(cls,X,Y):
2052 """
2053 Compute a matrix inner product in this algebra directly from
2054 its real embedding.
2055
2056 SETUP::
2057
2058 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2059
2060 TESTS:
2061
2062 This gives the same answer as the slow, default method implemented
2063 in :class:`MatrixEuclideanJordanAlgebra`::
2064
2065 sage: set_random_seed()
2066 sage: J = QuaternionHermitianEJA.random_instance()
2067 sage: x,y = J.random_elements(2)
2068 sage: Xe = x.to_matrix()
2069 sage: Ye = y.to_matrix()
2070 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
2071 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
2072 sage: expected = (X*Y).trace().coefficient_tuple()[0]
2073 sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye)
2074 sage: actual == expected
2075 True
2076
2077 """
2078 return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4
2079
2080
2081 class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
2082 QuaternionMatrixEuclideanJordanAlgebra):
2083 r"""
2084 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2085 matrices, the usual symmetric Jordan product, and the
2086 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2087 the reals.
2088
2089 SETUP::
2090
2091 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2092
2093 EXAMPLES:
2094
2095 In theory, our "field" can be any subfield of the reals::
2096
2097 sage: QuaternionHermitianEJA(2, RDF)
2098 Euclidean Jordan algebra of dimension 6 over Real Double Field
2099 sage: QuaternionHermitianEJA(2, RR)
2100 Euclidean Jordan algebra of dimension 6 over Real Field with
2101 53 bits of precision
2102
2103 TESTS:
2104
2105 The dimension of this algebra is `2*n^2 - n`::
2106
2107 sage: set_random_seed()
2108 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2109 sage: n = ZZ.random_element(1, n_max)
2110 sage: J = QuaternionHermitianEJA(n)
2111 sage: J.dimension() == 2*(n^2) - n
2112 True
2113
2114 The Jordan multiplication is what we think it is::
2115
2116 sage: set_random_seed()
2117 sage: J = QuaternionHermitianEJA.random_instance()
2118 sage: x,y = J.random_elements(2)
2119 sage: actual = (x*y).to_matrix()
2120 sage: X = x.to_matrix()
2121 sage: Y = y.to_matrix()
2122 sage: expected = (X*Y + Y*X)/2
2123 sage: actual == expected
2124 True
2125 sage: J(expected) == x*y
2126 True
2127
2128 We can change the generator prefix::
2129
2130 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2131 (a0, a1, a2, a3, a4, a5)
2132
2133 We can construct the (trivial) algebra of rank zero::
2134
2135 sage: QuaternionHermitianEJA(0)
2136 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2137
2138 """
2139 @classmethod
2140 def _denormalized_basis(cls, n, field):
2141 """
2142 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2143
2144 Why do we embed these? Basically, because all of numerical
2145 linear algebra assumes that you're working with vectors consisting
2146 of `n` entries from a field and scalars from the same field. There's
2147 no way to tell SageMath that (for example) the vectors contain
2148 complex numbers, while the scalar field is real.
2149
2150 SETUP::
2151
2152 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2153
2154 TESTS::
2155
2156 sage: set_random_seed()
2157 sage: n = ZZ.random_element(1,5)
2158 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
2159 sage: all( M.is_symmetric() for M in B )
2160 True
2161
2162 """
2163 Q = QuaternionAlgebra(QQ,-1,-1)
2164 I,J,K = Q.gens()
2165
2166 # This is like the symmetric case, but we need to be careful:
2167 #
2168 # * We want conjugate-symmetry, not just symmetry.
2169 # * The diagonal will (as a result) be real.
2170 #
2171 S = []
2172 for i in range(n):
2173 for j in range(i+1):
2174 Eij = matrix(Q, n, lambda k,l: k==i and l==j)
2175 if i == j:
2176 Sij = cls.real_embed(Eij)
2177 S.append(Sij)
2178 else:
2179 # The second, third, and fourth ones have a minus
2180 # because they're conjugated.
2181 Sij_real = cls.real_embed(Eij + Eij.transpose())
2182 S.append(Sij_real)
2183 Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
2184 S.append(Sij_I)
2185 Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
2186 S.append(Sij_J)
2187 Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
2188 S.append(Sij_K)
2189
2190 # Since we embedded these, we can drop back to the "field" that we
2191 # started with instead of the quaternion algebra "Q".
2192 return tuple( s.change_ring(field) for s in S )
2193
2194
2195 def __init__(self, n, field=AA, **kwargs):
2196 basis = self._denormalized_basis(n,field)
2197 super(QuaternionHermitianEJA, self).__init__(field,
2198 basis,
2199 self.jordan_product,
2200 self.trace_inner_product,
2201 **kwargs)
2202 self.rank.set_cache(n)
2203 # TODO: cache one()!
2204
2205 @staticmethod
2206 def _max_random_instance_size():
2207 r"""
2208 The maximum rank of a random QuaternionHermitianEJA.
2209 """
2210 return 2 # Dimension 6
2211
2212 @classmethod
2213 def random_instance(cls, field=AA, **kwargs):
2214 """
2215 Return a random instance of this type of algebra.
2216 """
2217 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2218 return cls(n, field, **kwargs)
2219
2220
2221 class HadamardEJA(ConcreteEuclideanJordanAlgebra):
2222 """
2223 Return the Euclidean Jordan Algebra corresponding to the set
2224 `R^n` under the Hadamard product.
2225
2226 Note: this is nothing more than the Cartesian product of ``n``
2227 copies of the spin algebra. Once Cartesian product algebras
2228 are implemented, this can go.
2229
2230 SETUP::
2231
2232 sage: from mjo.eja.eja_algebra import HadamardEJA
2233
2234 EXAMPLES:
2235
2236 This multiplication table can be verified by hand::
2237
2238 sage: J = HadamardEJA(3)
2239 sage: e0,e1,e2 = J.gens()
2240 sage: e0*e0
2241 e0
2242 sage: e0*e1
2243 0
2244 sage: e0*e2
2245 0
2246 sage: e1*e1
2247 e1
2248 sage: e1*e2
2249 0
2250 sage: e2*e2
2251 e2
2252
2253 TESTS:
2254
2255 We can change the generator prefix::
2256
2257 sage: HadamardEJA(3, prefix='r').gens()
2258 (r0, r1, r2)
2259
2260 """
2261 def __init__(self, n, field=AA, **kwargs):
2262 V = VectorSpace(field, n)
2263 basis = V.basis()
2264
2265 def jordan_product(x,y):
2266 return V([ xi*yi for (xi,yi) in zip(x,y) ])
2267 def inner_product(x,y):
2268 return x.inner_product(y)
2269
2270 super(HadamardEJA, self).__init__(field,
2271 basis,
2272 jordan_product,
2273 inner_product,
2274 **kwargs)
2275 self.rank.set_cache(n)
2276
2277 if n == 0:
2278 self.one.set_cache( self.zero() )
2279 else:
2280 self.one.set_cache( sum(self.gens()) )
2281
2282 @staticmethod
2283 def _max_random_instance_size():
2284 r"""
2285 The maximum dimension of a random HadamardEJA.
2286 """
2287 return 5
2288
2289 @classmethod
2290 def random_instance(cls, field=AA, **kwargs):
2291 """
2292 Return a random instance of this type of algebra.
2293 """
2294 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2295 return cls(n, field, **kwargs)
2296
2297
2298 class BilinearFormEJA(ConcreteEuclideanJordanAlgebra):
2299 r"""
2300 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2301 with the half-trace inner product and jordan product ``x*y =
2302 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2303 a symmetric positive-definite "bilinear form" matrix. Its
2304 dimension is the size of `B`, and it has rank two in dimensions
2305 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2306 the identity matrix of order ``n``.
2307
2308 We insist that the one-by-one upper-left identity block of `B` be
2309 passed in as well so that we can be passed a matrix of size zero
2310 to construct a trivial algebra.
2311
2312 SETUP::
2313
2314 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2315 ....: JordanSpinEJA)
2316
2317 EXAMPLES:
2318
2319 When no bilinear form is specified, the identity matrix is used,
2320 and the resulting algebra is the Jordan spin algebra::
2321
2322 sage: B = matrix.identity(AA,3)
2323 sage: J0 = BilinearFormEJA(B)
2324 sage: J1 = JordanSpinEJA(3)
2325 sage: J0.multiplication_table() == J0.multiplication_table()
2326 True
2327
2328 An error is raised if the matrix `B` does not correspond to a
2329 positive-definite bilinear form::
2330
2331 sage: B = matrix.random(QQ,2,3)
2332 sage: J = BilinearFormEJA(B)
2333 Traceback (most recent call last):
2334 ...
2335 ValueError: bilinear form is not positive-definite
2336 sage: B = matrix.zero(QQ,3)
2337 sage: J = BilinearFormEJA(B)
2338 Traceback (most recent call last):
2339 ...
2340 ValueError: bilinear form is not positive-definite
2341
2342 TESTS:
2343
2344 We can create a zero-dimensional algebra::
2345
2346 sage: B = matrix.identity(AA,0)
2347 sage: J = BilinearFormEJA(B)
2348 sage: J.basis()
2349 Finite family {}
2350
2351 We can check the multiplication condition given in the Jordan, von
2352 Neumann, and Wigner paper (and also discussed on my "On the
2353 symmetry..." paper). Note that this relies heavily on the standard
2354 choice of basis, as does anything utilizing the bilinear form
2355 matrix. We opt not to orthonormalize the basis, because if we
2356 did, we would have to normalize the `s_{i}` in a similar manner::
2357
2358 sage: set_random_seed()
2359 sage: n = ZZ.random_element(5)
2360 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2361 sage: B11 = matrix.identity(QQ,1)
2362 sage: B22 = M.transpose()*M
2363 sage: B = block_matrix(2,2,[ [B11,0 ],
2364 ....: [0, B22 ] ])
2365 sage: J = BilinearFormEJA(B, orthonormalize=False)
2366 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2367 sage: V = J.vector_space()
2368 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2369 ....: for ei in eis ]
2370 sage: actual = [ sis[i]*sis[j]
2371 ....: for i in range(n-1)
2372 ....: for j in range(n-1) ]
2373 sage: expected = [ J.one() if i == j else J.zero()
2374 ....: for i in range(n-1)
2375 ....: for j in range(n-1) ]
2376 sage: actual == expected
2377 True
2378 """
2379 def __init__(self, B, field=AA, **kwargs):
2380 if not B.is_positive_definite():
2381 raise ValueError("bilinear form is not positive-definite")
2382
2383 n = B.nrows()
2384 V = VectorSpace(field, n)
2385
2386 def inner_product(x,y):
2387 return (B*x).inner_product(y)
2388
2389 def jordan_product(x,y):
2390 x0 = x[0]
2391 xbar = x[1:]
2392 y0 = y[0]
2393 ybar = y[1:]
2394 z0 = inner_product(x,y)
2395 zbar = y0*xbar + x0*ybar
2396 return V([z0] + zbar.list())
2397
2398 super(BilinearFormEJA, self).__init__(field,
2399 V.basis(),
2400 jordan_product,
2401 inner_product,
2402 **kwargs)
2403
2404 # The rank of this algebra is two, unless we're in a
2405 # one-dimensional ambient space (because the rank is bounded
2406 # by the ambient dimension).
2407 self.rank.set_cache(min(n,2))
2408
2409 if n == 0:
2410 self.one.set_cache( self.zero() )
2411 else:
2412 self.one.set_cache( self.monomial(0) )
2413
2414 @staticmethod
2415 def _max_random_instance_size():
2416 r"""
2417 The maximum dimension of a random BilinearFormEJA.
2418 """
2419 return 5
2420
2421 @classmethod
2422 def random_instance(cls, field=AA, **kwargs):
2423 """
2424 Return a random instance of this algebra.
2425 """
2426 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2427 if n.is_zero():
2428 B = matrix.identity(field, n)
2429 return cls(B, field, **kwargs)
2430
2431 B11 = matrix.identity(field,1)
2432 M = matrix.random(field, n-1)
2433 I = matrix.identity(field, n-1)
2434 alpha = field.zero()
2435 while alpha.is_zero():
2436 alpha = field.random_element().abs()
2437 B22 = M.transpose()*M + alpha*I
2438
2439 from sage.matrix.special import block_matrix
2440 B = block_matrix(2,2, [ [B11, ZZ(0) ],
2441 [ZZ(0), B22 ] ])
2442
2443 return cls(B, field, **kwargs)
2444
2445
2446 class JordanSpinEJA(BilinearFormEJA):
2447 """
2448 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2449 with the usual inner product and jordan product ``x*y =
2450 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2451 the reals.
2452
2453 SETUP::
2454
2455 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2456
2457 EXAMPLES:
2458
2459 This multiplication table can be verified by hand::
2460
2461 sage: J = JordanSpinEJA(4)
2462 sage: e0,e1,e2,e3 = J.gens()
2463 sage: e0*e0
2464 e0
2465 sage: e0*e1
2466 e1
2467 sage: e0*e2
2468 e2
2469 sage: e0*e3
2470 e3
2471 sage: e1*e2
2472 0
2473 sage: e1*e3
2474 0
2475 sage: e2*e3
2476 0
2477
2478 We can change the generator prefix::
2479
2480 sage: JordanSpinEJA(2, prefix='B').gens()
2481 (B0, B1)
2482
2483 TESTS:
2484
2485 Ensure that we have the usual inner product on `R^n`::
2486
2487 sage: set_random_seed()
2488 sage: J = JordanSpinEJA.random_instance()
2489 sage: x,y = J.random_elements(2)
2490 sage: actual = x.inner_product(y)
2491 sage: expected = x.to_vector().inner_product(y.to_vector())
2492 sage: actual == expected
2493 True
2494
2495 """
2496 def __init__(self, n, field=AA, **kwargs):
2497 # This is a special case of the BilinearFormEJA with the identity
2498 # matrix as its bilinear form.
2499 B = matrix.identity(field, n)
2500 super(JordanSpinEJA, self).__init__(B, field, **kwargs)
2501
2502 @staticmethod
2503 def _max_random_instance_size():
2504 r"""
2505 The maximum dimension of a random JordanSpinEJA.
2506 """
2507 return 5
2508
2509 @classmethod
2510 def random_instance(cls, field=AA, **kwargs):
2511 """
2512 Return a random instance of this type of algebra.
2513
2514 Needed here to override the implementation for ``BilinearFormEJA``.
2515 """
2516 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2517 return cls(n, field, **kwargs)
2518
2519
2520 class TrivialEJA(ConcreteEuclideanJordanAlgebra):
2521 """
2522 The trivial Euclidean Jordan algebra consisting of only a zero element.
2523
2524 SETUP::
2525
2526 sage: from mjo.eja.eja_algebra import TrivialEJA
2527
2528 EXAMPLES::
2529
2530 sage: J = TrivialEJA()
2531 sage: J.dimension()
2532 0
2533 sage: J.zero()
2534 0
2535 sage: J.one()
2536 0
2537 sage: 7*J.one()*12*J.one()
2538 0
2539 sage: J.one().inner_product(J.one())
2540 0
2541 sage: J.one().norm()
2542 0
2543 sage: J.one().subalgebra_generated_by()
2544 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2545 sage: J.rank()
2546 0
2547
2548 """
2549 def __init__(self, field=AA, **kwargs):
2550 jordan_product = lambda x,y: x
2551 inner_product = lambda x,y: field(0)
2552 basis = ()
2553 super(TrivialEJA, self).__init__(field,
2554 basis,
2555 jordan_product,
2556 inner_product,
2557 **kwargs)
2558 # The rank is zero using my definition, namely the dimension of the
2559 # largest subalgebra generated by any element.
2560 self.rank.set_cache(0)
2561 self.one.set_cache( self.zero() )
2562
2563 @classmethod
2564 def random_instance(cls, field=AA, **kwargs):
2565 # We don't take a "size" argument so the superclass method is
2566 # inappropriate for us.
2567 return cls(field, **kwargs)
2568
2569 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
2570 r"""
2571 The external (orthogonal) direct sum of two other Euclidean Jordan
2572 algebras. Essentially the Cartesian product of its two factors.
2573 Every Euclidean Jordan algebra decomposes into an orthogonal
2574 direct sum of simple Euclidean Jordan algebras, so no generality
2575 is lost by providing only this construction.
2576
2577 SETUP::
2578
2579 sage: from mjo.eja.eja_algebra import (random_eja,
2580 ....: HadamardEJA,
2581 ....: RealSymmetricEJA,
2582 ....: DirectSumEJA)
2583
2584 EXAMPLES::
2585
2586 sage: J1 = HadamardEJA(2)
2587 sage: J2 = RealSymmetricEJA(3)
2588 sage: J = DirectSumEJA(J1,J2)
2589 sage: J.dimension()
2590 8
2591 sage: J.rank()
2592 5
2593
2594 TESTS:
2595
2596 The external direct sum construction is only valid when the two factors
2597 have the same base ring; an error is raised otherwise::
2598
2599 sage: set_random_seed()
2600 sage: J1 = random_eja(AA)
2601 sage: J2 = random_eja(QQ,orthonormalize=False)
2602 sage: J = DirectSumEJA(J1,J2)
2603 Traceback (most recent call last):
2604 ...
2605 ValueError: algebras must share the same base field
2606
2607 """
2608 def __init__(self, J1, J2, **kwargs):
2609 if J1.base_ring() != J2.base_ring():
2610 raise ValueError("algebras must share the same base field")
2611 field = J1.base_ring()
2612
2613 self._factors = (J1, J2)
2614 n1 = J1.dimension()
2615 n2 = J2.dimension()
2616 n = n1+n2
2617 V = VectorSpace(field, n)
2618 mult_table = [ [ V.zero() for j in range(i+1) ]
2619 for i in range(n) ]
2620 for i in range(n1):
2621 for j in range(i+1):
2622 p = (J1.monomial(i)*J1.monomial(j)).to_vector()
2623 mult_table[i][j] = V(p.list() + [field.zero()]*n2)
2624
2625 for i in range(n2):
2626 for j in range(i+1):
2627 p = (J2.monomial(i)*J2.monomial(j)).to_vector()
2628 mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list())
2629
2630 # TODO: build the IP table here from the two constituent IP
2631 # matrices (it'll be block diagonal, I think).
2632 ip_table = [ [ field.zero() for j in range(i+1) ]
2633 for i in range(n) ]
2634 super(DirectSumEJA, self).__init__(field,
2635 mult_table,
2636 ip_table,
2637 check_axioms=False,
2638 **kwargs)
2639 self.rank.set_cache(J1.rank() + J2.rank())
2640
2641
2642 def factors(self):
2643 r"""
2644 Return the pair of this algebra's factors.
2645
2646 SETUP::
2647
2648 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2649 ....: JordanSpinEJA,
2650 ....: DirectSumEJA)
2651
2652 EXAMPLES::
2653
2654 sage: J1 = HadamardEJA(2,QQ)
2655 sage: J2 = JordanSpinEJA(3,QQ)
2656 sage: J = DirectSumEJA(J1,J2)
2657 sage: J.factors()
2658 (Euclidean Jordan algebra of dimension 2 over Rational Field,
2659 Euclidean Jordan algebra of dimension 3 over Rational Field)
2660
2661 """
2662 return self._factors
2663
2664 def projections(self):
2665 r"""
2666 Return a pair of projections onto this algebra's factors.
2667
2668 SETUP::
2669
2670 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
2671 ....: ComplexHermitianEJA,
2672 ....: DirectSumEJA)
2673
2674 EXAMPLES::
2675
2676 sage: J1 = JordanSpinEJA(2)
2677 sage: J2 = ComplexHermitianEJA(2)
2678 sage: J = DirectSumEJA(J1,J2)
2679 sage: (pi_left, pi_right) = J.projections()
2680 sage: J.one().to_vector()
2681 (1, 0, 1, 0, 0, 1)
2682 sage: pi_left(J.one()).to_vector()
2683 (1, 0)
2684 sage: pi_right(J.one()).to_vector()
2685 (1, 0, 0, 1)
2686
2687 """
2688 (J1,J2) = self.factors()
2689 m = J1.dimension()
2690 n = J2.dimension()
2691 V_basis = self.vector_space().basis()
2692 # Need to specify the dimensions explicitly so that we don't
2693 # wind up with a zero-by-zero matrix when we want e.g. a
2694 # zero-by-two matrix (important for composing things).
2695 P1 = matrix(self.base_ring(), m, m+n, V_basis[:m])
2696 P2 = matrix(self.base_ring(), n, m+n, V_basis[m:])
2697 pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1)
2698 pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2)
2699 return (pi_left, pi_right)
2700
2701 def inclusions(self):
2702 r"""
2703 Return the pair of inclusion maps from our factors into us.
2704
2705 SETUP::
2706
2707 sage: from mjo.eja.eja_algebra import (random_eja,
2708 ....: JordanSpinEJA,
2709 ....: RealSymmetricEJA,
2710 ....: DirectSumEJA)
2711
2712 EXAMPLES::
2713
2714 sage: J1 = JordanSpinEJA(3)
2715 sage: J2 = RealSymmetricEJA(2)
2716 sage: J = DirectSumEJA(J1,J2)
2717 sage: (iota_left, iota_right) = J.inclusions()
2718 sage: iota_left(J1.zero()) == J.zero()
2719 True
2720 sage: iota_right(J2.zero()) == J.zero()
2721 True
2722 sage: J1.one().to_vector()
2723 (1, 0, 0)
2724 sage: iota_left(J1.one()).to_vector()
2725 (1, 0, 0, 0, 0, 0)
2726 sage: J2.one().to_vector()
2727 (1, 0, 1)
2728 sage: iota_right(J2.one()).to_vector()
2729 (0, 0, 0, 1, 0, 1)
2730 sage: J.one().to_vector()
2731 (1, 0, 0, 1, 0, 1)
2732
2733 TESTS:
2734
2735 Composing a projection with the corresponding inclusion should
2736 produce the identity map, and mismatching them should produce
2737 the zero map::
2738
2739 sage: set_random_seed()
2740 sage: J1 = random_eja()
2741 sage: J2 = random_eja()
2742 sage: J = DirectSumEJA(J1,J2)
2743 sage: (iota_left, iota_right) = J.inclusions()
2744 sage: (pi_left, pi_right) = J.projections()
2745 sage: pi_left*iota_left == J1.one().operator()
2746 True
2747 sage: pi_right*iota_right == J2.one().operator()
2748 True
2749 sage: (pi_left*iota_right).is_zero()
2750 True
2751 sage: (pi_right*iota_left).is_zero()
2752 True
2753
2754 """
2755 (J1,J2) = self.factors()
2756 m = J1.dimension()
2757 n = J2.dimension()
2758 V_basis = self.vector_space().basis()
2759 # Need to specify the dimensions explicitly so that we don't
2760 # wind up with a zero-by-zero matrix when we want e.g. a
2761 # two-by-zero matrix (important for composing things).
2762 I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m])
2763 I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:])
2764 iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1)
2765 iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2)
2766 return (iota_left, iota_right)
2767
2768 def inner_product(self, x, y):
2769 r"""
2770 The standard Cartesian inner-product.
2771
2772 We project ``x`` and ``y`` onto our factors, and add up the
2773 inner-products from the subalgebras.
2774
2775 SETUP::
2776
2777
2778 sage: from mjo.eja.eja_algebra import (HadamardEJA,
2779 ....: QuaternionHermitianEJA,
2780 ....: DirectSumEJA)
2781
2782 EXAMPLE::
2783
2784 sage: J1 = HadamardEJA(3,QQ)
2785 sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False)
2786 sage: J = DirectSumEJA(J1,J2)
2787 sage: x1 = J1.one()
2788 sage: x2 = x1
2789 sage: y1 = J2.one()
2790 sage: y2 = y1
2791 sage: x1.inner_product(x2)
2792 3
2793 sage: y1.inner_product(y2)
2794 2
2795 sage: J.one().inner_product(J.one())
2796 5
2797
2798 """
2799 (pi_left, pi_right) = self.projections()
2800 x1 = pi_left(x)
2801 x2 = pi_right(x)
2802 y1 = pi_left(y)
2803 y2 = pi_right(y)
2804
2805 return (x1.inner_product(y1) + x2.inner_product(y2))
2806
2807
2808
2809 random_eja = ConcreteEuclideanJordanAlgebra.random_instance