]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/eja/eja_algebra.py
eja: fix an octonion EJA test.
[sage.d.git] / mjo / eja / eja_algebra.py
1 """
2 Representations and constructions for Euclidean Jordan algebras.
3
4 A Euclidean Jordan algebra is a Jordan algebra that has some
5 additional properties:
6
7 1. It is finite-dimensional.
8 2. Its scalar field is the real numbers.
9 3a. An inner product is defined on it, and...
10 3b. That inner product is compatible with the Jordan product
11 in the sense that `<x*y,z> = <y,x*z>` for all elements
12 `x,y,z` in the algebra.
13
14 Every Euclidean Jordan algebra is formally-real: for any two elements
15 `x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
16 0`. Conversely, every finite-dimensional formally-real Jordan algebra
17 can be made into a Euclidean Jordan algebra with an appropriate choice
18 of inner-product.
19
20 Formally-real Jordan algebras were originally studied as a framework
21 for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
22 symmetric cone optimization, since every symmetric cone arises as the
23 cone of squares in some Euclidean Jordan algebra.
24
25 It is known that every Euclidean Jordan algebra decomposes into an
26 orthogonal direct sum (essentially, a Cartesian product) of simple
27 algebras, and that moreover, up to Jordan-algebra isomorphism, there
28 are only five families of simple algebras. We provide constructions
29 for these simple algebras:
30
31 * :class:`BilinearFormEJA`
32 * :class:`RealSymmetricEJA`
33 * :class:`ComplexHermitianEJA`
34 * :class:`QuaternionHermitianEJA`
35
36 Missing from this list is the algebra of three-by-three octononion
37 Hermitian matrices, as there is (as of yet) no implementation of the
38 octonions in SageMath. In addition to these, we provide two other
39 example constructions,
40
41 * :class:`HadamardEJA`
42 * :class:`TrivialEJA`
43
44 The Jordan spin algebra is a bilinear form algebra where the bilinear
45 form is the identity. The Hadamard EJA is simply a Cartesian product
46 of one-dimensional spin algebras. And last but not least, the trivial
47 EJA is exactly what you think. Cartesian products of these are also
48 supported using the usual ``cartesian_product()`` function; as a
49 result, we support (up to isomorphism) all Euclidean Jordan algebras
50 that don't involve octonions.
51
52 SETUP::
53
54 sage: from mjo.eja.eja_algebra import random_eja
55
56 EXAMPLES::
57
58 sage: random_eja()
59 Euclidean Jordan algebra of dimension...
60 """
61
62 from itertools import repeat
63
64 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
65 from sage.categories.magmatic_algebras import MagmaticAlgebras
66 from sage.categories.sets_cat import cartesian_product
67 from sage.combinat.free_module import CombinatorialFreeModule
68 from sage.matrix.constructor import matrix
69 from sage.matrix.matrix_space import MatrixSpace
70 from sage.misc.cachefunc import cached_method
71 from sage.misc.table import table
72 from sage.modules.free_module import FreeModule, VectorSpace
73 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
74 PolynomialRing,
75 QuadraticField)
76 from mjo.eja.eja_element import FiniteDimensionalEJAElement
77 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
78 from mjo.eja.eja_utils import _all2list, _mat2vec
79
80 class FiniteDimensionalEJA(CombinatorialFreeModule):
81 r"""
82 A finite-dimensional Euclidean Jordan algebra.
83
84 INPUT:
85
86 - ``basis`` -- a tuple; a tuple of basis elements in "matrix
87 form," which must be the same form as the arguments to
88 ``jordan_product`` and ``inner_product``. In reality, "matrix
89 form" can be either vectors, matrices, or a Cartesian product
90 (ordered tuple) of vectors or matrices. All of these would
91 ideally be vector spaces in sage with no special-casing
92 needed; but in reality we turn vectors into column-matrices
93 and Cartesian products `(a,b)` into column matrices
94 `(a,b)^{T}` after converting `a` and `b` themselves.
95
96 - ``jordan_product`` -- a function; afunction of two ``basis``
97 elements (in matrix form) that returns their jordan product,
98 also in matrix form; this will be applied to ``basis`` to
99 compute a multiplication table for the algebra.
100
101 - ``inner_product`` -- a function; a function of two ``basis``
102 elements (in matrix form) that returns their inner
103 product. This will be applied to ``basis`` to compute an
104 inner-product table (basically a matrix) for this algebra.
105
106 - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
107 field for the algebra.
108
109 - ``orthonormalize`` -- boolean (default: ``True``); whether or
110 not to orthonormalize the basis. Doing so is expensive and
111 generally rules out using the rationals as your ``field``, but
112 is required for spectral decompositions.
113
114 SETUP::
115
116 sage: from mjo.eja.eja_algebra import random_eja
117
118 TESTS:
119
120 We should compute that an element subalgebra is associative even
121 if we circumvent the element method::
122
123 sage: set_random_seed()
124 sage: J = random_eja(field=QQ,orthonormalize=False)
125 sage: x = J.random_element()
126 sage: A = x.subalgebra_generated_by(orthonormalize=False)
127 sage: basis = tuple(b.superalgebra_element() for b in A.basis())
128 sage: J.subalgebra(basis, orthonormalize=False).is_associative()
129 True
130
131 """
132 Element = FiniteDimensionalEJAElement
133
134 def __init__(self,
135 basis,
136 jordan_product,
137 inner_product,
138 field=AA,
139 orthonormalize=True,
140 associative=None,
141 cartesian_product=False,
142 check_field=True,
143 check_axioms=True,
144 prefix="b"):
145
146 n = len(basis)
147
148 if check_field:
149 if not field.is_subring(RR):
150 # Note: this does return true for the real algebraic
151 # field, the rationals, and any quadratic field where
152 # we've specified a real embedding.
153 raise ValueError("scalar field is not real")
154
155 if check_axioms:
156 # Check commutativity of the Jordan and inner-products.
157 # This has to be done before we build the multiplication
158 # and inner-product tables/matrices, because we take
159 # advantage of symmetry in the process.
160 if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
161 for bi in basis
162 for bj in basis ):
163 raise ValueError("Jordan product is not commutative")
164
165 if not all( inner_product(bi,bj) == inner_product(bj,bi)
166 for bi in basis
167 for bj in basis ):
168 raise ValueError("inner-product is not commutative")
169
170
171 category = MagmaticAlgebras(field).FiniteDimensional()
172 category = category.WithBasis().Unital().Commutative()
173
174 if associative is None:
175 # We should figure it out. As with check_axioms, we have to do
176 # this without the help of the _jordan_product_is_associative()
177 # method because we need to know the category before we
178 # initialize the algebra.
179 associative = all( jordan_product(jordan_product(bi,bj),bk)
180 ==
181 jordan_product(bi,jordan_product(bj,bk))
182 for bi in basis
183 for bj in basis
184 for bk in basis)
185
186 if associative:
187 # Element subalgebras can take advantage of this.
188 category = category.Associative()
189 if cartesian_product:
190 # Use join() here because otherwise we only get the
191 # "Cartesian product of..." and not the things themselves.
192 category = category.join([category,
193 category.CartesianProducts()])
194
195 # Call the superclass constructor so that we can use its from_vector()
196 # method to build our multiplication table.
197 CombinatorialFreeModule.__init__(self,
198 field,
199 range(n),
200 prefix=prefix,
201 category=category,
202 bracket=False)
203
204 # Now comes all of the hard work. We'll be constructing an
205 # ambient vector space V that our (vectorized) basis lives in,
206 # as well as a subspace W of V spanned by those (vectorized)
207 # basis elements. The W-coordinates are the coefficients that
208 # we see in things like x = 1*b1 + 2*b2.
209 vector_basis = basis
210
211 degree = 0
212 if n > 0:
213 degree = len(_all2list(basis[0]))
214
215 # Build an ambient space that fits our matrix basis when
216 # written out as "long vectors."
217 V = VectorSpace(field, degree)
218
219 # The matrix that will hole the orthonormal -> unorthonormal
220 # coordinate transformation.
221 self._deortho_matrix = None
222
223 if orthonormalize:
224 # Save a copy of the un-orthonormalized basis for later.
225 # Convert it to ambient V (vector) coordinates while we're
226 # at it, because we'd have to do it later anyway.
227 deortho_vector_basis = tuple( V(_all2list(b)) for b in basis )
228
229 from mjo.eja.eja_utils import gram_schmidt
230 basis = tuple(gram_schmidt(basis, inner_product))
231
232 # Save the (possibly orthonormalized) matrix basis for
233 # later...
234 self._matrix_basis = basis
235
236 # Now create the vector space for the algebra, which will have
237 # its own set of non-ambient coordinates (in terms of the
238 # supplied basis).
239 vector_basis = tuple( V(_all2list(b)) for b in basis )
240 W = V.span_of_basis( vector_basis, check=check_axioms)
241
242 if orthonormalize:
243 # Now "W" is the vector space of our algebra coordinates. The
244 # variables "X1", "X2",... refer to the entries of vectors in
245 # W. Thus to convert back and forth between the orthonormal
246 # coordinates and the given ones, we need to stick the original
247 # basis in W.
248 U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
249 self._deortho_matrix = matrix( U.coordinate_vector(q)
250 for q in vector_basis )
251
252
253 # Now we actually compute the multiplication and inner-product
254 # tables/matrices using the possibly-orthonormalized basis.
255 self._inner_product_matrix = matrix.identity(field, n)
256 self._multiplication_table = [ [0 for j in range(i+1)]
257 for i in range(n) ]
258
259 # Note: the Jordan and inner-products are defined in terms
260 # of the ambient basis. It's important that their arguments
261 # are in ambient coordinates as well.
262 for i in range(n):
263 for j in range(i+1):
264 # ortho basis w.r.t. ambient coords
265 q_i = basis[i]
266 q_j = basis[j]
267
268 # The jordan product returns a matrixy answer, so we
269 # have to convert it to the algebra coordinates.
270 elt = jordan_product(q_i, q_j)
271 elt = W.coordinate_vector(V(_all2list(elt)))
272 self._multiplication_table[i][j] = self.from_vector(elt)
273
274 if not orthonormalize:
275 # If we're orthonormalizing the basis with respect
276 # to an inner-product, then the inner-product
277 # matrix with respect to the resulting basis is
278 # just going to be the identity.
279 ip = inner_product(q_i, q_j)
280 self._inner_product_matrix[i,j] = ip
281 self._inner_product_matrix[j,i] = ip
282
283 self._inner_product_matrix._cache = {'hermitian': True}
284 self._inner_product_matrix.set_immutable()
285
286 if check_axioms:
287 if not self._is_jordanian():
288 raise ValueError("Jordan identity does not hold")
289 if not self._inner_product_is_associative():
290 raise ValueError("inner product is not associative")
291
292
293 def _coerce_map_from_base_ring(self):
294 """
295 Disable the map from the base ring into the algebra.
296
297 Performing a nonsense conversion like this automatically
298 is counterpedagogical. The fallback is to try the usual
299 element constructor, which should also fail.
300
301 SETUP::
302
303 sage: from mjo.eja.eja_algebra import random_eja
304
305 TESTS::
306
307 sage: set_random_seed()
308 sage: J = random_eja()
309 sage: J(1)
310 Traceback (most recent call last):
311 ...
312 ValueError: not an element of this algebra
313
314 """
315 return None
316
317
318 def product_on_basis(self, i, j):
319 r"""
320 Returns the Jordan product of the `i` and `j`th basis elements.
321
322 This completely defines the Jordan product on the algebra, and
323 is used direclty by our superclass machinery to implement
324 :meth:`product`.
325
326 SETUP::
327
328 sage: from mjo.eja.eja_algebra import random_eja
329
330 TESTS::
331
332 sage: set_random_seed()
333 sage: J = random_eja()
334 sage: n = J.dimension()
335 sage: bi = J.zero()
336 sage: bj = J.zero()
337 sage: bi_bj = J.zero()*J.zero()
338 sage: if n > 0:
339 ....: i = ZZ.random_element(n)
340 ....: j = ZZ.random_element(n)
341 ....: bi = J.monomial(i)
342 ....: bj = J.monomial(j)
343 ....: bi_bj = J.product_on_basis(i,j)
344 sage: bi*bj == bi_bj
345 True
346
347 """
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 inner_product(self, x, y):
356 """
357 The inner product associated with this Euclidean Jordan algebra.
358
359 Defaults to the trace inner product, but can be overridden by
360 subclasses if they are sure that the necessary properties are
361 satisfied.
362
363 SETUP::
364
365 sage: from mjo.eja.eja_algebra import (random_eja,
366 ....: HadamardEJA,
367 ....: BilinearFormEJA)
368
369 EXAMPLES:
370
371 Our inner product is "associative," which means the following for
372 a symmetric bilinear form::
373
374 sage: set_random_seed()
375 sage: J = random_eja()
376 sage: x,y,z = J.random_elements(3)
377 sage: (x*y).inner_product(z) == y.inner_product(x*z)
378 True
379
380 TESTS:
381
382 Ensure that this is the usual inner product for the algebras
383 over `R^n`::
384
385 sage: set_random_seed()
386 sage: J = HadamardEJA.random_instance()
387 sage: x,y = J.random_elements(2)
388 sage: actual = x.inner_product(y)
389 sage: expected = x.to_vector().inner_product(y.to_vector())
390 sage: actual == expected
391 True
392
393 Ensure that this is one-half of the trace inner-product in a
394 BilinearFormEJA that isn't just the reals (when ``n`` isn't
395 one). This is in Faraut and Koranyi, and also my "On the
396 symmetry..." paper::
397
398 sage: set_random_seed()
399 sage: J = BilinearFormEJA.random_instance()
400 sage: n = J.dimension()
401 sage: x = J.random_element()
402 sage: y = J.random_element()
403 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
404 True
405
406 """
407 B = self._inner_product_matrix
408 return (B*x.to_vector()).inner_product(y.to_vector())
409
410
411 def is_associative(self):
412 r"""
413 Return whether or not this algebra's Jordan product is associative.
414
415 SETUP::
416
417 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
418
419 EXAMPLES::
420
421 sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
422 sage: J.is_associative()
423 False
424 sage: x = sum(J.gens())
425 sage: A = x.subalgebra_generated_by(orthonormalize=False)
426 sage: A.is_associative()
427 True
428
429 """
430 return "Associative" in self.category().axioms()
431
432 def _is_commutative(self):
433 r"""
434 Whether or not this algebra's multiplication table is commutative.
435
436 This method should of course always return ``True``, unless
437 this algebra was constructed with ``check_axioms=False`` and
438 passed an invalid multiplication table.
439 """
440 return all( x*y == y*x for x in self.gens() for y in self.gens() )
441
442 def _is_jordanian(self):
443 r"""
444 Whether or not this algebra's multiplication table respects the
445 Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
446
447 We only check one arrangement of `x` and `y`, so for a
448 ``True`` result to be truly true, you should also check
449 :meth:`_is_commutative`. This method should of course always
450 return ``True``, unless this algebra was constructed with
451 ``check_axioms=False`` and passed an invalid multiplication table.
452 """
453 return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
454 ==
455 (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
456 for i in range(self.dimension())
457 for j in range(self.dimension()) )
458
459 def _jordan_product_is_associative(self):
460 r"""
461 Return whether or not this algebra's Jordan product is
462 associative; that is, whether or not `x*(y*z) = (x*y)*z`
463 for all `x,y,x`.
464
465 This method should agree with :meth:`is_associative` unless
466 you lied about the value of the ``associative`` parameter
467 when you constructed the algebra.
468
469 SETUP::
470
471 sage: from mjo.eja.eja_algebra import (random_eja,
472 ....: RealSymmetricEJA,
473 ....: ComplexHermitianEJA,
474 ....: QuaternionHermitianEJA)
475
476 EXAMPLES::
477
478 sage: J = RealSymmetricEJA(4, orthonormalize=False)
479 sage: J._jordan_product_is_associative()
480 False
481 sage: x = sum(J.gens())
482 sage: A = x.subalgebra_generated_by()
483 sage: A._jordan_product_is_associative()
484 True
485
486 ::
487
488 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
489 sage: J._jordan_product_is_associative()
490 False
491 sage: x = sum(J.gens())
492 sage: A = x.subalgebra_generated_by(orthonormalize=False)
493 sage: A._jordan_product_is_associative()
494 True
495
496 ::
497
498 sage: J = QuaternionHermitianEJA(2)
499 sage: J._jordan_product_is_associative()
500 False
501 sage: x = sum(J.gens())
502 sage: A = x.subalgebra_generated_by()
503 sage: A._jordan_product_is_associative()
504 True
505
506 TESTS:
507
508 The values we've presupplied to the constructors agree with
509 the computation::
510
511 sage: set_random_seed()
512 sage: J = random_eja()
513 sage: J.is_associative() == J._jordan_product_is_associative()
514 True
515
516 """
517 R = self.base_ring()
518
519 # Used to check whether or not something is zero.
520 epsilon = R.zero()
521 if not R.is_exact():
522 # I don't know of any examples that make this magnitude
523 # necessary because I don't know how to make an
524 # associative algebra when the element subalgebra
525 # construction is unreliable (as it is over RDF; we can't
526 # find the degree of an element because we can't compute
527 # the rank of a matrix). But even multiplication of floats
528 # is non-associative, so *some* epsilon is needed... let's
529 # just take the one from _inner_product_is_associative?
530 epsilon = 1e-15
531
532 for i in range(self.dimension()):
533 for j in range(self.dimension()):
534 for k in range(self.dimension()):
535 x = self.monomial(i)
536 y = self.monomial(j)
537 z = self.monomial(k)
538 diff = (x*y)*z - x*(y*z)
539
540 if diff.norm() > epsilon:
541 return False
542
543 return True
544
545 def _inner_product_is_associative(self):
546 r"""
547 Return whether or not this algebra's inner product `B` is
548 associative; that is, whether or not `B(xy,z) = B(x,yz)`.
549
550 This method should of course always return ``True``, unless
551 this algebra was constructed with ``check_axioms=False`` and
552 passed an invalid Jordan or inner-product.
553 """
554 R = self.base_ring()
555
556 # Used to check whether or not something is zero.
557 epsilon = R.zero()
558 if not R.is_exact():
559 # This choice is sufficient to allow the construction of
560 # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
561 epsilon = 1e-15
562
563 for i in range(self.dimension()):
564 for j in range(self.dimension()):
565 for k in range(self.dimension()):
566 x = self.monomial(i)
567 y = self.monomial(j)
568 z = self.monomial(k)
569 diff = (x*y).inner_product(z) - x.inner_product(y*z)
570
571 if diff.abs() > epsilon:
572 return False
573
574 return True
575
576 def _element_constructor_(self, elt):
577 """
578 Construct an element of this algebra from its vector or matrix
579 representation.
580
581 This gets called only after the parent element _call_ method
582 fails to find a coercion for the argument.
583
584 SETUP::
585
586 sage: from mjo.eja.eja_algebra import (random_eja,
587 ....: JordanSpinEJA,
588 ....: HadamardEJA,
589 ....: RealSymmetricEJA)
590
591 EXAMPLES:
592
593 The identity in `S^n` is converted to the identity in the EJA::
594
595 sage: J = RealSymmetricEJA(3)
596 sage: I = matrix.identity(QQ,3)
597 sage: J(I) == J.one()
598 True
599
600 This skew-symmetric matrix can't be represented in the EJA::
601
602 sage: J = RealSymmetricEJA(3)
603 sage: A = matrix(QQ,3, lambda i,j: i-j)
604 sage: J(A)
605 Traceback (most recent call last):
606 ...
607 ValueError: not an element of this algebra
608
609 Tuples work as well, provided that the matrix basis for the
610 algebra consists of them::
611
612 sage: J1 = HadamardEJA(3)
613 sage: J2 = RealSymmetricEJA(2)
614 sage: J = cartesian_product([J1,J2])
615 sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
616 b1 + b5
617
618 TESTS:
619
620 Ensure that we can convert any element back and forth
621 faithfully between its matrix and algebra representations::
622
623 sage: set_random_seed()
624 sage: J = random_eja()
625 sage: x = J.random_element()
626 sage: J(x.to_matrix()) == x
627 True
628
629 We cannot coerce elements between algebras just because their
630 matrix representations are compatible::
631
632 sage: J1 = HadamardEJA(3)
633 sage: J2 = JordanSpinEJA(3)
634 sage: J2(J1.one())
635 Traceback (most recent call last):
636 ...
637 ValueError: not an element of this algebra
638 sage: J1(J2.zero())
639 Traceback (most recent call last):
640 ...
641 ValueError: not an element of this algebra
642 """
643 msg = "not an element of this algebra"
644 if elt in self.base_ring():
645 # Ensure that no base ring -> algebra coercion is performed
646 # by this method. There's some stupidity in sage that would
647 # otherwise propagate to this method; for example, sage thinks
648 # that the integer 3 belongs to the space of 2-by-2 matrices.
649 raise ValueError(msg)
650
651 try:
652 # Try to convert a vector into a column-matrix...
653 elt = elt.column()
654 except (AttributeError, TypeError):
655 # and ignore failure, because we weren't really expecting
656 # a vector as an argument anyway.
657 pass
658
659 if elt not in self.matrix_space():
660 raise ValueError(msg)
661
662 # Thanks for nothing! Matrix spaces aren't vector spaces in
663 # Sage, so we have to figure out its matrix-basis coordinates
664 # ourselves. We use the basis space's ring instead of the
665 # element's ring because the basis space might be an algebraic
666 # closure whereas the base ring of the 3-by-3 identity matrix
667 # could be QQ instead of QQbar.
668 #
669 # And, we also have to handle Cartesian product bases (when
670 # the matrix basis consists of tuples) here. The "good news"
671 # is that we're already converting everything to long vectors,
672 # and that strategy works for tuples as well.
673 #
674 # We pass check=False because the matrix basis is "guaranteed"
675 # to be linearly independent... right? Ha ha.
676 elt = _all2list(elt)
677 V = VectorSpace(self.base_ring(), len(elt))
678 W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()),
679 check=False)
680
681 try:
682 coords = W.coordinate_vector(V(elt))
683 except ArithmeticError: # vector is not in free module
684 raise ValueError(msg)
685
686 return self.from_vector(coords)
687
688 def _repr_(self):
689 """
690 Return a string representation of ``self``.
691
692 SETUP::
693
694 sage: from mjo.eja.eja_algebra import JordanSpinEJA
695
696 TESTS:
697
698 Ensure that it says what we think it says::
699
700 sage: JordanSpinEJA(2, field=AA)
701 Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
702 sage: JordanSpinEJA(3, field=RDF)
703 Euclidean Jordan algebra of dimension 3 over Real Double Field
704
705 """
706 fmt = "Euclidean Jordan algebra of dimension {} over {}"
707 return fmt.format(self.dimension(), self.base_ring())
708
709
710 @cached_method
711 def characteristic_polynomial_of(self):
712 """
713 Return the algebra's "characteristic polynomial of" function,
714 which is itself a multivariate polynomial that, when evaluated
715 at the coordinates of some algebra element, returns that
716 element's characteristic polynomial.
717
718 The resulting polynomial has `n+1` variables, where `n` is the
719 dimension of this algebra. The first `n` variables correspond to
720 the coordinates of an algebra element: when evaluated at the
721 coordinates of an algebra element with respect to a certain
722 basis, the result is a univariate polynomial (in the one
723 remaining variable ``t``), namely the characteristic polynomial
724 of that element.
725
726 SETUP::
727
728 sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
729
730 EXAMPLES:
731
732 The characteristic polynomial in the spin algebra is given in
733 Alizadeh, Example 11.11::
734
735 sage: J = JordanSpinEJA(3)
736 sage: p = J.characteristic_polynomial_of(); p
737 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
738 sage: xvec = J.one().to_vector()
739 sage: p(*xvec)
740 t^2 - 2*t + 1
741
742 By definition, the characteristic polynomial is a monic
743 degree-zero polynomial in a rank-zero algebra. Note that
744 Cayley-Hamilton is indeed satisfied since the polynomial
745 ``1`` evaluates to the identity element of the algebra on
746 any argument::
747
748 sage: J = TrivialEJA()
749 sage: J.characteristic_polynomial_of()
750 1
751
752 """
753 r = self.rank()
754 n = self.dimension()
755
756 # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
757 a = self._charpoly_coefficients()
758
759 # We go to a bit of trouble here to reorder the
760 # indeterminates, so that it's easier to evaluate the
761 # characteristic polynomial at x's coordinates and get back
762 # something in terms of t, which is what we want.
763 S = PolynomialRing(self.base_ring(),'t')
764 t = S.gen(0)
765 if r > 0:
766 R = a[0].parent()
767 S = PolynomialRing(S, R.variable_names())
768 t = S(t)
769
770 return (t**r + sum( a[k]*(t**k) for k in range(r) ))
771
772 def coordinate_polynomial_ring(self):
773 r"""
774 The multivariate polynomial ring in which this algebra's
775 :meth:`characteristic_polynomial_of` lives.
776
777 SETUP::
778
779 sage: from mjo.eja.eja_algebra import (HadamardEJA,
780 ....: RealSymmetricEJA)
781
782 EXAMPLES::
783
784 sage: J = HadamardEJA(2)
785 sage: J.coordinate_polynomial_ring()
786 Multivariate Polynomial Ring in X1, X2...
787 sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
788 sage: J.coordinate_polynomial_ring()
789 Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
790
791 """
792 var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) )
793 return PolynomialRing(self.base_ring(), var_names)
794
795 def inner_product(self, x, y):
796 """
797 The inner product associated with this Euclidean Jordan algebra.
798
799 Defaults to the trace inner product, but can be overridden by
800 subclasses if they are sure that the necessary properties are
801 satisfied.
802
803 SETUP::
804
805 sage: from mjo.eja.eja_algebra import (random_eja,
806 ....: HadamardEJA,
807 ....: BilinearFormEJA)
808
809 EXAMPLES:
810
811 Our inner product is "associative," which means the following for
812 a symmetric bilinear form::
813
814 sage: set_random_seed()
815 sage: J = random_eja()
816 sage: x,y,z = J.random_elements(3)
817 sage: (x*y).inner_product(z) == y.inner_product(x*z)
818 True
819
820 TESTS:
821
822 Ensure that this is the usual inner product for the algebras
823 over `R^n`::
824
825 sage: set_random_seed()
826 sage: J = HadamardEJA.random_instance()
827 sage: x,y = J.random_elements(2)
828 sage: actual = x.inner_product(y)
829 sage: expected = x.to_vector().inner_product(y.to_vector())
830 sage: actual == expected
831 True
832
833 Ensure that this is one-half of the trace inner-product in a
834 BilinearFormEJA that isn't just the reals (when ``n`` isn't
835 one). This is in Faraut and Koranyi, and also my "On the
836 symmetry..." paper::
837
838 sage: set_random_seed()
839 sage: J = BilinearFormEJA.random_instance()
840 sage: n = J.dimension()
841 sage: x = J.random_element()
842 sage: y = J.random_element()
843 sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
844 True
845 """
846 B = self._inner_product_matrix
847 return (B*x.to_vector()).inner_product(y.to_vector())
848
849
850 def is_trivial(self):
851 """
852 Return whether or not this algebra is trivial.
853
854 A trivial algebra contains only the zero element.
855
856 SETUP::
857
858 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
859 ....: TrivialEJA)
860
861 EXAMPLES::
862
863 sage: J = ComplexHermitianEJA(3)
864 sage: J.is_trivial()
865 False
866
867 ::
868
869 sage: J = TrivialEJA()
870 sage: J.is_trivial()
871 True
872
873 """
874 return self.dimension() == 0
875
876
877 def multiplication_table(self):
878 """
879 Return a visual representation of this algebra's multiplication
880 table (on basis elements).
881
882 SETUP::
883
884 sage: from mjo.eja.eja_algebra import JordanSpinEJA
885
886 EXAMPLES::
887
888 sage: J = JordanSpinEJA(4)
889 sage: J.multiplication_table()
890 +----++----+----+----+----+
891 | * || b0 | b1 | b2 | b3 |
892 +====++====+====+====+====+
893 | b0 || b0 | b1 | b2 | b3 |
894 +----++----+----+----+----+
895 | b1 || b1 | b0 | 0 | 0 |
896 +----++----+----+----+----+
897 | b2 || b2 | 0 | b0 | 0 |
898 +----++----+----+----+----+
899 | b3 || b3 | 0 | 0 | b0 |
900 +----++----+----+----+----+
901
902 """
903 n = self.dimension()
904 # Prepend the header row.
905 M = [["*"] + list(self.gens())]
906
907 # And to each subsequent row, prepend an entry that belongs to
908 # the left-side "header column."
909 M += [ [self.monomial(i)] + [ self.monomial(i)*self.monomial(j)
910 for j in range(n) ]
911 for i in range(n) ]
912
913 return table(M, header_row=True, header_column=True, frame=True)
914
915
916 def matrix_basis(self):
917 """
918 Return an (often more natural) representation of this algebras
919 basis as an ordered tuple of matrices.
920
921 Every finite-dimensional Euclidean Jordan Algebra is a, up to
922 Jordan isomorphism, a direct sum of five simple
923 algebras---four of which comprise Hermitian matrices. And the
924 last type of algebra can of course be thought of as `n`-by-`1`
925 column matrices (ambiguusly called column vectors) to avoid
926 special cases. As a result, matrices (and column vectors) are
927 a natural representation format for Euclidean Jordan algebra
928 elements.
929
930 But, when we construct an algebra from a basis of matrices,
931 those matrix representations are lost in favor of coordinate
932 vectors *with respect to* that basis. We could eventually
933 convert back if we tried hard enough, but having the original
934 representations handy is valuable enough that we simply store
935 them and return them from this method.
936
937 Why implement this for non-matrix algebras? Avoiding special
938 cases for the :class:`BilinearFormEJA` pays with simplicity in
939 its own right. But mainly, we would like to be able to assume
940 that elements of a :class:`CartesianProductEJA` can be displayed
941 nicely, without having to have special classes for direct sums
942 one of whose components was a matrix algebra.
943
944 SETUP::
945
946 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
947 ....: RealSymmetricEJA)
948
949 EXAMPLES::
950
951 sage: J = RealSymmetricEJA(2)
952 sage: J.basis()
953 Finite family {0: b0, 1: b1, 2: b2}
954 sage: J.matrix_basis()
955 (
956 [1 0] [ 0 0.7071067811865475?] [0 0]
957 [0 0], [0.7071067811865475? 0], [0 1]
958 )
959
960 ::
961
962 sage: J = JordanSpinEJA(2)
963 sage: J.basis()
964 Finite family {0: b0, 1: b1}
965 sage: J.matrix_basis()
966 (
967 [1] [0]
968 [0], [1]
969 )
970 """
971 return self._matrix_basis
972
973
974 def matrix_space(self):
975 """
976 Return the matrix space in which this algebra's elements live, if
977 we think of them as matrices (including column vectors of the
978 appropriate size).
979
980 "By default" this will be an `n`-by-`1` column-matrix space,
981 except when the algebra is trivial. There it's `n`-by-`n`
982 (where `n` is zero), to ensure that two elements of the matrix
983 space (empty matrices) can be multiplied. For algebras of
984 matrices, this returns the space in which their
985 real embeddings live.
986
987 SETUP::
988
989 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
990 ....: JordanSpinEJA,
991 ....: QuaternionHermitianEJA,
992 ....: TrivialEJA)
993
994 EXAMPLES:
995
996 By default, the matrix representation is just a column-matrix
997 equivalent to the vector representation::
998
999 sage: J = JordanSpinEJA(3)
1000 sage: J.matrix_space()
1001 Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1002 Real Field
1003
1004 The matrix representation in the trivial algebra is
1005 zero-by-zero instead of the usual `n`-by-one::
1006
1007 sage: J = TrivialEJA()
1008 sage: J.matrix_space()
1009 Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1010 Real Field
1011
1012 The matrix space for complex/quaternion Hermitian matrix EJA
1013 is the space in which their real-embeddings live, not the
1014 original complex/quaternion matrix space::
1015
1016 sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1017 sage: J.matrix_space()
1018 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1019 sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1020 sage: J.matrix_space()
1021 Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1022
1023 """
1024 if self.is_trivial():
1025 return MatrixSpace(self.base_ring(), 0)
1026 else:
1027 return self.matrix_basis()[0].parent()
1028
1029
1030 @cached_method
1031 def one(self):
1032 """
1033 Return the unit element of this algebra.
1034
1035 SETUP::
1036
1037 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1038 ....: random_eja)
1039
1040 EXAMPLES:
1041
1042 We can compute unit element in the Hadamard EJA::
1043
1044 sage: J = HadamardEJA(5)
1045 sage: J.one()
1046 b0 + b1 + b2 + b3 + b4
1047
1048 The unit element in the Hadamard EJA is inherited in the
1049 subalgebras generated by its elements::
1050
1051 sage: J = HadamardEJA(5)
1052 sage: J.one()
1053 b0 + b1 + b2 + b3 + b4
1054 sage: x = sum(J.gens())
1055 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1056 sage: A.one()
1057 c0
1058 sage: A.one().superalgebra_element()
1059 b0 + b1 + b2 + b3 + b4
1060
1061 TESTS:
1062
1063 The identity element acts like the identity, regardless of
1064 whether or not we orthonormalize::
1065
1066 sage: set_random_seed()
1067 sage: J = random_eja()
1068 sage: x = J.random_element()
1069 sage: J.one()*x == x and x*J.one() == x
1070 True
1071 sage: A = x.subalgebra_generated_by()
1072 sage: y = A.random_element()
1073 sage: A.one()*y == y and y*A.one() == y
1074 True
1075
1076 ::
1077
1078 sage: set_random_seed()
1079 sage: J = random_eja(field=QQ, orthonormalize=False)
1080 sage: x = J.random_element()
1081 sage: J.one()*x == x and x*J.one() == x
1082 True
1083 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1084 sage: y = A.random_element()
1085 sage: A.one()*y == y and y*A.one() == y
1086 True
1087
1088 The matrix of the unit element's operator is the identity,
1089 regardless of the base field and whether or not we
1090 orthonormalize::
1091
1092 sage: set_random_seed()
1093 sage: J = random_eja()
1094 sage: actual = J.one().operator().matrix()
1095 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1096 sage: actual == expected
1097 True
1098 sage: x = J.random_element()
1099 sage: A = x.subalgebra_generated_by()
1100 sage: actual = A.one().operator().matrix()
1101 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1102 sage: actual == expected
1103 True
1104
1105 ::
1106
1107 sage: set_random_seed()
1108 sage: J = random_eja(field=QQ, orthonormalize=False)
1109 sage: actual = J.one().operator().matrix()
1110 sage: expected = matrix.identity(J.base_ring(), J.dimension())
1111 sage: actual == expected
1112 True
1113 sage: x = J.random_element()
1114 sage: A = x.subalgebra_generated_by(orthonormalize=False)
1115 sage: actual = A.one().operator().matrix()
1116 sage: expected = matrix.identity(A.base_ring(), A.dimension())
1117 sage: actual == expected
1118 True
1119
1120 Ensure that the cached unit element (often precomputed by
1121 hand) agrees with the computed one::
1122
1123 sage: set_random_seed()
1124 sage: J = random_eja()
1125 sage: cached = J.one()
1126 sage: J.one.clear_cache()
1127 sage: J.one() == cached
1128 True
1129
1130 ::
1131
1132 sage: set_random_seed()
1133 sage: J = random_eja(field=QQ, orthonormalize=False)
1134 sage: cached = J.one()
1135 sage: J.one.clear_cache()
1136 sage: J.one() == cached
1137 True
1138
1139 """
1140 # We can brute-force compute the matrices of the operators
1141 # that correspond to the basis elements of this algebra.
1142 # If some linear combination of those basis elements is the
1143 # algebra identity, then the same linear combination of
1144 # their matrices has to be the identity matrix.
1145 #
1146 # Of course, matrices aren't vectors in sage, so we have to
1147 # appeal to the "long vectors" isometry.
1148 oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
1149
1150 # Now we use basic linear algebra to find the coefficients,
1151 # of the matrices-as-vectors-linear-combination, which should
1152 # work for the original algebra basis too.
1153 A = matrix(self.base_ring(), oper_vecs)
1154
1155 # We used the isometry on the left-hand side already, but we
1156 # still need to do it for the right-hand side. Recall that we
1157 # wanted something that summed to the identity matrix.
1158 b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
1159
1160 # Now if there's an identity element in the algebra, this
1161 # should work. We solve on the left to avoid having to
1162 # transpose the matrix "A".
1163 return self.from_vector(A.solve_left(b))
1164
1165
1166 def peirce_decomposition(self, c):
1167 """
1168 The Peirce decomposition of this algebra relative to the
1169 idempotent ``c``.
1170
1171 In the future, this can be extended to a complete system of
1172 orthogonal idempotents.
1173
1174 INPUT:
1175
1176 - ``c`` -- an idempotent of this algebra.
1177
1178 OUTPUT:
1179
1180 A triple (J0, J5, J1) containing two subalgebras and one subspace
1181 of this algebra,
1182
1183 - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1184 corresponding to the eigenvalue zero.
1185
1186 - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1187 corresponding to the eigenvalue one-half.
1188
1189 - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1190 corresponding to the eigenvalue one.
1191
1192 These are the only possible eigenspaces for that operator, and this
1193 algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1194 orthogonal, and are subalgebras of this algebra with the appropriate
1195 restrictions.
1196
1197 SETUP::
1198
1199 sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1200
1201 EXAMPLES:
1202
1203 The canonical example comes from the symmetric matrices, which
1204 decompose into diagonal and off-diagonal parts::
1205
1206 sage: J = RealSymmetricEJA(3)
1207 sage: C = matrix(QQ, [ [1,0,0],
1208 ....: [0,1,0],
1209 ....: [0,0,0] ])
1210 sage: c = J(C)
1211 sage: J0,J5,J1 = J.peirce_decomposition(c)
1212 sage: J0
1213 Euclidean Jordan algebra of dimension 1...
1214 sage: J5
1215 Vector space of degree 6 and dimension 2...
1216 sage: J1
1217 Euclidean Jordan algebra of dimension 3...
1218 sage: J0.one().to_matrix()
1219 [0 0 0]
1220 [0 0 0]
1221 [0 0 1]
1222 sage: orig_df = AA.options.display_format
1223 sage: AA.options.display_format = 'radical'
1224 sage: J.from_vector(J5.basis()[0]).to_matrix()
1225 [ 0 0 1/2*sqrt(2)]
1226 [ 0 0 0]
1227 [1/2*sqrt(2) 0 0]
1228 sage: J.from_vector(J5.basis()[1]).to_matrix()
1229 [ 0 0 0]
1230 [ 0 0 1/2*sqrt(2)]
1231 [ 0 1/2*sqrt(2) 0]
1232 sage: AA.options.display_format = orig_df
1233 sage: J1.one().to_matrix()
1234 [1 0 0]
1235 [0 1 0]
1236 [0 0 0]
1237
1238 TESTS:
1239
1240 Every algebra decomposes trivially with respect to its identity
1241 element::
1242
1243 sage: set_random_seed()
1244 sage: J = random_eja()
1245 sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1246 sage: J0.dimension() == 0 and J5.dimension() == 0
1247 True
1248 sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1249 True
1250
1251 The decomposition is into eigenspaces, and its components are
1252 therefore necessarily orthogonal. Moreover, the identity
1253 elements in the two subalgebras are the projections onto their
1254 respective subspaces of the superalgebra's identity element::
1255
1256 sage: set_random_seed()
1257 sage: J = random_eja()
1258 sage: x = J.random_element()
1259 sage: if not J.is_trivial():
1260 ....: while x.is_nilpotent():
1261 ....: x = J.random_element()
1262 sage: c = x.subalgebra_idempotent()
1263 sage: J0,J5,J1 = J.peirce_decomposition(c)
1264 sage: ipsum = 0
1265 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1266 ....: w = w.superalgebra_element()
1267 ....: y = J.from_vector(y)
1268 ....: z = z.superalgebra_element()
1269 ....: ipsum += w.inner_product(y).abs()
1270 ....: ipsum += w.inner_product(z).abs()
1271 ....: ipsum += y.inner_product(z).abs()
1272 sage: ipsum
1273 0
1274 sage: J1(c) == J1.one()
1275 True
1276 sage: J0(J.one() - c) == J0.one()
1277 True
1278
1279 """
1280 if not c.is_idempotent():
1281 raise ValueError("element is not idempotent: %s" % c)
1282
1283 # Default these to what they should be if they turn out to be
1284 # trivial, because eigenspaces_left() won't return eigenvalues
1285 # corresponding to trivial spaces (e.g. it returns only the
1286 # eigenspace corresponding to lambda=1 if you take the
1287 # decomposition relative to the identity element).
1288 trivial = self.subalgebra(())
1289 J0 = trivial # eigenvalue zero
1290 J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
1291 J1 = trivial # eigenvalue one
1292
1293 for (eigval, eigspace) in c.operator().matrix().right_eigenspaces():
1294 if eigval == ~(self.base_ring()(2)):
1295 J5 = eigspace
1296 else:
1297 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
1298 subalg = self.subalgebra(gens, check_axioms=False)
1299 if eigval == 0:
1300 J0 = subalg
1301 elif eigval == 1:
1302 J1 = subalg
1303 else:
1304 raise ValueError("unexpected eigenvalue: %s" % eigval)
1305
1306 return (J0, J5, J1)
1307
1308
1309 def random_element(self, thorough=False):
1310 r"""
1311 Return a random element of this algebra.
1312
1313 Our algebra superclass method only returns a linear
1314 combination of at most two basis elements. We instead
1315 want the vector space "random element" method that
1316 returns a more diverse selection.
1317
1318 INPUT:
1319
1320 - ``thorough`` -- (boolean; default False) whether or not we
1321 should generate irrational coefficients for the random
1322 element when our base ring is irrational; this slows the
1323 algebra operations to a crawl, but any truly random method
1324 should include them
1325
1326 """
1327 # For a general base ring... maybe we can trust this to do the
1328 # right thing? Unlikely, but.
1329 V = self.vector_space()
1330 v = V.random_element()
1331
1332 if self.base_ring() is AA:
1333 # The "random element" method of the algebraic reals is
1334 # stupid at the moment, and only returns integers between
1335 # -2 and 2, inclusive:
1336 #
1337 # https://trac.sagemath.org/ticket/30875
1338 #
1339 # Instead, we implement our own "random vector" method,
1340 # and then coerce that into the algebra. We use the vector
1341 # space degree here instead of the dimension because a
1342 # subalgebra could (for example) be spanned by only two
1343 # vectors, each with five coordinates. We need to
1344 # generate all five coordinates.
1345 if thorough:
1346 v *= QQbar.random_element().real()
1347 else:
1348 v *= QQ.random_element()
1349
1350 return self.from_vector(V.coordinate_vector(v))
1351
1352 def random_elements(self, count, thorough=False):
1353 """
1354 Return ``count`` random elements as a tuple.
1355
1356 INPUT:
1357
1358 - ``thorough`` -- (boolean; default False) whether or not we
1359 should generate irrational coefficients for the random
1360 elements when our base ring is irrational; this slows the
1361 algebra operations to a crawl, but any truly random method
1362 should include them
1363
1364 SETUP::
1365
1366 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1367
1368 EXAMPLES::
1369
1370 sage: J = JordanSpinEJA(3)
1371 sage: x,y,z = J.random_elements(3)
1372 sage: all( [ x in J, y in J, z in J ])
1373 True
1374 sage: len( J.random_elements(10) ) == 10
1375 True
1376
1377 """
1378 return tuple( self.random_element(thorough)
1379 for idx in range(count) )
1380
1381
1382 @cached_method
1383 def _charpoly_coefficients(self):
1384 r"""
1385 The `r` polynomial coefficients of the "characteristic polynomial
1386 of" function.
1387
1388 SETUP::
1389
1390 sage: from mjo.eja.eja_algebra import random_eja
1391
1392 TESTS:
1393
1394 The theory shows that these are all homogeneous polynomials of
1395 a known degree::
1396
1397 sage: set_random_seed()
1398 sage: J = random_eja()
1399 sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1400 True
1401
1402 """
1403 n = self.dimension()
1404 R = self.coordinate_polynomial_ring()
1405 vars = R.gens()
1406 F = R.fraction_field()
1407
1408 def L_x_i_j(i,j):
1409 # From a result in my book, these are the entries of the
1410 # basis representation of L_x.
1411 return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
1412 for k in range(n) )
1413
1414 L_x = matrix(F, n, n, L_x_i_j)
1415
1416 r = None
1417 if self.rank.is_in_cache():
1418 r = self.rank()
1419 # There's no need to pad the system with redundant
1420 # columns if we *know* they'll be redundant.
1421 n = r
1422
1423 # Compute an extra power in case the rank is equal to
1424 # the dimension (otherwise, we would stop at x^(r-1)).
1425 x_powers = [ (L_x**k)*self.one().to_vector()
1426 for k in range(n+1) ]
1427 A = matrix.column(F, x_powers[:n])
1428 AE = A.extended_echelon_form()
1429 E = AE[:,n:]
1430 A_rref = AE[:,:n]
1431 if r is None:
1432 r = A_rref.rank()
1433 b = x_powers[r]
1434
1435 # The theory says that only the first "r" coefficients are
1436 # nonzero, and they actually live in the original polynomial
1437 # ring and not the fraction field. We negate them because in
1438 # the actual characteristic polynomial, they get moved to the
1439 # other side where x^r lives. We don't bother to trim A_rref
1440 # down to a square matrix and solve the resulting system,
1441 # because the upper-left r-by-r portion of A_rref is
1442 # guaranteed to be the identity matrix, so e.g.
1443 #
1444 # A_rref.solve_right(Y)
1445 #
1446 # would just be returning Y.
1447 return (-E*b)[:r].change_ring(R)
1448
1449 @cached_method
1450 def rank(self):
1451 r"""
1452 Return the rank of this EJA.
1453
1454 This is a cached method because we know the rank a priori for
1455 all of the algebras we can construct. Thus we can avoid the
1456 expensive ``_charpoly_coefficients()`` call unless we truly
1457 need to compute the whole characteristic polynomial.
1458
1459 SETUP::
1460
1461 sage: from mjo.eja.eja_algebra import (HadamardEJA,
1462 ....: JordanSpinEJA,
1463 ....: RealSymmetricEJA,
1464 ....: ComplexHermitianEJA,
1465 ....: QuaternionHermitianEJA,
1466 ....: random_eja)
1467
1468 EXAMPLES:
1469
1470 The rank of the Jordan spin algebra is always two::
1471
1472 sage: JordanSpinEJA(2).rank()
1473 2
1474 sage: JordanSpinEJA(3).rank()
1475 2
1476 sage: JordanSpinEJA(4).rank()
1477 2
1478
1479 The rank of the `n`-by-`n` Hermitian real, complex, or
1480 quaternion matrices is `n`::
1481
1482 sage: RealSymmetricEJA(4).rank()
1483 4
1484 sage: ComplexHermitianEJA(3).rank()
1485 3
1486 sage: QuaternionHermitianEJA(2).rank()
1487 2
1488
1489 TESTS:
1490
1491 Ensure that every EJA that we know how to construct has a
1492 positive integer rank, unless the algebra is trivial in
1493 which case its rank will be zero::
1494
1495 sage: set_random_seed()
1496 sage: J = random_eja()
1497 sage: r = J.rank()
1498 sage: r in ZZ
1499 True
1500 sage: r > 0 or (r == 0 and J.is_trivial())
1501 True
1502
1503 Ensure that computing the rank actually works, since the ranks
1504 of all simple algebras are known and will be cached by default::
1505
1506 sage: set_random_seed() # long time
1507 sage: J = random_eja() # long time
1508 sage: cached = J.rank() # long time
1509 sage: J.rank.clear_cache() # long time
1510 sage: J.rank() == cached # long time
1511 True
1512
1513 """
1514 return len(self._charpoly_coefficients())
1515
1516
1517 def subalgebra(self, basis, **kwargs):
1518 r"""
1519 Create a subalgebra of this algebra from the given basis.
1520 """
1521 from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
1522 return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
1523
1524
1525 def vector_space(self):
1526 """
1527 Return the vector space that underlies this algebra.
1528
1529 SETUP::
1530
1531 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1532
1533 EXAMPLES::
1534
1535 sage: J = RealSymmetricEJA(2)
1536 sage: J.vector_space()
1537 Vector space of dimension 3 over...
1538
1539 """
1540 return self.zero().to_vector().parent().ambient_vector_space()
1541
1542
1543
1544 class RationalBasisEJA(FiniteDimensionalEJA):
1545 r"""
1546 New class for algebras whose supplied basis elements have all rational entries.
1547
1548 SETUP::
1549
1550 sage: from mjo.eja.eja_algebra import BilinearFormEJA
1551
1552 EXAMPLES:
1553
1554 The supplied basis is orthonormalized by default::
1555
1556 sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1557 sage: J = BilinearFormEJA(B)
1558 sage: J.matrix_basis()
1559 (
1560 [1] [ 0] [ 0]
1561 [0] [1/5] [32/5]
1562 [0], [ 0], [ 5]
1563 )
1564
1565 """
1566 def __init__(self,
1567 basis,
1568 jordan_product,
1569 inner_product,
1570 field=AA,
1571 check_field=True,
1572 **kwargs):
1573
1574 if check_field:
1575 # Abuse the check_field parameter to check that the entries of
1576 # out basis (in ambient coordinates) are in the field QQ.
1577 if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
1578 raise TypeError("basis not rational")
1579
1580 super().__init__(basis,
1581 jordan_product,
1582 inner_product,
1583 field=field,
1584 check_field=check_field,
1585 **kwargs)
1586
1587 self._rational_algebra = None
1588 if field is not QQ:
1589 # There's no point in constructing the extra algebra if this
1590 # one is already rational.
1591 #
1592 # Note: the same Jordan and inner-products work here,
1593 # because they are necessarily defined with respect to
1594 # ambient coordinates and not any particular basis.
1595 self._rational_algebra = FiniteDimensionalEJA(
1596 basis,
1597 jordan_product,
1598 inner_product,
1599 field=QQ,
1600 associative=self.is_associative(),
1601 orthonormalize=False,
1602 check_field=False,
1603 check_axioms=False)
1604
1605 @cached_method
1606 def _charpoly_coefficients(self):
1607 r"""
1608 SETUP::
1609
1610 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1611 ....: JordanSpinEJA)
1612
1613 EXAMPLES:
1614
1615 The base ring of the resulting polynomial coefficients is what
1616 it should be, and not the rationals (unless the algebra was
1617 already over the rationals)::
1618
1619 sage: J = JordanSpinEJA(3)
1620 sage: J._charpoly_coefficients()
1621 (X1^2 - X2^2 - X3^2, -2*X1)
1622 sage: a0 = J._charpoly_coefficients()[0]
1623 sage: J.base_ring()
1624 Algebraic Real Field
1625 sage: a0.base_ring()
1626 Algebraic Real Field
1627
1628 """
1629 if self._rational_algebra is None:
1630 # There's no need to construct *another* algebra over the
1631 # rationals if this one is already over the
1632 # rationals. Likewise, if we never orthonormalized our
1633 # basis, we might as well just use the given one.
1634 return super()._charpoly_coefficients()
1635
1636 # Do the computation over the rationals. The answer will be
1637 # the same, because all we've done is a change of basis.
1638 # Then, change back from QQ to our real base ring
1639 a = ( a_i.change_ring(self.base_ring())
1640 for a_i in self._rational_algebra._charpoly_coefficients() )
1641
1642 if self._deortho_matrix is None:
1643 # This can happen if our base ring was, say, AA and we
1644 # chose not to (or didn't need to) orthonormalize. It's
1645 # still faster to do the computations over QQ even if
1646 # the numbers in the boxes stay the same.
1647 return tuple(a)
1648
1649 # Otherwise, convert the coordinate variables back to the
1650 # deorthonormalized ones.
1651 R = self.coordinate_polynomial_ring()
1652 from sage.modules.free_module_element import vector
1653 X = vector(R, R.gens())
1654 BX = self._deortho_matrix*X
1655
1656 subs_dict = { X[i]: BX[i] for i in range(len(X)) }
1657 return tuple( a_i.subs(subs_dict) for a_i in a )
1658
1659 class ConcreteEJA(RationalBasisEJA):
1660 r"""
1661 A class for the Euclidean Jordan algebras that we know by name.
1662
1663 These are the Jordan algebras whose basis, multiplication table,
1664 rank, and so on are known a priori. More to the point, they are
1665 the Euclidean Jordan algebras for which we are able to conjure up
1666 a "random instance."
1667
1668 SETUP::
1669
1670 sage: from mjo.eja.eja_algebra import ConcreteEJA
1671
1672 TESTS:
1673
1674 Our basis is normalized with respect to the algebra's inner
1675 product, unless we specify otherwise::
1676
1677 sage: set_random_seed()
1678 sage: J = ConcreteEJA.random_instance()
1679 sage: all( b.norm() == 1 for b in J.gens() )
1680 True
1681
1682 Since our basis is orthonormal with respect to the algebra's inner
1683 product, and since we know that this algebra is an EJA, any
1684 left-multiplication operator's matrix will be symmetric because
1685 natural->EJA basis representation is an isometry and within the
1686 EJA the operator is self-adjoint by the Jordan axiom::
1687
1688 sage: set_random_seed()
1689 sage: J = ConcreteEJA.random_instance()
1690 sage: x = J.random_element()
1691 sage: x.operator().is_self_adjoint()
1692 True
1693 """
1694
1695 @staticmethod
1696 def _max_random_instance_size():
1697 """
1698 Return an integer "size" that is an upper bound on the size of
1699 this algebra when it is used in a random test
1700 case. Unfortunately, the term "size" is ambiguous -- when
1701 dealing with `R^n` under either the Hadamard or Jordan spin
1702 product, the "size" refers to the dimension `n`. When dealing
1703 with a matrix algebra (real symmetric or complex/quaternion
1704 Hermitian), it refers to the size of the matrix, which is far
1705 less than the dimension of the underlying vector space.
1706
1707 This method must be implemented in each subclass.
1708 """
1709 raise NotImplementedError
1710
1711 @classmethod
1712 def random_instance(cls, *args, **kwargs):
1713 """
1714 Return a random instance of this type of algebra.
1715
1716 This method should be implemented in each subclass.
1717 """
1718 from sage.misc.prandom import choice
1719 eja_class = choice(cls.__subclasses__())
1720
1721 # These all bubble up to the RationalBasisEJA superclass
1722 # constructor, so any (kw)args valid there are also valid
1723 # here.
1724 return eja_class.random_instance(*args, **kwargs)
1725
1726
1727 class MatrixEJA:
1728 @staticmethod
1729 def jordan_product(X,Y):
1730 return (X*Y + Y*X)/2
1731
1732 @staticmethod
1733 def trace_inner_product(X,Y):
1734 r"""
1735 A trace inner-product for matrices that aren't embedded in the
1736 reals.
1737 """
1738 # We take the norm (absolute value) because Octonions() isn't
1739 # smart enough yet to coerce its one() into the base field.
1740 return (X*Y).trace().real().abs()
1741
1742 class RealEmbeddedMatrixEJA(MatrixEJA):
1743 @staticmethod
1744 def dimension_over_reals():
1745 r"""
1746 The dimension of this matrix's base ring over the reals.
1747
1748 The reals are dimension one over themselves, obviously; that's
1749 just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1750 have dimension two. Finally, the quaternions have dimension
1751 four over the reals.
1752
1753 This is used to determine the size of the matrix returned from
1754 :meth:`real_embed`, among other things.
1755 """
1756 raise NotImplementedError
1757
1758 @classmethod
1759 def real_embed(cls,M):
1760 """
1761 Embed the matrix ``M`` into a space of real matrices.
1762
1763 The matrix ``M`` can have entries in any field at the moment:
1764 the real numbers, complex numbers, or quaternions. And although
1765 they are not a field, we can probably support octonions at some
1766 point, too. This function returns a real matrix that "acts like"
1767 the original with respect to matrix multiplication; i.e.
1768
1769 real_embed(M*N) = real_embed(M)*real_embed(N)
1770
1771 """
1772 if M.ncols() != M.nrows():
1773 raise ValueError("the matrix 'M' must be square")
1774 return M
1775
1776
1777 @classmethod
1778 def real_unembed(cls,M):
1779 """
1780 The inverse of :meth:`real_embed`.
1781 """
1782 if M.ncols() != M.nrows():
1783 raise ValueError("the matrix 'M' must be square")
1784 if not ZZ(M.nrows()).mod(cls.dimension_over_reals()).is_zero():
1785 raise ValueError("the matrix 'M' must be a real embedding")
1786 return M
1787
1788
1789 @classmethod
1790 def trace_inner_product(cls,X,Y):
1791 r"""
1792 Compute the trace inner-product of two real-embeddings.
1793
1794 SETUP::
1795
1796 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1797 ....: QuaternionHermitianEJA)
1798
1799 EXAMPLES::
1800
1801 sage: set_random_seed()
1802 sage: J = ComplexHermitianEJA.random_instance()
1803 sage: x,y = J.random_elements(2)
1804 sage: Xe = x.to_matrix()
1805 sage: Ye = y.to_matrix()
1806 sage: X = J.real_unembed(Xe)
1807 sage: Y = J.real_unembed(Ye)
1808 sage: expected = (X*Y).trace().real()
1809 sage: actual = J.trace_inner_product(Xe,Ye)
1810 sage: actual == expected
1811 True
1812
1813 ::
1814
1815 sage: set_random_seed()
1816 sage: J = QuaternionHermitianEJA.random_instance()
1817 sage: x,y = J.random_elements(2)
1818 sage: Xe = x.to_matrix()
1819 sage: Ye = y.to_matrix()
1820 sage: X = J.real_unembed(Xe)
1821 sage: Y = J.real_unembed(Ye)
1822 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1823 sage: actual = J.trace_inner_product(Xe,Ye)
1824 sage: actual == expected
1825 True
1826
1827 """
1828 # This does in fact compute the real part of the trace.
1829 # If we compute the trace of e.g. a complex matrix M,
1830 # then we do so by adding up its diagonal entries --
1831 # call them z_1 through z_n. The real embedding of z_1
1832 # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
1833 # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
1834 return (X*Y).trace()/cls.dimension_over_reals()
1835
1836 class RealSymmetricEJA(ConcreteEJA, MatrixEJA):
1837 """
1838 The rank-n simple EJA consisting of real symmetric n-by-n
1839 matrices, the usual symmetric Jordan product, and the trace inner
1840 product. It has dimension `(n^2 + n)/2` over the reals.
1841
1842 SETUP::
1843
1844 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1845
1846 EXAMPLES::
1847
1848 sage: J = RealSymmetricEJA(2)
1849 sage: b0, b1, b2 = J.gens()
1850 sage: b0*b0
1851 b0
1852 sage: b1*b1
1853 1/2*b0 + 1/2*b2
1854 sage: b2*b2
1855 b2
1856
1857 In theory, our "field" can be any subfield of the reals::
1858
1859 sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1860 Euclidean Jordan algebra of dimension 3 over Real Double Field
1861 sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1862 Euclidean Jordan algebra of dimension 3 over Real Field with
1863 53 bits of precision
1864
1865 TESTS:
1866
1867 The dimension of this algebra is `(n^2 + n) / 2`::
1868
1869 sage: set_random_seed()
1870 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1871 sage: n = ZZ.random_element(1, n_max)
1872 sage: J = RealSymmetricEJA(n)
1873 sage: J.dimension() == (n^2 + n)/2
1874 True
1875
1876 The Jordan multiplication is what we think it is::
1877
1878 sage: set_random_seed()
1879 sage: J = RealSymmetricEJA.random_instance()
1880 sage: x,y = J.random_elements(2)
1881 sage: actual = (x*y).to_matrix()
1882 sage: X = x.to_matrix()
1883 sage: Y = y.to_matrix()
1884 sage: expected = (X*Y + Y*X)/2
1885 sage: actual == expected
1886 True
1887 sage: J(expected) == x*y
1888 True
1889
1890 We can change the generator prefix::
1891
1892 sage: RealSymmetricEJA(3, prefix='q').gens()
1893 (q0, q1, q2, q3, q4, q5)
1894
1895 We can construct the (trivial) algebra of rank zero::
1896
1897 sage: RealSymmetricEJA(0)
1898 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1899
1900 """
1901 @classmethod
1902 def _denormalized_basis(cls, n, field):
1903 """
1904 Return a basis for the space of real symmetric n-by-n matrices.
1905
1906 SETUP::
1907
1908 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1909
1910 TESTS::
1911
1912 sage: set_random_seed()
1913 sage: n = ZZ.random_element(1,5)
1914 sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)
1915 sage: all( M.is_symmetric() for M in B)
1916 True
1917
1918 """
1919 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1920 # coordinates.
1921 S = []
1922 for i in range(n):
1923 for j in range(i+1):
1924 Eij = matrix(field, n, lambda k,l: k==i and l==j)
1925 if i == j:
1926 Sij = Eij
1927 else:
1928 Sij = Eij + Eij.transpose()
1929 S.append(Sij)
1930 return tuple(S)
1931
1932
1933 @staticmethod
1934 def _max_random_instance_size():
1935 return 4 # Dimension 10
1936
1937 @classmethod
1938 def random_instance(cls, **kwargs):
1939 """
1940 Return a random instance of this type of algebra.
1941 """
1942 n = ZZ.random_element(cls._max_random_instance_size() + 1)
1943 return cls(n, **kwargs)
1944
1945 def __init__(self, n, field=AA, **kwargs):
1946 # We know this is a valid EJA, but will double-check
1947 # if the user passes check_axioms=True.
1948 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
1949
1950 associative = False
1951 if n <= 1:
1952 associative = True
1953
1954 super().__init__(self._denormalized_basis(n,field),
1955 self.jordan_product,
1956 self.trace_inner_product,
1957 field=field,
1958 associative=associative,
1959 **kwargs)
1960
1961 # TODO: this could be factored out somehow, but is left here
1962 # because the MatrixEJA is not presently a subclass of the
1963 # FDEJA class that defines rank() and one().
1964 self.rank.set_cache(n)
1965 idV = self.matrix_space().one()
1966 self.one.set_cache(self(idV))
1967
1968
1969
1970 class ComplexMatrixEJA(RealEmbeddedMatrixEJA):
1971 # A manual dictionary-cache for the complex_extension() method,
1972 # since apparently @classmethods can't also be @cached_methods.
1973 _complex_extension = {}
1974
1975 @classmethod
1976 def complex_extension(cls,field):
1977 r"""
1978 The complex field that we embed/unembed, as an extension
1979 of the given ``field``.
1980 """
1981 if field in cls._complex_extension:
1982 return cls._complex_extension[field]
1983
1984 # Sage doesn't know how to adjoin the complex "i" (the root of
1985 # x^2 + 1) to a field in a general way. Here, we just enumerate
1986 # all of the cases that I have cared to support so far.
1987 if field is AA:
1988 # Sage doesn't know how to embed AA into QQbar, i.e. how
1989 # to adjoin sqrt(-1) to AA.
1990 F = QQbar
1991 elif not field.is_exact():
1992 # RDF or RR
1993 F = field.complex_field()
1994 else:
1995 # Works for QQ and... maybe some other fields.
1996 R = PolynomialRing(field, 'z')
1997 z = R.gen()
1998 F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
1999
2000 cls._complex_extension[field] = F
2001 return F
2002
2003 @staticmethod
2004 def dimension_over_reals():
2005 return 2
2006
2007 @classmethod
2008 def real_embed(cls,M):
2009 """
2010 Embed the n-by-n complex matrix ``M`` into the space of real
2011 matrices of size 2n-by-2n via the map the sends each entry `z = a +
2012 bi` to the block matrix ``[[a,b],[-b,a]]``.
2013
2014 SETUP::
2015
2016 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2017
2018 EXAMPLES::
2019
2020 sage: F = QuadraticField(-1, 'I')
2021 sage: x1 = F(4 - 2*i)
2022 sage: x2 = F(1 + 2*i)
2023 sage: x3 = F(-i)
2024 sage: x4 = F(6)
2025 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2026 sage: ComplexMatrixEJA.real_embed(M)
2027 [ 4 -2| 1 2]
2028 [ 2 4|-2 1]
2029 [-----+-----]
2030 [ 0 -1| 6 0]
2031 [ 1 0| 0 6]
2032
2033 TESTS:
2034
2035 Embedding is a homomorphism (isomorphism, in fact)::
2036
2037 sage: set_random_seed()
2038 sage: n = ZZ.random_element(3)
2039 sage: F = QuadraticField(-1, 'I')
2040 sage: X = random_matrix(F, n)
2041 sage: Y = random_matrix(F, n)
2042 sage: Xe = ComplexMatrixEJA.real_embed(X)
2043 sage: Ye = ComplexMatrixEJA.real_embed(Y)
2044 sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
2045 sage: Xe*Ye == XYe
2046 True
2047
2048 """
2049 super().real_embed(M)
2050 n = M.nrows()
2051
2052 # We don't need any adjoined elements...
2053 field = M.base_ring().base_ring()
2054
2055 blocks = []
2056 for z in M.list():
2057 a = z.real()
2058 b = z.imag()
2059 blocks.append(matrix(field, 2, [ [ a, b],
2060 [-b, a] ]))
2061
2062 return matrix.block(field, n, blocks)
2063
2064
2065 @classmethod
2066 def real_unembed(cls,M):
2067 """
2068 The inverse of _embed_complex_matrix().
2069
2070 SETUP::
2071
2072 sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2073
2074 EXAMPLES::
2075
2076 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
2077 ....: [-2, 1, -4, 3],
2078 ....: [ 9, 10, 11, 12],
2079 ....: [-10, 9, -12, 11] ])
2080 sage: ComplexMatrixEJA.real_unembed(A)
2081 [ 2*I + 1 4*I + 3]
2082 [ 10*I + 9 12*I + 11]
2083
2084 TESTS:
2085
2086 Unembedding is the inverse of embedding::
2087
2088 sage: set_random_seed()
2089 sage: F = QuadraticField(-1, 'I')
2090 sage: M = random_matrix(F, 3)
2091 sage: Me = ComplexMatrixEJA.real_embed(M)
2092 sage: ComplexMatrixEJA.real_unembed(Me) == M
2093 True
2094
2095 """
2096 super().real_unembed(M)
2097 n = ZZ(M.nrows())
2098 d = cls.dimension_over_reals()
2099 F = cls.complex_extension(M.base_ring())
2100 i = F.gen()
2101
2102 # Go top-left to bottom-right (reading order), converting every
2103 # 2-by-2 block we see to a single complex element.
2104 elements = []
2105 for k in range(n/d):
2106 for j in range(n/d):
2107 submat = M[d*k:d*k+d,d*j:d*j+d]
2108 if submat[0,0] != submat[1,1]:
2109 raise ValueError('bad on-diagonal submatrix')
2110 if submat[0,1] != -submat[1,0]:
2111 raise ValueError('bad off-diagonal submatrix')
2112 z = submat[0,0] + submat[0,1]*i
2113 elements.append(z)
2114
2115 return matrix(F, n/d, elements)
2116
2117
2118 class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
2119 """
2120 The rank-n simple EJA consisting of complex Hermitian n-by-n
2121 matrices over the real numbers, the usual symmetric Jordan product,
2122 and the real-part-of-trace inner product. It has dimension `n^2` over
2123 the reals.
2124
2125 SETUP::
2126
2127 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2128
2129 EXAMPLES:
2130
2131 In theory, our "field" can be any subfield of the reals::
2132
2133 sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2134 Euclidean Jordan algebra of dimension 4 over Real Double Field
2135 sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2136 Euclidean Jordan algebra of dimension 4 over Real Field with
2137 53 bits of precision
2138
2139 TESTS:
2140
2141 The dimension of this algebra is `n^2`::
2142
2143 sage: set_random_seed()
2144 sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2145 sage: n = ZZ.random_element(1, n_max)
2146 sage: J = ComplexHermitianEJA(n)
2147 sage: J.dimension() == n^2
2148 True
2149
2150 The Jordan multiplication is what we think it is::
2151
2152 sage: set_random_seed()
2153 sage: J = ComplexHermitianEJA.random_instance()
2154 sage: x,y = J.random_elements(2)
2155 sage: actual = (x*y).to_matrix()
2156 sage: X = x.to_matrix()
2157 sage: Y = y.to_matrix()
2158 sage: expected = (X*Y + Y*X)/2
2159 sage: actual == expected
2160 True
2161 sage: J(expected) == x*y
2162 True
2163
2164 We can change the generator prefix::
2165
2166 sage: ComplexHermitianEJA(2, prefix='z').gens()
2167 (z0, z1, z2, z3)
2168
2169 We can construct the (trivial) algebra of rank zero::
2170
2171 sage: ComplexHermitianEJA(0)
2172 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2173
2174 """
2175
2176 @classmethod
2177 def _denormalized_basis(cls, n, field):
2178 """
2179 Returns a basis for the space of complex Hermitian n-by-n matrices.
2180
2181 Why do we embed these? Basically, because all of numerical linear
2182 algebra assumes that you're working with vectors consisting of `n`
2183 entries from a field and scalars from the same field. There's no way
2184 to tell SageMath that (for example) the vectors contain complex
2185 numbers, while the scalar field is real.
2186
2187 SETUP::
2188
2189 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2190
2191 TESTS::
2192
2193 sage: set_random_seed()
2194 sage: n = ZZ.random_element(1,5)
2195 sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ)
2196 sage: all( M.is_symmetric() for M in B)
2197 True
2198
2199 """
2200 R = PolynomialRing(ZZ, 'z')
2201 z = R.gen()
2202 F = ZZ.extension(z**2 + 1, 'I')
2203 I = F.gen(1)
2204
2205 # This is like the symmetric case, but we need to be careful:
2206 #
2207 # * We want conjugate-symmetry, not just symmetry.
2208 # * The diagonal will (as a result) be real.
2209 #
2210 S = []
2211 Eij = matrix.zero(F,n)
2212 for i in range(n):
2213 for j in range(i+1):
2214 # "build" E_ij
2215 Eij[i,j] = 1
2216 if i == j:
2217 Sij = cls.real_embed(Eij)
2218 S.append(Sij)
2219 else:
2220 # The second one has a minus because it's conjugated.
2221 Eij[j,i] = 1 # Eij = Eij + Eij.transpose()
2222 Sij_real = cls.real_embed(Eij)
2223 S.append(Sij_real)
2224 # Eij = I*Eij - I*Eij.transpose()
2225 Eij[i,j] = I
2226 Eij[j,i] = -I
2227 Sij_imag = cls.real_embed(Eij)
2228 S.append(Sij_imag)
2229 Eij[j,i] = 0
2230 # "erase" E_ij
2231 Eij[i,j] = 0
2232
2233 # Since we embedded the entries, we can drop back to the
2234 # desired real "field" instead of the extension "F".
2235 return tuple( s.change_ring(field) for s in S )
2236
2237
2238 def __init__(self, n, field=AA, **kwargs):
2239 # We know this is a valid EJA, but will double-check
2240 # if the user passes check_axioms=True.
2241 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2242
2243 associative = False
2244 if n <= 1:
2245 associative = True
2246
2247 super().__init__(self._denormalized_basis(n,field),
2248 self.jordan_product,
2249 self.trace_inner_product,
2250 field=field,
2251 associative=associative,
2252 **kwargs)
2253 # TODO: this could be factored out somehow, but is left here
2254 # because the MatrixEJA is not presently a subclass of the
2255 # FDEJA class that defines rank() and one().
2256 self.rank.set_cache(n)
2257 idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
2258 self.one.set_cache(self(idV))
2259
2260 @staticmethod
2261 def _max_random_instance_size():
2262 return 3 # Dimension 9
2263
2264 @classmethod
2265 def random_instance(cls, **kwargs):
2266 """
2267 Return a random instance of this type of algebra.
2268 """
2269 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2270 return cls(n, **kwargs)
2271
2272 class QuaternionMatrixEJA(RealEmbeddedMatrixEJA):
2273
2274 # A manual dictionary-cache for the quaternion_extension() method,
2275 # since apparently @classmethods can't also be @cached_methods.
2276 _quaternion_extension = {}
2277
2278 @classmethod
2279 def quaternion_extension(cls,field):
2280 r"""
2281 The quaternion field that we embed/unembed, as an extension
2282 of the given ``field``.
2283 """
2284 if field in cls._quaternion_extension:
2285 return cls._quaternion_extension[field]
2286
2287 Q = QuaternionAlgebra(field,-1,-1)
2288
2289 cls._quaternion_extension[field] = Q
2290 return Q
2291
2292 @staticmethod
2293 def dimension_over_reals():
2294 return 4
2295
2296 @classmethod
2297 def real_embed(cls,M):
2298 """
2299 Embed the n-by-n quaternion matrix ``M`` into the space of real
2300 matrices of size 4n-by-4n by first sending each quaternion entry `z
2301 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2302 c+di],[-c + di, a-bi]]`, and then embedding those into a real
2303 matrix.
2304
2305 SETUP::
2306
2307 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2308
2309 EXAMPLES::
2310
2311 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2312 sage: i,j,k = Q.gens()
2313 sage: x = 1 + 2*i + 3*j + 4*k
2314 sage: M = matrix(Q, 1, [[x]])
2315 sage: QuaternionMatrixEJA.real_embed(M)
2316 [ 1 2 3 4]
2317 [-2 1 -4 3]
2318 [-3 4 1 -2]
2319 [-4 -3 2 1]
2320
2321 Embedding is a homomorphism (isomorphism, in fact)::
2322
2323 sage: set_random_seed()
2324 sage: n = ZZ.random_element(2)
2325 sage: Q = QuaternionAlgebra(QQ,-1,-1)
2326 sage: X = random_matrix(Q, n)
2327 sage: Y = random_matrix(Q, n)
2328 sage: Xe = QuaternionMatrixEJA.real_embed(X)
2329 sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2330 sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2331 sage: Xe*Ye == XYe
2332 True
2333
2334 """
2335 super().real_embed(M)
2336 quaternions = M.base_ring()
2337 n = M.nrows()
2338
2339 F = QuadraticField(-1, 'I')
2340 i = F.gen()
2341
2342 blocks = []
2343 for z in M.list():
2344 t = z.coefficient_tuple()
2345 a = t[0]
2346 b = t[1]
2347 c = t[2]
2348 d = t[3]
2349 cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
2350 [-c + d*i, a - b*i]])
2351 realM = ComplexMatrixEJA.real_embed(cplxM)
2352 blocks.append(realM)
2353
2354 # We should have real entries by now, so use the realest field
2355 # we've got for the return value.
2356 return matrix.block(quaternions.base_ring(), n, blocks)
2357
2358
2359
2360 @classmethod
2361 def real_unembed(cls,M):
2362 """
2363 The inverse of _embed_quaternion_matrix().
2364
2365 SETUP::
2366
2367 sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2368
2369 EXAMPLES::
2370
2371 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
2372 ....: [-2, 1, -4, 3],
2373 ....: [-3, 4, 1, -2],
2374 ....: [-4, -3, 2, 1]])
2375 sage: QuaternionMatrixEJA.real_unembed(M)
2376 [1 + 2*i + 3*j + 4*k]
2377
2378 TESTS:
2379
2380 Unembedding is the inverse of embedding::
2381
2382 sage: set_random_seed()
2383 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2384 sage: M = random_matrix(Q, 3)
2385 sage: Me = QuaternionMatrixEJA.real_embed(M)
2386 sage: QuaternionMatrixEJA.real_unembed(Me) == M
2387 True
2388
2389 """
2390 super().real_unembed(M)
2391 n = ZZ(M.nrows())
2392 d = cls.dimension_over_reals()
2393
2394 # Use the base ring of the matrix to ensure that its entries can be
2395 # multiplied by elements of the quaternion algebra.
2396 Q = cls.quaternion_extension(M.base_ring())
2397 i,j,k = Q.gens()
2398
2399 # Go top-left to bottom-right (reading order), converting every
2400 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2401 # quaternion block.
2402 elements = []
2403 for l in range(n/d):
2404 for m in range(n/d):
2405 submat = ComplexMatrixEJA.real_unembed(
2406 M[d*l:d*l+d,d*m:d*m+d] )
2407 if submat[0,0] != submat[1,1].conjugate():
2408 raise ValueError('bad on-diagonal submatrix')
2409 if submat[0,1] != -submat[1,0].conjugate():
2410 raise ValueError('bad off-diagonal submatrix')
2411 z = submat[0,0].real()
2412 z += submat[0,0].imag()*i
2413 z += submat[0,1].real()*j
2414 z += submat[0,1].imag()*k
2415 elements.append(z)
2416
2417 return matrix(Q, n/d, elements)
2418
2419
2420 class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
2421 r"""
2422 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2423 matrices, the usual symmetric Jordan product, and the
2424 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2425 the reals.
2426
2427 SETUP::
2428
2429 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2430
2431 EXAMPLES:
2432
2433 In theory, our "field" can be any subfield of the reals::
2434
2435 sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2436 Euclidean Jordan algebra of dimension 6 over Real Double Field
2437 sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2438 Euclidean Jordan algebra of dimension 6 over Real Field with
2439 53 bits of precision
2440
2441 TESTS:
2442
2443 The dimension of this algebra is `2*n^2 - n`::
2444
2445 sage: set_random_seed()
2446 sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2447 sage: n = ZZ.random_element(1, n_max)
2448 sage: J = QuaternionHermitianEJA(n)
2449 sage: J.dimension() == 2*(n^2) - n
2450 True
2451
2452 The Jordan multiplication is what we think it is::
2453
2454 sage: set_random_seed()
2455 sage: J = QuaternionHermitianEJA.random_instance()
2456 sage: x,y = J.random_elements(2)
2457 sage: actual = (x*y).to_matrix()
2458 sage: X = x.to_matrix()
2459 sage: Y = y.to_matrix()
2460 sage: expected = (X*Y + Y*X)/2
2461 sage: actual == expected
2462 True
2463 sage: J(expected) == x*y
2464 True
2465
2466 We can change the generator prefix::
2467
2468 sage: QuaternionHermitianEJA(2, prefix='a').gens()
2469 (a0, a1, a2, a3, a4, a5)
2470
2471 We can construct the (trivial) algebra of rank zero::
2472
2473 sage: QuaternionHermitianEJA(0)
2474 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2475
2476 """
2477 @classmethod
2478 def _denormalized_basis(cls, n, field):
2479 """
2480 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2481
2482 Why do we embed these? Basically, because all of numerical
2483 linear algebra assumes that you're working with vectors consisting
2484 of `n` entries from a field and scalars from the same field. There's
2485 no way to tell SageMath that (for example) the vectors contain
2486 complex numbers, while the scalar field is real.
2487
2488 SETUP::
2489
2490 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2491
2492 TESTS::
2493
2494 sage: set_random_seed()
2495 sage: n = ZZ.random_element(1,5)
2496 sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)
2497 sage: all( M.is_symmetric() for M in B )
2498 True
2499
2500 """
2501 Q = QuaternionAlgebra(QQ,-1,-1)
2502 I,J,K = Q.gens()
2503
2504 # This is like the symmetric case, but we need to be careful:
2505 #
2506 # * We want conjugate-symmetry, not just symmetry.
2507 # * The diagonal will (as a result) be real.
2508 #
2509 S = []
2510 Eij = matrix.zero(Q,n)
2511 for i in range(n):
2512 for j in range(i+1):
2513 # "build" E_ij
2514 Eij[i,j] = 1
2515 if i == j:
2516 Sij = cls.real_embed(Eij)
2517 S.append(Sij)
2518 else:
2519 # The second, third, and fourth ones have a minus
2520 # because they're conjugated.
2521 # Eij = Eij + Eij.transpose()
2522 Eij[j,i] = 1
2523 Sij_real = cls.real_embed(Eij)
2524 S.append(Sij_real)
2525 # Eij = I*(Eij - Eij.transpose())
2526 Eij[i,j] = I
2527 Eij[j,i] = -I
2528 Sij_I = cls.real_embed(Eij)
2529 S.append(Sij_I)
2530 # Eij = J*(Eij - Eij.transpose())
2531 Eij[i,j] = J
2532 Eij[j,i] = -J
2533 Sij_J = cls.real_embed(Eij)
2534 S.append(Sij_J)
2535 # Eij = K*(Eij - Eij.transpose())
2536 Eij[i,j] = K
2537 Eij[j,i] = -K
2538 Sij_K = cls.real_embed(Eij)
2539 S.append(Sij_K)
2540 Eij[j,i] = 0
2541 # "erase" E_ij
2542 Eij[i,j] = 0
2543
2544 # Since we embedded the entries, we can drop back to the
2545 # desired real "field" instead of the quaternion algebra "Q".
2546 return tuple( s.change_ring(field) for s in S )
2547
2548
2549 def __init__(self, n, field=AA, **kwargs):
2550 # We know this is a valid EJA, but will double-check
2551 # if the user passes check_axioms=True.
2552 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2553
2554 associative = False
2555 if n <= 1:
2556 associative = True
2557
2558 super().__init__(self._denormalized_basis(n,field),
2559 self.jordan_product,
2560 self.trace_inner_product,
2561 field=field,
2562 associative=associative,
2563 **kwargs)
2564
2565 # TODO: this could be factored out somehow, but is left here
2566 # because the MatrixEJA is not presently a subclass of the
2567 # FDEJA class that defines rank() and one().
2568 self.rank.set_cache(n)
2569 idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
2570 self.one.set_cache(self(idV))
2571
2572
2573 @staticmethod
2574 def _max_random_instance_size():
2575 r"""
2576 The maximum rank of a random QuaternionHermitianEJA.
2577 """
2578 return 2 # Dimension 6
2579
2580 @classmethod
2581 def random_instance(cls, **kwargs):
2582 """
2583 Return a random instance of this type of algebra.
2584 """
2585 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2586 return cls(n, **kwargs)
2587
2588 class OctonionHermitianEJA(FiniteDimensionalEJA, MatrixEJA):
2589 r"""
2590
2591 SETUP::
2592
2593 sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2594
2595 EXAMPLES:
2596
2597 The 3-by-3 algebra satisfies the axioms of an EJA::
2598
2599 sage: OctonionHermitianEJA(3, # long time
2600 ....: field=QQ, # long time
2601 ....: orthonormalize=False, # long time
2602 ....: check_axioms=True) # long time
2603 Euclidean Jordan algebra of dimension 27 over Rational Field
2604
2605 TESTS::
2606
2607 sage: J = OctonionHermitianEJA(3, # long time
2608 ....: field=QQ, # long time
2609 ....: orthonormalize=False) # long time
2610 sage: J.one.clear_cache() # long time
2611 sage: J.one() # long time
2612 b0 + b9 + b26
2613 sage: J.one().to_matrix() # long time
2614 +----+----+----+
2615 | e0 | 0 | 0 |
2616 +----+----+----+
2617 | 0 | e0 | 0 |
2618 +----+----+----+
2619 | 0 | 0 | e0 |
2620 +----+----+----+
2621
2622 """
2623 def __init__(self, n, field=AA, **kwargs):
2624 if n > 3:
2625 # Otherwise we don't get an EJA.
2626 raise ValueError("n cannot exceed 3")
2627
2628 # We know this is a valid EJA, but will double-check
2629 # if the user passes check_axioms=True.
2630 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2631
2632 super().__init__(self._denormalized_basis(n,field),
2633 self.jordan_product,
2634 self.trace_inner_product,
2635 field=field,
2636 **kwargs)
2637
2638 # TODO: this could be factored out somehow, but is left here
2639 # because the MatrixEJA is not presently a subclass of the
2640 # FDEJA class that defines rank() and one().
2641 self.rank.set_cache(n)
2642 idV = self.matrix_space().one()
2643 self.one.set_cache(self(idV))
2644
2645
2646 @classmethod
2647 def _denormalized_basis(cls, n, field):
2648 """
2649 Returns a basis for the space of octonion Hermitian n-by-n
2650 matrices.
2651
2652 SETUP::
2653
2654 sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2655
2656 EXAMPLES::
2657
2658 sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ)
2659 sage: all( M.is_hermitian() for M in B )
2660 True
2661 sage: len(B)
2662 27
2663
2664 """
2665 from mjo.octonions import OctonionMatrixAlgebra
2666 MS = OctonionMatrixAlgebra(n, scalars=field)
2667 es = MS.entry_algebra().gens()
2668
2669 basis = []
2670 for i in range(n):
2671 for j in range(i+1):
2672 if i == j:
2673 E_ii = MS.monomial( (i,j,es[0]) )
2674 basis.append(E_ii)
2675 else:
2676 for e in es:
2677 E_ij = MS.monomial( (i,j,e) )
2678 ec = e.conjugate()
2679 # If the conjugate has a negative sign in front
2680 # of it, (j,i,ec) won't be a monomial!
2681 if (j,i,ec) in MS.indices():
2682 E_ij += MS.monomial( (j,i,ec) )
2683 else:
2684 E_ij -= MS.monomial( (j,i,-ec) )
2685 basis.append(E_ij)
2686
2687 return tuple( basis )
2688
2689
2690 class HadamardEJA(ConcreteEJA):
2691 """
2692 Return the Euclidean Jordan Algebra corresponding to the set
2693 `R^n` under the Hadamard product.
2694
2695 Note: this is nothing more than the Cartesian product of ``n``
2696 copies of the spin algebra. Once Cartesian product algebras
2697 are implemented, this can go.
2698
2699 SETUP::
2700
2701 sage: from mjo.eja.eja_algebra import HadamardEJA
2702
2703 EXAMPLES:
2704
2705 This multiplication table can be verified by hand::
2706
2707 sage: J = HadamardEJA(3)
2708 sage: b0,b1,b2 = J.gens()
2709 sage: b0*b0
2710 b0
2711 sage: b0*b1
2712 0
2713 sage: b0*b2
2714 0
2715 sage: b1*b1
2716 b1
2717 sage: b1*b2
2718 0
2719 sage: b2*b2
2720 b2
2721
2722 TESTS:
2723
2724 We can change the generator prefix::
2725
2726 sage: HadamardEJA(3, prefix='r').gens()
2727 (r0, r1, r2)
2728
2729 """
2730 def __init__(self, n, field=AA, **kwargs):
2731 if n == 0:
2732 jordan_product = lambda x,y: x
2733 inner_product = lambda x,y: x
2734 else:
2735 def jordan_product(x,y):
2736 P = x.parent()
2737 return P( xi*yi for (xi,yi) in zip(x,y) )
2738
2739 def inner_product(x,y):
2740 return (x.T*y)[0,0]
2741
2742 # New defaults for keyword arguments. Don't orthonormalize
2743 # because our basis is already orthonormal with respect to our
2744 # inner-product. Don't check the axioms, because we know this
2745 # is a valid EJA... but do double-check if the user passes
2746 # check_axioms=True. Note: we DON'T override the "check_field"
2747 # default here, because the user can pass in a field!
2748 if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
2749 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2750
2751 column_basis = tuple( b.column()
2752 for b in FreeModule(field, n).basis() )
2753 super().__init__(column_basis,
2754 jordan_product,
2755 inner_product,
2756 field=field,
2757 associative=True,
2758 **kwargs)
2759 self.rank.set_cache(n)
2760
2761 if n == 0:
2762 self.one.set_cache( self.zero() )
2763 else:
2764 self.one.set_cache( sum(self.gens()) )
2765
2766 @staticmethod
2767 def _max_random_instance_size():
2768 r"""
2769 The maximum dimension of a random HadamardEJA.
2770 """
2771 return 5
2772
2773 @classmethod
2774 def random_instance(cls, **kwargs):
2775 """
2776 Return a random instance of this type of algebra.
2777 """
2778 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2779 return cls(n, **kwargs)
2780
2781
2782 class BilinearFormEJA(ConcreteEJA):
2783 r"""
2784 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2785 with the half-trace inner product and jordan product ``x*y =
2786 (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2787 a symmetric positive-definite "bilinear form" matrix. Its
2788 dimension is the size of `B`, and it has rank two in dimensions
2789 larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2790 the identity matrix of order ``n``.
2791
2792 We insist that the one-by-one upper-left identity block of `B` be
2793 passed in as well so that we can be passed a matrix of size zero
2794 to construct a trivial algebra.
2795
2796 SETUP::
2797
2798 sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2799 ....: JordanSpinEJA)
2800
2801 EXAMPLES:
2802
2803 When no bilinear form is specified, the identity matrix is used,
2804 and the resulting algebra is the Jordan spin algebra::
2805
2806 sage: B = matrix.identity(AA,3)
2807 sage: J0 = BilinearFormEJA(B)
2808 sage: J1 = JordanSpinEJA(3)
2809 sage: J0.multiplication_table() == J0.multiplication_table()
2810 True
2811
2812 An error is raised if the matrix `B` does not correspond to a
2813 positive-definite bilinear form::
2814
2815 sage: B = matrix.random(QQ,2,3)
2816 sage: J = BilinearFormEJA(B)
2817 Traceback (most recent call last):
2818 ...
2819 ValueError: bilinear form is not positive-definite
2820 sage: B = matrix.zero(QQ,3)
2821 sage: J = BilinearFormEJA(B)
2822 Traceback (most recent call last):
2823 ...
2824 ValueError: bilinear form is not positive-definite
2825
2826 TESTS:
2827
2828 We can create a zero-dimensional algebra::
2829
2830 sage: B = matrix.identity(AA,0)
2831 sage: J = BilinearFormEJA(B)
2832 sage: J.basis()
2833 Finite family {}
2834
2835 We can check the multiplication condition given in the Jordan, von
2836 Neumann, and Wigner paper (and also discussed on my "On the
2837 symmetry..." paper). Note that this relies heavily on the standard
2838 choice of basis, as does anything utilizing the bilinear form
2839 matrix. We opt not to orthonormalize the basis, because if we
2840 did, we would have to normalize the `s_{i}` in a similar manner::
2841
2842 sage: set_random_seed()
2843 sage: n = ZZ.random_element(5)
2844 sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2845 sage: B11 = matrix.identity(QQ,1)
2846 sage: B22 = M.transpose()*M
2847 sage: B = block_matrix(2,2,[ [B11,0 ],
2848 ....: [0, B22 ] ])
2849 sage: J = BilinearFormEJA(B, orthonormalize=False)
2850 sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2851 sage: V = J.vector_space()
2852 sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2853 ....: for ei in eis ]
2854 sage: actual = [ sis[i]*sis[j]
2855 ....: for i in range(n-1)
2856 ....: for j in range(n-1) ]
2857 sage: expected = [ J.one() if i == j else J.zero()
2858 ....: for i in range(n-1)
2859 ....: for j in range(n-1) ]
2860 sage: actual == expected
2861 True
2862
2863 """
2864 def __init__(self, B, field=AA, **kwargs):
2865 # The matrix "B" is supplied by the user in most cases,
2866 # so it makes sense to check whether or not its positive-
2867 # definite unless we are specifically asked not to...
2868 if ("check_axioms" not in kwargs) or kwargs["check_axioms"]:
2869 if not B.is_positive_definite():
2870 raise ValueError("bilinear form is not positive-definite")
2871
2872 # However, all of the other data for this EJA is computed
2873 # by us in manner that guarantees the axioms are
2874 # satisfied. So, again, unless we are specifically asked to
2875 # verify things, we'll skip the rest of the checks.
2876 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2877
2878 def inner_product(x,y):
2879 return (y.T*B*x)[0,0]
2880
2881 def jordan_product(x,y):
2882 P = x.parent()
2883 x0 = x[0,0]
2884 xbar = x[1:,0]
2885 y0 = y[0,0]
2886 ybar = y[1:,0]
2887 z0 = inner_product(y,x)
2888 zbar = y0*xbar + x0*ybar
2889 return P([z0] + zbar.list())
2890
2891 n = B.nrows()
2892 column_basis = tuple( b.column()
2893 for b in FreeModule(field, n).basis() )
2894
2895 # TODO: I haven't actually checked this, but it seems legit.
2896 associative = False
2897 if n <= 2:
2898 associative = True
2899
2900 super().__init__(column_basis,
2901 jordan_product,
2902 inner_product,
2903 field=field,
2904 associative=associative,
2905 **kwargs)
2906
2907 # The rank of this algebra is two, unless we're in a
2908 # one-dimensional ambient space (because the rank is bounded
2909 # by the ambient dimension).
2910 self.rank.set_cache(min(n,2))
2911
2912 if n == 0:
2913 self.one.set_cache( self.zero() )
2914 else:
2915 self.one.set_cache( self.monomial(0) )
2916
2917 @staticmethod
2918 def _max_random_instance_size():
2919 r"""
2920 The maximum dimension of a random BilinearFormEJA.
2921 """
2922 return 5
2923
2924 @classmethod
2925 def random_instance(cls, **kwargs):
2926 """
2927 Return a random instance of this algebra.
2928 """
2929 n = ZZ.random_element(cls._max_random_instance_size() + 1)
2930 if n.is_zero():
2931 B = matrix.identity(ZZ, n)
2932 return cls(B, **kwargs)
2933
2934 B11 = matrix.identity(ZZ, 1)
2935 M = matrix.random(ZZ, n-1)
2936 I = matrix.identity(ZZ, n-1)
2937 alpha = ZZ.zero()
2938 while alpha.is_zero():
2939 alpha = ZZ.random_element().abs()
2940 B22 = M.transpose()*M + alpha*I
2941
2942 from sage.matrix.special import block_matrix
2943 B = block_matrix(2,2, [ [B11, ZZ(0) ],
2944 [ZZ(0), B22 ] ])
2945
2946 return cls(B, **kwargs)
2947
2948
2949 class JordanSpinEJA(BilinearFormEJA):
2950 """
2951 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2952 with the usual inner product and jordan product ``x*y =
2953 (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2954 the reals.
2955
2956 SETUP::
2957
2958 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2959
2960 EXAMPLES:
2961
2962 This multiplication table can be verified by hand::
2963
2964 sage: J = JordanSpinEJA(4)
2965 sage: b0,b1,b2,b3 = J.gens()
2966 sage: b0*b0
2967 b0
2968 sage: b0*b1
2969 b1
2970 sage: b0*b2
2971 b2
2972 sage: b0*b3
2973 b3
2974 sage: b1*b2
2975 0
2976 sage: b1*b3
2977 0
2978 sage: b2*b3
2979 0
2980
2981 We can change the generator prefix::
2982
2983 sage: JordanSpinEJA(2, prefix='B').gens()
2984 (B0, B1)
2985
2986 TESTS:
2987
2988 Ensure that we have the usual inner product on `R^n`::
2989
2990 sage: set_random_seed()
2991 sage: J = JordanSpinEJA.random_instance()
2992 sage: x,y = J.random_elements(2)
2993 sage: actual = x.inner_product(y)
2994 sage: expected = x.to_vector().inner_product(y.to_vector())
2995 sage: actual == expected
2996 True
2997
2998 """
2999 def __init__(self, n, *args, **kwargs):
3000 # This is a special case of the BilinearFormEJA with the
3001 # identity matrix as its bilinear form.
3002 B = matrix.identity(ZZ, n)
3003
3004 # Don't orthonormalize because our basis is already
3005 # orthonormal with respect to our inner-product.
3006 if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
3007
3008 # But also don't pass check_field=False here, because the user
3009 # can pass in a field!
3010 super().__init__(B, *args, **kwargs)
3011
3012 @staticmethod
3013 def _max_random_instance_size():
3014 r"""
3015 The maximum dimension of a random JordanSpinEJA.
3016 """
3017 return 5
3018
3019 @classmethod
3020 def random_instance(cls, **kwargs):
3021 """
3022 Return a random instance of this type of algebra.
3023
3024 Needed here to override the implementation for ``BilinearFormEJA``.
3025 """
3026 n = ZZ.random_element(cls._max_random_instance_size() + 1)
3027 return cls(n, **kwargs)
3028
3029
3030 class TrivialEJA(ConcreteEJA):
3031 """
3032 The trivial Euclidean Jordan algebra consisting of only a zero element.
3033
3034 SETUP::
3035
3036 sage: from mjo.eja.eja_algebra import TrivialEJA
3037
3038 EXAMPLES::
3039
3040 sage: J = TrivialEJA()
3041 sage: J.dimension()
3042 0
3043 sage: J.zero()
3044 0
3045 sage: J.one()
3046 0
3047 sage: 7*J.one()*12*J.one()
3048 0
3049 sage: J.one().inner_product(J.one())
3050 0
3051 sage: J.one().norm()
3052 0
3053 sage: J.one().subalgebra_generated_by()
3054 Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
3055 sage: J.rank()
3056 0
3057
3058 """
3059 def __init__(self, **kwargs):
3060 jordan_product = lambda x,y: x
3061 inner_product = lambda x,y: 0
3062 basis = ()
3063
3064 # New defaults for keyword arguments
3065 if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
3066 if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
3067
3068 super().__init__(basis,
3069 jordan_product,
3070 inner_product,
3071 associative=True,
3072 **kwargs)
3073
3074 # The rank is zero using my definition, namely the dimension of the
3075 # largest subalgebra generated by any element.
3076 self.rank.set_cache(0)
3077 self.one.set_cache( self.zero() )
3078
3079 @classmethod
3080 def random_instance(cls, **kwargs):
3081 # We don't take a "size" argument so the superclass method is
3082 # inappropriate for us.
3083 return cls(**kwargs)
3084
3085
3086 class CartesianProductEJA(FiniteDimensionalEJA):
3087 r"""
3088 The external (orthogonal) direct sum of two or more Euclidean
3089 Jordan algebras. Every Euclidean Jordan algebra decomposes into an
3090 orthogonal direct sum of simple Euclidean Jordan algebras which is
3091 then isometric to a Cartesian product, so no generality is lost by
3092 providing only this construction.
3093
3094 SETUP::
3095
3096 sage: from mjo.eja.eja_algebra import (random_eja,
3097 ....: CartesianProductEJA,
3098 ....: HadamardEJA,
3099 ....: JordanSpinEJA,
3100 ....: RealSymmetricEJA)
3101
3102 EXAMPLES:
3103
3104 The Jordan product is inherited from our factors and implemented by
3105 our CombinatorialFreeModule Cartesian product superclass::
3106
3107 sage: set_random_seed()
3108 sage: J1 = HadamardEJA(2)
3109 sage: J2 = RealSymmetricEJA(2)
3110 sage: J = cartesian_product([J1,J2])
3111 sage: x,y = J.random_elements(2)
3112 sage: x*y in J
3113 True
3114
3115 The ability to retrieve the original factors is implemented by our
3116 CombinatorialFreeModule Cartesian product superclass::
3117
3118 sage: J1 = HadamardEJA(2, field=QQ)
3119 sage: J2 = JordanSpinEJA(3, field=QQ)
3120 sage: J = cartesian_product([J1,J2])
3121 sage: J.cartesian_factors()
3122 (Euclidean Jordan algebra of dimension 2 over Rational Field,
3123 Euclidean Jordan algebra of dimension 3 over Rational Field)
3124
3125 You can provide more than two factors::
3126
3127 sage: J1 = HadamardEJA(2)
3128 sage: J2 = JordanSpinEJA(3)
3129 sage: J3 = RealSymmetricEJA(3)
3130 sage: cartesian_product([J1,J2,J3])
3131 Euclidean Jordan algebra of dimension 2 over Algebraic Real
3132 Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3133 Real Field (+) Euclidean Jordan algebra of dimension 6 over
3134 Algebraic Real Field
3135
3136 Rank is additive on a Cartesian product::
3137
3138 sage: J1 = HadamardEJA(1)
3139 sage: J2 = RealSymmetricEJA(2)
3140 sage: J = cartesian_product([J1,J2])
3141 sage: J1.rank.clear_cache()
3142 sage: J2.rank.clear_cache()
3143 sage: J.rank.clear_cache()
3144 sage: J.rank()
3145 3
3146 sage: J.rank() == J1.rank() + J2.rank()
3147 True
3148
3149 The same rank computation works over the rationals, with whatever
3150 basis you like::
3151
3152 sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3153 sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3154 sage: J = cartesian_product([J1,J2])
3155 sage: J1.rank.clear_cache()
3156 sage: J2.rank.clear_cache()
3157 sage: J.rank.clear_cache()
3158 sage: J.rank()
3159 3
3160 sage: J.rank() == J1.rank() + J2.rank()
3161 True
3162
3163 The product algebra will be associative if and only if all of its
3164 components are associative::
3165
3166 sage: J1 = HadamardEJA(2)
3167 sage: J1.is_associative()
3168 True
3169 sage: J2 = HadamardEJA(3)
3170 sage: J2.is_associative()
3171 True
3172 sage: J3 = RealSymmetricEJA(3)
3173 sage: J3.is_associative()
3174 False
3175 sage: CP1 = cartesian_product([J1,J2])
3176 sage: CP1.is_associative()
3177 True
3178 sage: CP2 = cartesian_product([J1,J3])
3179 sage: CP2.is_associative()
3180 False
3181
3182 Cartesian products of Cartesian products work::
3183
3184 sage: J1 = JordanSpinEJA(1)
3185 sage: J2 = JordanSpinEJA(1)
3186 sage: J3 = JordanSpinEJA(1)
3187 sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3188 sage: J.multiplication_table()
3189 +----++----+----+----+
3190 | * || b0 | b1 | b2 |
3191 +====++====+====+====+
3192 | b0 || b0 | 0 | 0 |
3193 +----++----+----+----+
3194 | b1 || 0 | b1 | 0 |
3195 +----++----+----+----+
3196 | b2 || 0 | 0 | b2 |
3197 +----++----+----+----+
3198 sage: HadamardEJA(3).multiplication_table()
3199 +----++----+----+----+
3200 | * || b0 | b1 | b2 |
3201 +====++====+====+====+
3202 | b0 || b0 | 0 | 0 |
3203 +----++----+----+----+
3204 | b1 || 0 | b1 | 0 |
3205 +----++----+----+----+
3206 | b2 || 0 | 0 | b2 |
3207 +----++----+----+----+
3208
3209 TESTS:
3210
3211 All factors must share the same base field::
3212
3213 sage: J1 = HadamardEJA(2, field=QQ)
3214 sage: J2 = RealSymmetricEJA(2)
3215 sage: CartesianProductEJA((J1,J2))
3216 Traceback (most recent call last):
3217 ...
3218 ValueError: all factors must share the same base field
3219
3220 The cached unit element is the same one that would be computed::
3221
3222 sage: set_random_seed() # long time
3223 sage: J1 = random_eja() # long time
3224 sage: J2 = random_eja() # long time
3225 sage: J = cartesian_product([J1,J2]) # long time
3226 sage: actual = J.one() # long time
3227 sage: J.one.clear_cache() # long time
3228 sage: expected = J.one() # long time
3229 sage: actual == expected # long time
3230 True
3231
3232 """
3233 Element = FiniteDimensionalEJAElement
3234
3235
3236 def __init__(self, factors, **kwargs):
3237 m = len(factors)
3238 if m == 0:
3239 return TrivialEJA()
3240
3241 self._sets = factors
3242
3243 field = factors[0].base_ring()
3244 if not all( J.base_ring() == field for J in factors ):
3245 raise ValueError("all factors must share the same base field")
3246
3247 associative = all( f.is_associative() for f in factors )
3248
3249 MS = self.matrix_space()
3250 basis = []
3251 zero = MS.zero()
3252 for i in range(m):
3253 for b in factors[i].matrix_basis():
3254 z = list(zero)
3255 z[i] = b
3256 basis.append(z)
3257
3258 basis = tuple( MS(b) for b in basis )
3259
3260 # Define jordan/inner products that operate on that matrix_basis.
3261 def jordan_product(x,y):
3262 return MS(tuple(
3263 (factors[i](x[i])*factors[i](y[i])).to_matrix()
3264 for i in range(m)
3265 ))
3266
3267 def inner_product(x, y):
3268 return sum(
3269 factors[i](x[i]).inner_product(factors[i](y[i]))
3270 for i in range(m)
3271 )
3272
3273 # There's no need to check the field since it already came
3274 # from an EJA. Likewise the axioms are guaranteed to be
3275 # satisfied, unless the guy writing this class sucks.
3276 #
3277 # If you want the basis to be orthonormalized, orthonormalize
3278 # the factors.
3279 FiniteDimensionalEJA.__init__(self,
3280 basis,
3281 jordan_product,
3282 inner_product,
3283 field=field,
3284 orthonormalize=False,
3285 associative=associative,
3286 cartesian_product=True,
3287 check_field=False,
3288 check_axioms=False)
3289
3290 ones = tuple(J.one().to_matrix() for J in factors)
3291 self.one.set_cache(self(ones))
3292 self.rank.set_cache(sum(J.rank() for J in factors))
3293
3294 def cartesian_factors(self):
3295 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3296 return self._sets
3297
3298 def cartesian_factor(self, i):
3299 r"""
3300 Return the ``i``th factor of this algebra.
3301 """
3302 return self._sets[i]
3303
3304 def _repr_(self):
3305 # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3306 from sage.categories.cartesian_product import cartesian_product
3307 return cartesian_product.symbol.join("%s" % factor
3308 for factor in self._sets)
3309
3310 def matrix_space(self):
3311 r"""
3312 Return the space that our matrix basis lives in as a Cartesian
3313 product.
3314
3315 SETUP::
3316
3317 sage: from mjo.eja.eja_algebra import (HadamardEJA,
3318 ....: RealSymmetricEJA)
3319
3320 EXAMPLES::
3321
3322 sage: J1 = HadamardEJA(1)
3323 sage: J2 = RealSymmetricEJA(2)
3324 sage: J = cartesian_product([J1,J2])
3325 sage: J.matrix_space()
3326 The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3327 matrices over Algebraic Real Field, Full MatrixSpace of 2
3328 by 2 dense matrices over Algebraic Real Field)
3329
3330 """
3331 from sage.categories.cartesian_product import cartesian_product
3332 return cartesian_product( [J.matrix_space()
3333 for J in self.cartesian_factors()] )
3334
3335 @cached_method
3336 def cartesian_projection(self, i):
3337 r"""
3338 SETUP::
3339
3340 sage: from mjo.eja.eja_algebra import (random_eja,
3341 ....: JordanSpinEJA,
3342 ....: HadamardEJA,
3343 ....: RealSymmetricEJA,
3344 ....: ComplexHermitianEJA)
3345
3346 EXAMPLES:
3347
3348 The projection morphisms are Euclidean Jordan algebra
3349 operators::
3350
3351 sage: J1 = HadamardEJA(2)
3352 sage: J2 = RealSymmetricEJA(2)
3353 sage: J = cartesian_product([J1,J2])
3354 sage: J.cartesian_projection(0)
3355 Linear operator between finite-dimensional Euclidean Jordan
3356 algebras represented by the matrix:
3357 [1 0 0 0 0]
3358 [0 1 0 0 0]
3359 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3360 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3361 Algebraic Real Field
3362 Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3363 Real Field
3364 sage: J.cartesian_projection(1)
3365 Linear operator between finite-dimensional Euclidean Jordan
3366 algebras represented by the matrix:
3367 [0 0 1 0 0]
3368 [0 0 0 1 0]
3369 [0 0 0 0 1]
3370 Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3371 Real Field (+) Euclidean Jordan algebra of dimension 3 over
3372 Algebraic Real Field
3373 Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3374 Real Field
3375
3376 The projections work the way you'd expect on the vector
3377 representation of an element::
3378
3379 sage: J1 = JordanSpinEJA(2)
3380 sage: J2 = ComplexHermitianEJA(2)
3381 sage: J = cartesian_product([J1,J2])
3382 sage: pi_left = J.cartesian_projection(0)
3383 sage: pi_right = J.cartesian_projection(1)
3384 sage: pi_left(J.one()).to_vector()
3385 (1, 0)
3386 sage: pi_right(J.one()).to_vector()
3387 (1, 0, 0, 1)
3388 sage: J.one().to_vector()
3389 (1, 0, 1, 0, 0, 1)
3390
3391 TESTS:
3392
3393 The answer never changes::
3394
3395 sage: set_random_seed()
3396 sage: J1 = random_eja()
3397 sage: J2 = random_eja()
3398 sage: J = cartesian_product([J1,J2])
3399 sage: P0 = J.cartesian_projection(0)
3400 sage: P1 = J.cartesian_projection(0)
3401 sage: P0 == P1
3402 True
3403
3404 """
3405 offset = sum( self.cartesian_factor(k).dimension()
3406 for k in range(i) )
3407 Ji = self.cartesian_factor(i)
3408 Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
3409 codomain=Ji)
3410
3411 return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
3412
3413 @cached_method
3414 def cartesian_embedding(self, i):
3415 r"""
3416 SETUP::
3417
3418 sage: from mjo.eja.eja_algebra import (random_eja,
3419 ....: JordanSpinEJA,
3420 ....: HadamardEJA,
3421 ....: RealSymmetricEJA)
3422
3423 EXAMPLES:
3424
3425 The embedding morphisms are Euclidean Jordan algebra
3426 operators::
3427
3428 sage: J1 = HadamardEJA(2)
3429 sage: J2 = RealSymmetricEJA(2)
3430 sage: J = cartesian_product([J1,J2])
3431 sage: J.cartesian_embedding(0)
3432 Linear operator between finite-dimensional Euclidean Jordan
3433 algebras represented by the matrix:
3434 [1 0]
3435 [0 1]
3436 [0 0]
3437 [0 0]
3438 [0 0]
3439 Domain: Euclidean Jordan algebra of dimension 2 over
3440 Algebraic Real Field
3441 Codomain: Euclidean Jordan algebra of dimension 2 over
3442 Algebraic Real Field (+) Euclidean Jordan algebra of
3443 dimension 3 over Algebraic Real Field
3444 sage: J.cartesian_embedding(1)
3445 Linear operator between finite-dimensional Euclidean Jordan
3446 algebras represented by the matrix:
3447 [0 0 0]
3448 [0 0 0]
3449 [1 0 0]
3450 [0 1 0]
3451 [0 0 1]
3452 Domain: Euclidean Jordan algebra of dimension 3 over
3453 Algebraic Real Field
3454 Codomain: Euclidean Jordan algebra of dimension 2 over
3455 Algebraic Real Field (+) Euclidean Jordan algebra of
3456 dimension 3 over Algebraic Real Field
3457
3458 The embeddings work the way you'd expect on the vector
3459 representation of an element::
3460
3461 sage: J1 = JordanSpinEJA(3)
3462 sage: J2 = RealSymmetricEJA(2)
3463 sage: J = cartesian_product([J1,J2])
3464 sage: iota_left = J.cartesian_embedding(0)
3465 sage: iota_right = J.cartesian_embedding(1)
3466 sage: iota_left(J1.zero()) == J.zero()
3467 True
3468 sage: iota_right(J2.zero()) == J.zero()
3469 True
3470 sage: J1.one().to_vector()
3471 (1, 0, 0)
3472 sage: iota_left(J1.one()).to_vector()
3473 (1, 0, 0, 0, 0, 0)
3474 sage: J2.one().to_vector()
3475 (1, 0, 1)
3476 sage: iota_right(J2.one()).to_vector()
3477 (0, 0, 0, 1, 0, 1)
3478 sage: J.one().to_vector()
3479 (1, 0, 0, 1, 0, 1)
3480
3481 TESTS:
3482
3483 The answer never changes::
3484
3485 sage: set_random_seed()
3486 sage: J1 = random_eja()
3487 sage: J2 = random_eja()
3488 sage: J = cartesian_product([J1,J2])
3489 sage: E0 = J.cartesian_embedding(0)
3490 sage: E1 = J.cartesian_embedding(0)
3491 sage: E0 == E1
3492 True
3493
3494 Composing a projection with the corresponding inclusion should
3495 produce the identity map, and mismatching them should produce
3496 the zero map::
3497
3498 sage: set_random_seed()
3499 sage: J1 = random_eja()
3500 sage: J2 = random_eja()
3501 sage: J = cartesian_product([J1,J2])
3502 sage: iota_left = J.cartesian_embedding(0)
3503 sage: iota_right = J.cartesian_embedding(1)
3504 sage: pi_left = J.cartesian_projection(0)
3505 sage: pi_right = J.cartesian_projection(1)
3506 sage: pi_left*iota_left == J1.one().operator()
3507 True
3508 sage: pi_right*iota_right == J2.one().operator()
3509 True
3510 sage: (pi_left*iota_right).is_zero()
3511 True
3512 sage: (pi_right*iota_left).is_zero()
3513 True
3514
3515 """
3516 offset = sum( self.cartesian_factor(k).dimension()
3517 for k in range(i) )
3518 Ji = self.cartesian_factor(i)
3519 Ei = Ji._module_morphism(lambda j: self.monomial(j + offset),
3520 codomain=self)
3521 return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
3522
3523
3524
3525 FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
3526
3527 class RationalBasisCartesianProductEJA(CartesianProductEJA,
3528 RationalBasisEJA):
3529 r"""
3530 A separate class for products of algebras for which we know a
3531 rational basis.
3532
3533 SETUP::
3534
3535 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3536 ....: RealSymmetricEJA)
3537
3538 EXAMPLES:
3539
3540 This gives us fast characteristic polynomial computations in
3541 product algebras, too::
3542
3543
3544 sage: J1 = JordanSpinEJA(2)
3545 sage: J2 = RealSymmetricEJA(3)
3546 sage: J = cartesian_product([J1,J2])
3547 sage: J.characteristic_polynomial_of().degree()
3548 5
3549 sage: J.rank()
3550 5
3551
3552 """
3553 def __init__(self, algebras, **kwargs):
3554 CartesianProductEJA.__init__(self, algebras, **kwargs)
3555
3556 self._rational_algebra = None
3557 if self.vector_space().base_field() is not QQ:
3558 self._rational_algebra = cartesian_product([
3559 r._rational_algebra for r in algebras
3560 ])
3561
3562
3563 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
3564
3565 def random_eja(*args, **kwargs):
3566 J1 = ConcreteEJA.random_instance(*args, **kwargs)
3567
3568 # This might make Cartesian products appear roughly as often as
3569 # any other ConcreteEJA.
3570 if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
3571 # Use random_eja() again so we can get more than two factors.
3572 J2 = random_eja(*args, **kwargs)
3573 J = cartesian_product([J1,J2])
3574 return J
3575 else:
3576 return J1