]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/eja/eja_algebra.py
projects / sage.d.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
history | raw | HEAD
eja: clean up imports.
[sage.d.git] / mjo / eja / eja_algebra.py
1 """
2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
6 """
7
8
9
10 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra
11 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element import FiniteDimensionalAlgebraElement
12 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
13 from sage.categories.finite_dimensional_algebras_with_basis import FiniteDimensionalAlgebrasWithBasis
14 from sage.functions.other import sqrt
15 from sage.matrix.constructor import matrix
16 from sage.misc.cachefunc import cached_method
17 from sage.misc.prandom import choice
18 from sage.modules.free_module import VectorSpace
19 from sage.modules.free_module_element import vector
20 from sage.rings.integer_ring import ZZ
21 from sage.rings.number_field.number_field import QuadraticField
22 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
23 from sage.rings.rational_field import QQ
24 from sage.structure.element import is_Matrix
25 from sage.structure.category_object import normalize_names
26
27 from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
28
29
30 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
31 @staticmethod
32 def __classcall_private__(cls,
33 field,
34 mult_table,
35 names='e',
36 assume_associative=False,
37 category=None,
38 rank=None,
39 natural_basis=None):
40 n = len(mult_table)
41 mult_table = [b.base_extend(field) for b in mult_table]
42 for b in mult_table:
43 b.set_immutable()
44 if not (is_Matrix(b) and b.dimensions() == (n, n)):
45 raise ValueError("input is not a multiplication table")
46 mult_table = tuple(mult_table)
47
48 cat = FiniteDimensionalAlgebrasWithBasis(field)
49 cat.or_subcategory(category)
50 if assume_associative:
51 cat = cat.Associative()
52
53 names = normalize_names(n, names)
54
55 fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls)
56 return fda.__classcall__(cls,
57 field,
58 mult_table,
59 assume_associative=assume_associative,
60 names=names,
61 category=cat,
62 rank=rank,
63 natural_basis=natural_basis)
64
65
66 def __init__(self,
67 field,
68 mult_table,
69 names='e',
70 assume_associative=False,
71 category=None,
72 rank=None,
73 natural_basis=None):
74 """
75 SETUP::
76
77 sage: from mjo.eja.eja_algebra import random_eja
78
79 EXAMPLES:
80
81 By definition, Jordan multiplication commutes::
82
83 sage: set_random_seed()
84 sage: J = random_eja()
85 sage: x = J.random_element()
86 sage: y = J.random_element()
87 sage: x*y == y*x
88 True
89
90 """
91 self._rank = rank
92 self._natural_basis = natural_basis
93 self._multiplication_table = mult_table
94 fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
95 fda.__init__(field,
96 mult_table,
97 names=names,
98 category=category)
99
100
101 def _repr_(self):
102 """
103 Return a string representation of ``self``.
104 """
105 fmt = "Euclidean Jordan algebra of degree {} over {}"
106 return fmt.format(self.degree(), self.base_ring())
107
108
109 def _a_regular_element(self):
110 """
111 Guess a regular element. Needed to compute the basis for our
112 characteristic polynomial coefficients.
113 """
114 gs = self.gens()
115 z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
116 if not z.is_regular():
117 raise ValueError("don't know a regular element")
118 return z
119
120
121 @cached_method
122 def _charpoly_basis_space(self):
123 """
124 Return the vector space spanned by the basis used in our
125 characteristic polynomial coefficients. This is used not only to
126 compute those coefficients, but also any time we need to
127 evaluate the coefficients (like when we compute the trace or
128 determinant).
129 """
130 z = self._a_regular_element()
131 V = self.vector_space()
132 V1 = V.span_of_basis( (z**k).vector() for k in range(self.rank()) )
133 b = (V1.basis() + V1.complement().basis())
134 return V.span_of_basis(b)
135
136
137 @cached_method
138 def _charpoly_coeff(self, i):
139 """
140 Return the coefficient polynomial "a_{i}" of this algebra's
141 general characteristic polynomial.
142
143 Having this be a separate cached method lets us compute and
144 store the trace/determinant (a_{r-1} and a_{0} respectively)
145 separate from the entire characteristic polynomial.
146 """
147 (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
148 R = A_of_x.base_ring()
149 if i >= self.rank():
150 # Guaranteed by theory
151 return R.zero()
152
153 # Danger: the in-place modification is done for performance
154 # reasons (reconstructing a matrix with huge polynomial
155 # entries is slow), but I don't know how cached_method works,
156 # so it's highly possible that we're modifying some global
157 # list variable by reference, here. In other words, you
158 # probably shouldn't call this method twice on the same
159 # algebra, at the same time, in two threads
160 Ai_orig = A_of_x.column(i)
161 A_of_x.set_column(i,xr)
162 numerator = A_of_x.det()
163 A_of_x.set_column(i,Ai_orig)
164
165 # We're relying on the theory here to ensure that each a_i is
166 # indeed back in R, and the added negative signs are to make
167 # the whole charpoly expression sum to zero.
168 return R(-numerator/detA)
169
170
171 @cached_method
172 def _charpoly_matrix_system(self):
173 """
174 Compute the matrix whose entries A_ij are polynomials in
175 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
176 corresponding to `x^r` and the determinent of the matrix A =
177 [A_ij]. In other words, all of the fixed (cachable) data needed
178 to compute the coefficients of the characteristic polynomial.
179 """
180 r = self.rank()
181 n = self.dimension()
182
183 # Construct a new algebra over a multivariate polynomial ring...
184 names = ['X' + str(i) for i in range(1,n+1)]
185 R = PolynomialRing(self.base_ring(), names)
186 J = FiniteDimensionalEuclideanJordanAlgebra(R,
187 self._multiplication_table,
188 rank=r)
189
190 idmat = matrix.identity(J.base_ring(), n)
191
192 W = self._charpoly_basis_space()
193 W = W.change_ring(R.fraction_field())
194
195 # Starting with the standard coordinates x = (X1,X2,...,Xn)
196 # and then converting the entries to W-coordinates allows us
197 # to pass in the standard coordinates to the charpoly and get
198 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
199 # we have
200 #
201 # W.coordinates(x^2) eval'd at (standard z-coords)
202 # =
203 # W-coords of (z^2)
204 # =
205 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
206 #
207 # We want the middle equivalent thing in our matrix, but use
208 # the first equivalent thing instead so that we can pass in
209 # standard coordinates.
210 x = J(W(R.gens()))
211 l1 = [matrix.column(W.coordinates((x**k).vector())) for k in range(r)]
212 l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)]
213 A_of_x = matrix.block(R, 1, n, (l1 + l2))
214 xr = W.coordinates((x**r).vector())
215 return (A_of_x, x, xr, A_of_x.det())
216
217
218 @cached_method
219 def characteristic_polynomial(self):
220 """
221
222 .. WARNING::
223
224 This implementation doesn't guarantee that the polynomial
225 denominator in the coefficients is not identically zero, so
226 theoretically it could crash. The way that this is handled
227 in e.g. Faraut and Koranyi is to use a basis that guarantees
228 the denominator is non-zero. But, doing so requires knowledge
229 of at least one regular element, and we don't even know how
230 to do that. The trade-off is that, if we use the standard basis,
231 the resulting polynomial will accept the "usual" coordinates. In
232 other words, we don't have to do a change of basis before e.g.
233 computing the trace or determinant.
234
235 SETUP::
236
237 sage: from mjo.eja.eja_algebra import JordanSpinEJA
238
239 EXAMPLES:
240
241 The characteristic polynomial in the spin algebra is given in
242 Alizadeh, Example 11.11::
243
244 sage: J = JordanSpinEJA(3)
245 sage: p = J.characteristic_polynomial(); p
246 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
247 sage: xvec = J.one().vector()
248 sage: p(*xvec)
249 t^2 - 2*t + 1
250
251 """
252 r = self.rank()
253 n = self.dimension()
254
255 # The list of coefficient polynomials a_1, a_2, ..., a_n.
256 a = [ self._charpoly_coeff(i) for i in range(n) ]
257
258 # We go to a bit of trouble here to reorder the
259 # indeterminates, so that it's easier to evaluate the
260 # characteristic polynomial at x's coordinates and get back
261 # something in terms of t, which is what we want.
262 R = a[0].parent()
263 S = PolynomialRing(self.base_ring(),'t')
264 t = S.gen(0)
265 S = PolynomialRing(S, R.variable_names())
266 t = S(t)
267
268 # Note: all entries past the rth should be zero. The
269 # coefficient of the highest power (x^r) is 1, but it doesn't
270 # appear in the solution vector which contains coefficients
271 # for the other powers (to make them sum to x^r).
272 if (r < n):
273 a[r] = 1 # corresponds to x^r
274 else:
275 # When the rank is equal to the dimension, trying to
276 # assign a[r] goes out-of-bounds.
277 a.append(1) # corresponds to x^r
278
279 return sum( a[k]*(t**k) for k in range(len(a)) )
280
281
282 def inner_product(self, x, y):
283 """
284 The inner product associated with this Euclidean Jordan algebra.
285
286 Defaults to the trace inner product, but can be overridden by
287 subclasses if they are sure that the necessary properties are
288 satisfied.
289
290 SETUP::
291
292 sage: from mjo.eja.eja_algebra import random_eja
293
294 EXAMPLES:
295
296 The inner product must satisfy its axiom for this algebra to truly
297 be a Euclidean Jordan Algebra::
298
299 sage: set_random_seed()
300 sage: J = random_eja()
301 sage: x = J.random_element()
302 sage: y = J.random_element()
303 sage: z = J.random_element()
304 sage: (x*y).inner_product(z) == y.inner_product(x*z)
305 True
306
307 """
308 if (not x in self) or (not y in self):
309 raise TypeError("arguments must live in this algebra")
310 return x.trace_inner_product(y)
311
312
313 def natural_basis(self):
314 """
315 Return a more-natural representation of this algebra's basis.
316
317 Every finite-dimensional Euclidean Jordan Algebra is a direct
318 sum of five simple algebras, four of which comprise Hermitian
319 matrices. This method returns the original "natural" basis
320 for our underlying vector space. (Typically, the natural basis
321 is used to construct the multiplication table in the first place.)
322
323 Note that this will always return a matrix. The standard basis
324 in `R^n` will be returned as `n`-by-`1` column matrices.
325
326 SETUP::
327
328 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
329 ....: RealSymmetricEJA)
330
331 EXAMPLES::
332
333 sage: J = RealSymmetricEJA(2)
334 sage: J.basis()
335 Family (e0, e1, e2)
336 sage: J.natural_basis()
337 (
338 [1 0] [0 1] [0 0]
339 [0 0], [1 0], [0 1]
340 )
341
342 ::
343
344 sage: J = JordanSpinEJA(2)
345 sage: J.basis()
346 Family (e0, e1)
347 sage: J.natural_basis()
348 (
349 [1] [0]
350 [0], [1]
351 )
352
353 """
354 if self._natural_basis is None:
355 return tuple( b.vector().column() for b in self.basis() )
356 else:
357 return self._natural_basis
358
359
360 def rank(self):
361 """
362 Return the rank of this EJA.
363 """
364 if self._rank is None:
365 raise ValueError("no rank specified at genesis")
366 else:
367 return self._rank
368
369
370 def vector_space(self):
371 """
372 Return the vector space that underlies this algebra.
373
374 SETUP::
375
376 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
377
378 EXAMPLES::
379
380 sage: J = RealSymmetricEJA(2)
381 sage: J.vector_space()
382 Vector space of dimension 3 over Rational Field
383
384 """
385 return self.zero().vector().parent().ambient_vector_space()
386
387
388 class Element(FiniteDimensionalAlgebraElement):
389 """
390 An element of a Euclidean Jordan algebra.
391 """
392
393 def __dir__(self):
394 """
395 Oh man, I should not be doing this. This hides the "disabled"
396 methods ``left_matrix`` and ``matrix`` from introspection;
397 in particular it removes them from tab-completion.
398 """
399 return filter(lambda s: s not in ['left_matrix', 'matrix'],
400 dir(self.__class__) )
401
402
403 def __init__(self, A, elt=None):
404 """
405
406 SETUP::
407
408 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
409
410 EXAMPLES:
411
412 The identity in `S^n` is converted to the identity in the EJA::
413
414 sage: J = RealSymmetricEJA(3)
415 sage: I = matrix.identity(QQ,3)
416 sage: J(I) == J.one()
417 True
418
419 This skew-symmetric matrix can't be represented in the EJA::
420
421 sage: J = RealSymmetricEJA(3)
422 sage: A = matrix(QQ,3, lambda i,j: i-j)
423 sage: J(A)
424 Traceback (most recent call last):
425 ...
426 ArithmeticError: vector is not in free module
427
428 """
429 # Goal: if we're given a matrix, and if it lives in our
430 # parent algebra's "natural ambient space," convert it
431 # into an algebra element.
432 #
433 # The catch is, we make a recursive call after converting
434 # the given matrix into a vector that lives in the algebra.
435 # This we need to try the parent class initializer first,
436 # to avoid recursing forever if we're given something that
437 # already fits into the algebra, but also happens to live
438 # in the parent's "natural ambient space" (this happens with
439 # vectors in R^n).
440 try:
441 FiniteDimensionalAlgebraElement.__init__(self, A, elt)
442 except ValueError:
443 natural_basis = A.natural_basis()
444 if elt in natural_basis[0].matrix_space():
445 # Thanks for nothing! Matrix spaces aren't vector
446 # spaces in Sage, so we have to figure out its
447 # natural-basis coordinates ourselves.
448 V = VectorSpace(elt.base_ring(), elt.nrows()**2)
449 W = V.span( _mat2vec(s) for s in natural_basis )
450 coords = W.coordinates(_mat2vec(elt))
451 FiniteDimensionalAlgebraElement.__init__(self, A, coords)
452
453 def __pow__(self, n):
454 """
455 Return ``self`` raised to the power ``n``.
456
457 Jordan algebras are always power-associative; see for
458 example Faraut and Koranyi, Proposition II.1.2 (ii).
459
460 .. WARNING:
461
462 We have to override this because our superclass uses row vectors
463 instead of column vectors! We, on the other hand, assume column
464 vectors everywhere.
465
466 SETUP::
467
468 sage: from mjo.eja.eja_algebra import random_eja
469
470 EXAMPLES::
471
472 sage: set_random_seed()
473 sage: x = random_eja().random_element()
474 sage: x.operator()(x) == (x^2)
475 True
476
477 A few examples of power-associativity::
478
479 sage: set_random_seed()
480 sage: x = random_eja().random_element()
481 sage: x*(x*x)*(x*x) == x^5
482 True
483 sage: (x*x)*(x*x*x) == x^5
484 True
485
486 We also know that powers operator-commute (Koecher, Chapter
487 III, Corollary 1)::
488
489 sage: set_random_seed()
490 sage: x = random_eja().random_element()
491 sage: m = ZZ.random_element(0,10)
492 sage: n = ZZ.random_element(0,10)
493 sage: Lxm = (x^m).operator()
494 sage: Lxn = (x^n).operator()
495 sage: Lxm*Lxn == Lxn*Lxm
496 True
497
498 """
499 if n == 0:
500 return self.parent().one()
501 elif n == 1:
502 return self
503 else:
504 return (self.operator()**(n-1))(self)
505
506
507 def apply_univariate_polynomial(self, p):
508 """
509 Apply the univariate polynomial ``p`` to this element.
510
511 A priori, SageMath won't allow us to apply a univariate
512 polynomial to an element of an EJA, because we don't know
513 that EJAs are rings (they are usually not associative). Of
514 course, we know that EJAs are power-associative, so the
515 operation is ultimately kosher. This function sidesteps
516 the CAS to get the answer we want and expect.
517
518 SETUP::
519
520 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
521 ....: random_eja)
522
523 EXAMPLES::
524
525 sage: R = PolynomialRing(QQ, 't')
526 sage: t = R.gen(0)
527 sage: p = t^4 - t^3 + 5*t - 2
528 sage: J = RealCartesianProductEJA(5)
529 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
530 True
531
532 TESTS:
533
534 We should always get back an element of the algebra::
535
536 sage: set_random_seed()
537 sage: p = PolynomialRing(QQ, 't').random_element()
538 sage: J = random_eja()
539 sage: x = J.random_element()
540 sage: x.apply_univariate_polynomial(p) in J
541 True
542
543 """
544 if len(p.variables()) > 1:
545 raise ValueError("not a univariate polynomial")
546 P = self.parent()
547 R = P.base_ring()
548 # Convert the coeficcients to the parent's base ring,
549 # because a priori they might live in an (unnecessarily)
550 # larger ring for which P.sum() would fail below.
551 cs = [ R(c) for c in p.coefficients(sparse=False) ]
552 return P.sum( cs[k]*(self**k) for k in range(len(cs)) )
553
554
555 def characteristic_polynomial(self):
556 """
557 Return the characteristic polynomial of this element.
558
559 SETUP::
560
561 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
562
563 EXAMPLES:
564
565 The rank of `R^3` is three, and the minimal polynomial of
566 the identity element is `(t-1)` from which it follows that
567 the characteristic polynomial should be `(t-1)^3`::
568
569 sage: J = RealCartesianProductEJA(3)
570 sage: J.one().characteristic_polynomial()
571 t^3 - 3*t^2 + 3*t - 1
572
573 Likewise, the characteristic of the zero element in the
574 rank-three algebra `R^{n}` should be `t^{3}`::
575
576 sage: J = RealCartesianProductEJA(3)
577 sage: J.zero().characteristic_polynomial()
578 t^3
579
580 The characteristic polynomial of an element should evaluate
581 to zero on that element::
582
583 sage: set_random_seed()
584 sage: x = RealCartesianProductEJA(3).random_element()
585 sage: p = x.characteristic_polynomial()
586 sage: x.apply_univariate_polynomial(p)
587 0
588
589 """
590 p = self.parent().characteristic_polynomial()
591 return p(*self.vector())
592
593
594 def inner_product(self, other):
595 """
596 Return the parent algebra's inner product of myself and ``other``.
597
598 SETUP::
599
600 sage: from mjo.eja.eja_algebra import (
601 ....: ComplexHermitianEJA,
602 ....: JordanSpinEJA,
603 ....: QuaternionHermitianEJA,
604 ....: RealSymmetricEJA,
605 ....: random_eja)
606
607 EXAMPLES:
608
609 The inner product in the Jordan spin algebra is the usual
610 inner product on `R^n` (this example only works because the
611 basis for the Jordan algebra is the standard basis in `R^n`)::
612
613 sage: J = JordanSpinEJA(3)
614 sage: x = vector(QQ,[1,2,3])
615 sage: y = vector(QQ,[4,5,6])
616 sage: x.inner_product(y)
617 32
618 sage: J(x).inner_product(J(y))
619 32
620
621 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
622 multiplication is the usual matrix multiplication in `S^n`,
623 so the inner product of the identity matrix with itself
624 should be the `n`::
625
626 sage: J = RealSymmetricEJA(3)
627 sage: J.one().inner_product(J.one())
628 3
629
630 Likewise, the inner product on `C^n` is `<X,Y> =
631 Re(trace(X*Y))`, where we must necessarily take the real
632 part because the product of Hermitian matrices may not be
633 Hermitian::
634
635 sage: J = ComplexHermitianEJA(3)
636 sage: J.one().inner_product(J.one())
637 3
638
639 Ditto for the quaternions::
640
641 sage: J = QuaternionHermitianEJA(3)
642 sage: J.one().inner_product(J.one())
643 3
644
645 TESTS:
646
647 Ensure that we can always compute an inner product, and that
648 it gives us back a real number::
649
650 sage: set_random_seed()
651 sage: J = random_eja()
652 sage: x = J.random_element()
653 sage: y = J.random_element()
654 sage: x.inner_product(y) in RR
655 True
656
657 """
658 P = self.parent()
659 if not other in P:
660 raise TypeError("'other' must live in the same algebra")
661
662 return P.inner_product(self, other)
663
664
665 def operator_commutes_with(self, other):
666 """
667 Return whether or not this element operator-commutes
668 with ``other``.
669
670 SETUP::
671
672 sage: from mjo.eja.eja_algebra import random_eja
673
674 EXAMPLES:
675
676 The definition of a Jordan algebra says that any element
677 operator-commutes with its square::
678
679 sage: set_random_seed()
680 sage: x = random_eja().random_element()
681 sage: x.operator_commutes_with(x^2)
682 True
683
684 TESTS:
685
686 Test Lemma 1 from Chapter III of Koecher::
687
688 sage: set_random_seed()
689 sage: J = random_eja()
690 sage: u = J.random_element()
691 sage: v = J.random_element()
692 sage: lhs = u.operator_commutes_with(u*v)
693 sage: rhs = v.operator_commutes_with(u^2)
694 sage: lhs == rhs
695 True
696
697 Test the first polarization identity from my notes, Koecher Chapter
698 III, or from Baes (2.3)::
699
700 sage: set_random_seed()
701 sage: J = random_eja()
702 sage: x = J.random_element()
703 sage: y = J.random_element()
704 sage: Lx = x.operator()
705 sage: Ly = y.operator()
706 sage: Lxx = (x*x).operator()
707 sage: Lxy = (x*y).operator()
708 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
709 True
710
711 Test the second polarization identity from my notes or from
712 Baes (2.4)::
713
714 sage: set_random_seed()
715 sage: J = random_eja()
716 sage: x = J.random_element()
717 sage: y = J.random_element()
718 sage: z = J.random_element()
719 sage: Lx = x.operator()
720 sage: Ly = y.operator()
721 sage: Lz = z.operator()
722 sage: Lzy = (z*y).operator()
723 sage: Lxy = (x*y).operator()
724 sage: Lxz = (x*z).operator()
725 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
726 True
727
728 Test the third polarization identity from my notes or from
729 Baes (2.5)::
730
731 sage: set_random_seed()
732 sage: J = random_eja()
733 sage: u = J.random_element()
734 sage: y = J.random_element()
735 sage: z = J.random_element()
736 sage: Lu = u.operator()
737 sage: Ly = y.operator()
738 sage: Lz = z.operator()
739 sage: Lzy = (z*y).operator()
740 sage: Luy = (u*y).operator()
741 sage: Luz = (u*z).operator()
742 sage: Luyz = (u*(y*z)).operator()
743 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
744 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
745 sage: bool(lhs == rhs)
746 True
747
748 """
749 if not other in self.parent():
750 raise TypeError("'other' must live in the same algebra")
751
752 A = self.operator()
753 B = other.operator()
754 return (A*B == B*A)
755
756
757 def det(self):
758 """
759 Return my determinant, the product of my eigenvalues.
760
761 SETUP::
762
763 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
764 ....: random_eja)
765
766 EXAMPLES::
767
768 sage: J = JordanSpinEJA(2)
769 sage: e0,e1 = J.gens()
770 sage: x = sum( J.gens() )
771 sage: x.det()
772 0
773
774 ::
775
776 sage: J = JordanSpinEJA(3)
777 sage: e0,e1,e2 = J.gens()
778 sage: x = sum( J.gens() )
779 sage: x.det()
780 -1
781
782 TESTS:
783
784 An element is invertible if and only if its determinant is
785 non-zero::
786
787 sage: set_random_seed()
788 sage: x = random_eja().random_element()
789 sage: x.is_invertible() == (x.det() != 0)
790 True
791
792 """
793 P = self.parent()
794 r = P.rank()
795 p = P._charpoly_coeff(0)
796 # The _charpoly_coeff function already adds the factor of
797 # -1 to ensure that _charpoly_coeff(0) is really what
798 # appears in front of t^{0} in the charpoly. However,
799 # we want (-1)^r times THAT for the determinant.
800 return ((-1)**r)*p(*self.vector())
801
802
803 def inverse(self):
804 """
805 Return the Jordan-multiplicative inverse of this element.
806
807 ALGORITHM:
808
809 We appeal to the quadratic representation as in Koecher's
810 Theorem 12 in Chapter III, Section 5.
811
812 SETUP::
813
814 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
815 ....: random_eja)
816
817 EXAMPLES:
818
819 The inverse in the spin factor algebra is given in Alizadeh's
820 Example 11.11::
821
822 sage: set_random_seed()
823 sage: n = ZZ.random_element(1,10)
824 sage: J = JordanSpinEJA(n)
825 sage: x = J.random_element()
826 sage: while not x.is_invertible():
827 ....: x = J.random_element()
828 sage: x_vec = x.vector()
829 sage: x0 = x_vec[0]
830 sage: x_bar = x_vec[1:]
831 sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
832 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
833 sage: x_inverse = coeff*inv_vec
834 sage: x.inverse() == J(x_inverse)
835 True
836
837 TESTS:
838
839 The identity element is its own inverse::
840
841 sage: set_random_seed()
842 sage: J = random_eja()
843 sage: J.one().inverse() == J.one()
844 True
845
846 If an element has an inverse, it acts like one::
847
848 sage: set_random_seed()
849 sage: J = random_eja()
850 sage: x = J.random_element()
851 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
852 True
853
854 The inverse of the inverse is what we started with::
855
856 sage: set_random_seed()
857 sage: J = random_eja()
858 sage: x = J.random_element()
859 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
860 True
861
862 The zero element is never invertible::
863
864 sage: set_random_seed()
865 sage: J = random_eja().zero().inverse()
866 Traceback (most recent call last):
867 ...
868 ValueError: element is not invertible
869
870 """
871 if not self.is_invertible():
872 raise ValueError("element is not invertible")
873
874 return (~self.quadratic_representation())(self)
875
876
877 def is_invertible(self):
878 """
879 Return whether or not this element is invertible.
880
881 We can't use the superclass method because it relies on
882 the algebra being associative.
883
884 ALGORITHM:
885
886 The usual way to do this is to check if the determinant is
887 zero, but we need the characteristic polynomial for the
888 determinant. The minimal polynomial is a lot easier to get,
889 so we use Corollary 2 in Chapter V of Koecher to check
890 whether or not the paren't algebra's zero element is a root
891 of this element's minimal polynomial.
892
893 SETUP::
894
895 sage: from mjo.eja.eja_algebra import random_eja
896
897 TESTS:
898
899 The identity element is always invertible::
900
901 sage: set_random_seed()
902 sage: J = random_eja()
903 sage: J.one().is_invertible()
904 True
905
906 The zero element is never invertible::
907
908 sage: set_random_seed()
909 sage: J = random_eja()
910 sage: J.zero().is_invertible()
911 False
912
913 """
914 zero = self.parent().zero()
915 p = self.minimal_polynomial()
916 return not (p(zero) == zero)
917
918
919 def is_nilpotent(self):
920 """
921 Return whether or not some power of this element is zero.
922
923 The superclass method won't work unless we're in an
924 associative algebra, and we aren't. However, we generate
925 an assocoative subalgebra and we're nilpotent there if and
926 only if we're nilpotent here (probably).
927
928 SETUP::
929
930 sage: from mjo.eja.eja_algebra import random_eja
931
932 TESTS:
933
934 The identity element is never nilpotent::
935
936 sage: set_random_seed()
937 sage: random_eja().one().is_nilpotent()
938 False
939
940 The additive identity is always nilpotent::
941
942 sage: set_random_seed()
943 sage: random_eja().zero().is_nilpotent()
944 True
945
946 """
947 # The element we're going to call "is_nilpotent()" on.
948 # Either myself, interpreted as an element of a finite-
949 # dimensional algebra, or an element of an associative
950 # subalgebra.
951 elt = None
952
953 if self.parent().is_associative():
954 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
955 else:
956 V = self.span_of_powers()
957 assoc_subalg = self.subalgebra_generated_by()
958 # Mis-design warning: the basis used for span_of_powers()
959 # and subalgebra_generated_by() must be the same, and in
960 # the same order!
961 elt = assoc_subalg(V.coordinates(self.vector()))
962
963 # Recursive call, but should work since elt lives in an
964 # associative algebra.
965 return elt.is_nilpotent()
966
967
968 def is_regular(self):
969 """
970 Return whether or not this is a regular element.
971
972 SETUP::
973
974 sage: from mjo.eja.eja_algebra import JordanSpinEJA
975
976 EXAMPLES:
977
978 The identity element always has degree one, but any element
979 linearly-independent from it is regular::
980
981 sage: J = JordanSpinEJA(5)
982 sage: J.one().is_regular()
983 False
984 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
985 sage: for x in J.gens():
986 ....: (J.one() + x).is_regular()
987 False
988 True
989 True
990 True
991 True
992
993 """
994 return self.degree() == self.parent().rank()
995
996
997 def degree(self):
998 """
999 Compute the degree of this element the straightforward way
1000 according to the definition; by appending powers to a list
1001 and figuring out its dimension (that is, whether or not
1002 they're linearly dependent).
1003
1004 SETUP::
1005
1006 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1007
1008 EXAMPLES::
1009
1010 sage: J = JordanSpinEJA(4)
1011 sage: J.one().degree()
1012 1
1013 sage: e0,e1,e2,e3 = J.gens()
1014 sage: (e0 - e1).degree()
1015 2
1016
1017 In the spin factor algebra (of rank two), all elements that
1018 aren't multiples of the identity are regular::
1019
1020 sage: set_random_seed()
1021 sage: n = ZZ.random_element(1,10)
1022 sage: J = JordanSpinEJA(n)
1023 sage: x = J.random_element()
1024 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
1025 True
1026
1027 """
1028 return self.span_of_powers().dimension()
1029
1030
1031 def left_matrix(self):
1032 """
1033 Our parent class defines ``left_matrix`` and ``matrix``
1034 methods whose names are misleading. We don't want them.
1035 """
1036 raise NotImplementedError("use operator().matrix() instead")
1037
1038 matrix = left_matrix
1039
1040
1041 def minimal_polynomial(self):
1042 """
1043 Return the minimal polynomial of this element,
1044 as a function of the variable `t`.
1045
1046 ALGORITHM:
1047
1048 We restrict ourselves to the associative subalgebra
1049 generated by this element, and then return the minimal
1050 polynomial of this element's operator matrix (in that
1051 subalgebra). This works by Baes Proposition 2.3.16.
1052
1053 SETUP::
1054
1055 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1056 ....: random_eja)
1057
1058 TESTS:
1059
1060 The minimal polynomial of the identity and zero elements are
1061 always the same::
1062
1063 sage: set_random_seed()
1064 sage: J = random_eja()
1065 sage: J.one().minimal_polynomial()
1066 t - 1
1067 sage: J.zero().minimal_polynomial()
1068 t
1069
1070 The degree of an element is (by one definition) the degree
1071 of its minimal polynomial::
1072
1073 sage: set_random_seed()
1074 sage: x = random_eja().random_element()
1075 sage: x.degree() == x.minimal_polynomial().degree()
1076 True
1077
1078 The minimal polynomial and the characteristic polynomial coincide
1079 and are known (see Alizadeh, Example 11.11) for all elements of
1080 the spin factor algebra that aren't scalar multiples of the
1081 identity::
1082
1083 sage: set_random_seed()
1084 sage: n = ZZ.random_element(2,10)
1085 sage: J = JordanSpinEJA(n)
1086 sage: y = J.random_element()
1087 sage: while y == y.coefficient(0)*J.one():
1088 ....: y = J.random_element()
1089 sage: y0 = y.vector()[0]
1090 sage: y_bar = y.vector()[1:]
1091 sage: actual = y.minimal_polynomial()
1092 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1093 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1094 sage: bool(actual == expected)
1095 True
1096
1097 The minimal polynomial should always kill its element::
1098
1099 sage: set_random_seed()
1100 sage: x = random_eja().random_element()
1101 sage: p = x.minimal_polynomial()
1102 sage: x.apply_univariate_polynomial(p)
1103 0
1104
1105 """
1106 V = self.span_of_powers()
1107 assoc_subalg = self.subalgebra_generated_by()
1108 # Mis-design warning: the basis used for span_of_powers()
1109 # and subalgebra_generated_by() must be the same, and in
1110 # the same order!
1111 elt = assoc_subalg(V.coordinates(self.vector()))
1112 return elt.operator().minimal_polynomial()
1113
1114
1115
1116 def natural_representation(self):
1117 """
1118 Return a more-natural representation of this element.
1119
1120 Every finite-dimensional Euclidean Jordan Algebra is a
1121 direct sum of five simple algebras, four of which comprise
1122 Hermitian matrices. This method returns the original
1123 "natural" representation of this element as a Hermitian
1124 matrix, if it has one. If not, you get the usual representation.
1125
1126 SETUP::
1127
1128 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1129 ....: QuaternionHermitianEJA)
1130
1131 EXAMPLES::
1132
1133 sage: J = ComplexHermitianEJA(3)
1134 sage: J.one()
1135 e0 + e5 + e8
1136 sage: J.one().natural_representation()
1137 [1 0 0 0 0 0]
1138 [0 1 0 0 0 0]
1139 [0 0 1 0 0 0]
1140 [0 0 0 1 0 0]
1141 [0 0 0 0 1 0]
1142 [0 0 0 0 0 1]
1143
1144 ::
1145
1146 sage: J = QuaternionHermitianEJA(3)
1147 sage: J.one()
1148 e0 + e9 + e14
1149 sage: J.one().natural_representation()
1150 [1 0 0 0 0 0 0 0 0 0 0 0]
1151 [0 1 0 0 0 0 0 0 0 0 0 0]
1152 [0 0 1 0 0 0 0 0 0 0 0 0]
1153 [0 0 0 1 0 0 0 0 0 0 0 0]
1154 [0 0 0 0 1 0 0 0 0 0 0 0]
1155 [0 0 0 0 0 1 0 0 0 0 0 0]
1156 [0 0 0 0 0 0 1 0 0 0 0 0]
1157 [0 0 0 0 0 0 0 1 0 0 0 0]
1158 [0 0 0 0 0 0 0 0 1 0 0 0]
1159 [0 0 0 0 0 0 0 0 0 1 0 0]
1160 [0 0 0 0 0 0 0 0 0 0 1 0]
1161 [0 0 0 0 0 0 0 0 0 0 0 1]
1162
1163 """
1164 B = self.parent().natural_basis()
1165 W = B[0].matrix_space()
1166 return W.linear_combination(zip(self.vector(), B))
1167
1168
1169 def operator(self):
1170 """
1171 Return the left-multiplication-by-this-element
1172 operator on the ambient algebra.
1173
1174 SETUP::
1175
1176 sage: from mjo.eja.eja_algebra import random_eja
1177
1178 TESTS::
1179
1180 sage: set_random_seed()
1181 sage: J = random_eja()
1182 sage: x = J.random_element()
1183 sage: y = J.random_element()
1184 sage: x.operator()(y) == x*y
1185 True
1186 sage: y.operator()(x) == x*y
1187 True
1188
1189 """
1190 P = self.parent()
1191 fda_elt = FiniteDimensionalAlgebraElement(P, self)
1192 return FiniteDimensionalEuclideanJordanAlgebraOperator(
1193 P,
1194 P,
1195 fda_elt.matrix().transpose() )
1196
1197
1198 def quadratic_representation(self, other=None):
1199 """
1200 Return the quadratic representation of this element.
1201
1202 SETUP::
1203
1204 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1205 ....: random_eja)
1206
1207 EXAMPLES:
1208
1209 The explicit form in the spin factor algebra is given by
1210 Alizadeh's Example 11.12::
1211
1212 sage: set_random_seed()
1213 sage: n = ZZ.random_element(1,10)
1214 sage: J = JordanSpinEJA(n)
1215 sage: x = J.random_element()
1216 sage: x_vec = x.vector()
1217 sage: x0 = x_vec[0]
1218 sage: x_bar = x_vec[1:]
1219 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1220 sage: B = 2*x0*x_bar.row()
1221 sage: C = 2*x0*x_bar.column()
1222 sage: D = matrix.identity(QQ, n-1)
1223 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1224 sage: D = D + 2*x_bar.tensor_product(x_bar)
1225 sage: Q = matrix.block(2,2,[A,B,C,D])
1226 sage: Q == x.quadratic_representation().matrix()
1227 True
1228
1229 Test all of the properties from Theorem 11.2 in Alizadeh::
1230
1231 sage: set_random_seed()
1232 sage: J = random_eja()
1233 sage: x = J.random_element()
1234 sage: y = J.random_element()
1235 sage: Lx = x.operator()
1236 sage: Lxx = (x*x).operator()
1237 sage: Qx = x.quadratic_representation()
1238 sage: Qy = y.quadratic_representation()
1239 sage: Qxy = x.quadratic_representation(y)
1240 sage: Qex = J.one().quadratic_representation(x)
1241 sage: n = ZZ.random_element(10)
1242 sage: Qxn = (x^n).quadratic_representation()
1243
1244 Property 1:
1245
1246 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1247 True
1248
1249 Property 2 (multiply on the right for :trac:`28272`):
1250
1251 sage: alpha = QQ.random_element()
1252 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1253 True
1254
1255 Property 3:
1256
1257 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1258 True
1259
1260 sage: not x.is_invertible() or (
1261 ....: ~Qx
1262 ....: ==
1263 ....: x.inverse().quadratic_representation() )
1264 True
1265
1266 sage: Qxy(J.one()) == x*y
1267 True
1268
1269 Property 4:
1270
1271 sage: not x.is_invertible() or (
1272 ....: x.quadratic_representation(x.inverse())*Qx
1273 ....: == Qx*x.quadratic_representation(x.inverse()) )
1274 True
1275
1276 sage: not x.is_invertible() or (
1277 ....: x.quadratic_representation(x.inverse())*Qx
1278 ....: ==
1279 ....: 2*x.operator()*Qex - Qx )
1280 True
1281
1282 sage: 2*x.operator()*Qex - Qx == Lxx
1283 True
1284
1285 Property 5:
1286
1287 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1288 True
1289
1290 Property 6:
1291
1292 sage: Qxn == (Qx)^n
1293 True
1294
1295 Property 7:
1296
1297 sage: not x.is_invertible() or (
1298 ....: Qx*x.inverse().operator() == Lx )
1299 True
1300
1301 Property 8:
1302
1303 sage: not x.operator_commutes_with(y) or (
1304 ....: Qx(y)^n == Qxn(y^n) )
1305 True
1306
1307 """
1308 if other is None:
1309 other=self
1310 elif not other in self.parent():
1311 raise TypeError("'other' must live in the same algebra")
1312
1313 L = self.operator()
1314 M = other.operator()
1315 return ( L*M + M*L - (self*other).operator() )
1316
1317
1318 def span_of_powers(self):
1319 """
1320 Return the vector space spanned by successive powers of
1321 this element.
1322 """
1323 # The dimension of the subalgebra can't be greater than
1324 # the big algebra, so just put everything into a list
1325 # and let span() get rid of the excess.
1326 #
1327 # We do the extra ambient_vector_space() in case we're messing
1328 # with polynomials and the direct parent is a module.
1329 V = self.parent().vector_space()
1330 return V.span( (self**d).vector() for d in xrange(V.dimension()) )
1331
1332
1333 def subalgebra_generated_by(self):
1334 """
1335 Return the associative subalgebra of the parent EJA generated
1336 by this element.
1337
1338 SETUP::
1339
1340 sage: from mjo.eja.eja_algebra import random_eja
1341
1342 TESTS::
1343
1344 sage: set_random_seed()
1345 sage: x = random_eja().random_element()
1346 sage: x.subalgebra_generated_by().is_associative()
1347 True
1348
1349 Squaring in the subalgebra should work the same as in
1350 the superalgebra::
1351
1352 sage: set_random_seed()
1353 sage: x = random_eja().random_element()
1354 sage: u = x.subalgebra_generated_by().random_element()
1355 sage: u.operator()(u) == u^2
1356 True
1357
1358 """
1359 # First get the subspace spanned by the powers of myself...
1360 V = self.span_of_powers()
1361 F = self.base_ring()
1362
1363 # Now figure out the entries of the right-multiplication
1364 # matrix for the successive basis elements b0, b1,... of
1365 # that subspace.
1366 mats = []
1367 for b_right in V.basis():
1368 eja_b_right = self.parent()(b_right)
1369 b_right_rows = []
1370 # The first row of the right-multiplication matrix by
1371 # b1 is what we get if we apply that matrix to b1. The
1372 # second row of the right multiplication matrix by b1
1373 # is what we get when we apply that matrix to b2...
1374 #
1375 # IMPORTANT: this assumes that all vectors are COLUMN
1376 # vectors, unlike our superclass (which uses row vectors).
1377 for b_left in V.basis():
1378 eja_b_left = self.parent()(b_left)
1379 # Multiply in the original EJA, but then get the
1380 # coordinates from the subalgebra in terms of its
1381 # basis.
1382 this_row = V.coordinates((eja_b_left*eja_b_right).vector())
1383 b_right_rows.append(this_row)
1384 b_right_matrix = matrix(F, b_right_rows)
1385 mats.append(b_right_matrix)
1386
1387 # It's an algebra of polynomials in one element, and EJAs
1388 # are power-associative.
1389 #
1390 # TODO: choose generator names intelligently.
1391 return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
1392
1393
1394 def subalgebra_idempotent(self):
1395 """
1396 Find an idempotent in the associative subalgebra I generate
1397 using Proposition 2.3.5 in Baes.
1398
1399 SETUP::
1400
1401 sage: from mjo.eja.eja_algebra import random_eja
1402
1403 TESTS::
1404
1405 sage: set_random_seed()
1406 sage: J = random_eja()
1407 sage: x = J.random_element()
1408 sage: while x.is_nilpotent():
1409 ....: x = J.random_element()
1410 sage: c = x.subalgebra_idempotent()
1411 sage: c^2 == c
1412 True
1413
1414 """
1415 if self.is_nilpotent():
1416 raise ValueError("this only works with non-nilpotent elements!")
1417
1418 V = self.span_of_powers()
1419 J = self.subalgebra_generated_by()
1420 # Mis-design warning: the basis used for span_of_powers()
1421 # and subalgebra_generated_by() must be the same, and in
1422 # the same order!
1423 u = J(V.coordinates(self.vector()))
1424
1425 # The image of the matrix of left-u^m-multiplication
1426 # will be minimal for some natural number s...
1427 s = 0
1428 minimal_dim = V.dimension()
1429 for i in xrange(1, V.dimension()):
1430 this_dim = (u**i).operator().matrix().image().dimension()
1431 if this_dim < minimal_dim:
1432 minimal_dim = this_dim
1433 s = i
1434
1435 # Now minimal_matrix should correspond to the smallest
1436 # non-zero subspace in Baes's (or really, Koecher's)
1437 # proposition.
1438 #
1439 # However, we need to restrict the matrix to work on the
1440 # subspace... or do we? Can't we just solve, knowing that
1441 # A(c) = u^(s+1) should have a solution in the big space,
1442 # too?
1443 #
1444 # Beware, solve_right() means that we're using COLUMN vectors.
1445 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1446 u_next = u**(s+1)
1447 A = u_next.operator().matrix()
1448 c_coordinates = A.solve_right(u_next.vector())
1449
1450 # Now c_coordinates is the idempotent we want, but it's in
1451 # the coordinate system of the subalgebra.
1452 #
1453 # We need the basis for J, but as elements of the parent algebra.
1454 #
1455 basis = [self.parent(v) for v in V.basis()]
1456 return self.parent().linear_combination(zip(c_coordinates, basis))
1457
1458
1459 def trace(self):
1460 """
1461 Return my trace, the sum of my eigenvalues.
1462
1463 SETUP::
1464
1465 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1466 ....: RealCartesianProductEJA,
1467 ....: random_eja)
1468
1469 EXAMPLES::
1470
1471 sage: J = JordanSpinEJA(3)
1472 sage: x = sum(J.gens())
1473 sage: x.trace()
1474 2
1475
1476 ::
1477
1478 sage: J = RealCartesianProductEJA(5)
1479 sage: J.one().trace()
1480 5
1481
1482 TESTS:
1483
1484 The trace of an element is a real number::
1485
1486 sage: set_random_seed()
1487 sage: J = random_eja()
1488 sage: J.random_element().trace() in J.base_ring()
1489 True
1490
1491 """
1492 P = self.parent()
1493 r = P.rank()
1494 p = P._charpoly_coeff(r-1)
1495 # The _charpoly_coeff function already adds the factor of
1496 # -1 to ensure that _charpoly_coeff(r-1) is really what
1497 # appears in front of t^{r-1} in the charpoly. However,
1498 # we want the negative of THAT for the trace.
1499 return -p(*self.vector())
1500
1501
1502 def trace_inner_product(self, other):
1503 """
1504 Return the trace inner product of myself and ``other``.
1505
1506 SETUP::
1507
1508 sage: from mjo.eja.eja_algebra import random_eja
1509
1510 TESTS:
1511
1512 The trace inner product is commutative::
1513
1514 sage: set_random_seed()
1515 sage: J = random_eja()
1516 sage: x = J.random_element(); y = J.random_element()
1517 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1518 True
1519
1520 The trace inner product is bilinear::
1521
1522 sage: set_random_seed()
1523 sage: J = random_eja()
1524 sage: x = J.random_element()
1525 sage: y = J.random_element()
1526 sage: z = J.random_element()
1527 sage: a = QQ.random_element();
1528 sage: actual = (a*(x+z)).trace_inner_product(y)
1529 sage: expected = ( a*x.trace_inner_product(y) +
1530 ....: a*z.trace_inner_product(y) )
1531 sage: actual == expected
1532 True
1533 sage: actual = x.trace_inner_product(a*(y+z))
1534 sage: expected = ( a*x.trace_inner_product(y) +
1535 ....: a*x.trace_inner_product(z) )
1536 sage: actual == expected
1537 True
1538
1539 The trace inner product satisfies the compatibility
1540 condition in the definition of a Euclidean Jordan algebra::
1541
1542 sage: set_random_seed()
1543 sage: J = random_eja()
1544 sage: x = J.random_element()
1545 sage: y = J.random_element()
1546 sage: z = J.random_element()
1547 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1548 True
1549
1550 """
1551 if not other in self.parent():
1552 raise TypeError("'other' must live in the same algebra")
1553
1554 return (self*other).trace()
1555
1556
1557 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
1558 """
1559 Return the Euclidean Jordan Algebra corresponding to the set
1560 `R^n` under the Hadamard product.
1561
1562 Note: this is nothing more than the Cartesian product of ``n``
1563 copies of the spin algebra. Once Cartesian product algebras
1564 are implemented, this can go.
1565
1566 SETUP::
1567
1568 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
1569
1570 EXAMPLES:
1571
1572 This multiplication table can be verified by hand::
1573
1574 sage: J = RealCartesianProductEJA(3)
1575 sage: e0,e1,e2 = J.gens()
1576 sage: e0*e0
1577 e0
1578 sage: e0*e1
1579 0
1580 sage: e0*e2
1581 0
1582 sage: e1*e1
1583 e1
1584 sage: e1*e2
1585 0
1586 sage: e2*e2
1587 e2
1588
1589 """
1590 @staticmethod
1591 def __classcall_private__(cls, n, field=QQ):
1592 # The FiniteDimensionalAlgebra constructor takes a list of
1593 # matrices, the ith representing right multiplication by the ith
1594 # basis element in the vector space. So if e_1 = (1,0,0), then
1595 # right (Hadamard) multiplication of x by e_1 picks out the first
1596 # component of x; and likewise for the ith basis element e_i.
1597 Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i))
1598 for i in xrange(n) ]
1599
1600 fdeja = super(RealCartesianProductEJA, cls)
1601 return fdeja.__classcall_private__(cls, field, Qs, rank=n)
1602
1603 def inner_product(self, x, y):
1604 return _usual_ip(x,y)
1605
1606
1607 def random_eja():
1608 """
1609 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1610
1611 ALGORITHM:
1612
1613 For now, we choose a random natural number ``n`` (greater than zero)
1614 and then give you back one of the following:
1615
1616 * The cartesian product of the rational numbers ``n`` times; this is
1617 ``QQ^n`` with the Hadamard product.
1618
1619 * The Jordan spin algebra on ``QQ^n``.
1620
1621 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1622 product.
1623
1624 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1625 in the space of ``2n``-by-``2n`` real symmetric matrices.
1626
1627 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1628 in the space of ``4n``-by-``4n`` real symmetric matrices.
1629
1630 Later this might be extended to return Cartesian products of the
1631 EJAs above.
1632
1633 SETUP::
1634
1635 sage: from mjo.eja.eja_algebra import random_eja
1636
1637 TESTS::
1638
1639 sage: random_eja()
1640 Euclidean Jordan algebra of degree...
1641
1642 """
1643
1644 # The max_n component lets us choose different upper bounds on the
1645 # value "n" that gets passed to the constructor. This is needed
1646 # because e.g. R^{10} is reasonable to test, while the Hermitian
1647 # 10-by-10 quaternion matrices are not.
1648 (constructor, max_n) = choice([(RealCartesianProductEJA, 6),
1649 (JordanSpinEJA, 6),
1650 (RealSymmetricEJA, 5),
1651 (ComplexHermitianEJA, 4),
1652 (QuaternionHermitianEJA, 3)])
1653 n = ZZ.random_element(1, max_n)
1654 return constructor(n, field=QQ)
1655
1656
1657
1658 def _real_symmetric_basis(n, field=QQ):
1659 """
1660 Return a basis for the space of real symmetric n-by-n matrices.
1661 """
1662 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1663 # coordinates.
1664 S = []
1665 for i in xrange(n):
1666 for j in xrange(i+1):
1667 Eij = matrix(field, n, lambda k,l: k==i and l==j)
1668 if i == j:
1669 Sij = Eij
1670 else:
1671 # Beware, orthogonal but not normalized!
1672 Sij = Eij + Eij.transpose()
1673 S.append(Sij)
1674 return tuple(S)
1675
1676
1677 def _complex_hermitian_basis(n, field=QQ):
1678 """
1679 Returns a basis for the space of complex Hermitian n-by-n matrices.
1680
1681 SETUP::
1682
1683 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
1684
1685 TESTS::
1686
1687 sage: set_random_seed()
1688 sage: n = ZZ.random_element(1,5)
1689 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1690 True
1691
1692 """
1693 F = QuadraticField(-1, 'I')
1694 I = F.gen()
1695
1696 # This is like the symmetric case, but we need to be careful:
1697 #
1698 # * We want conjugate-symmetry, not just symmetry.
1699 # * The diagonal will (as a result) be real.
1700 #
1701 S = []
1702 for i in xrange(n):
1703 for j in xrange(i+1):
1704 Eij = matrix(field, n, lambda k,l: k==i and l==j)
1705 if i == j:
1706 Sij = _embed_complex_matrix(Eij)
1707 S.append(Sij)
1708 else:
1709 # Beware, orthogonal but not normalized! The second one
1710 # has a minus because it's conjugated.
1711 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1712 S.append(Sij_real)
1713 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1714 S.append(Sij_imag)
1715 return tuple(S)
1716
1717
1718 def _quaternion_hermitian_basis(n, field=QQ):
1719 """
1720 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1721
1722 SETUP::
1723
1724 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
1725
1726 TESTS::
1727
1728 sage: set_random_seed()
1729 sage: n = ZZ.random_element(1,5)
1730 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1731 True
1732
1733 """
1734 Q = QuaternionAlgebra(QQ,-1,-1)
1735 I,J,K = Q.gens()
1736
1737 # This is like the symmetric case, but we need to be careful:
1738 #
1739 # * We want conjugate-symmetry, not just symmetry.
1740 # * The diagonal will (as a result) be real.
1741 #
1742 S = []
1743 for i in xrange(n):
1744 for j in xrange(i+1):
1745 Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1746 if i == j:
1747 Sij = _embed_quaternion_matrix(Eij)
1748 S.append(Sij)
1749 else:
1750 # Beware, orthogonal but not normalized! The second,
1751 # third, and fourth ones have a minus because they're
1752 # conjugated.
1753 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1754 S.append(Sij_real)
1755 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1756 S.append(Sij_I)
1757 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1758 S.append(Sij_J)
1759 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1760 S.append(Sij_K)
1761 return tuple(S)
1762
1763
1764 def _mat2vec(m):
1765 return vector(m.base_ring(), m.list())
1766
1767 def _vec2mat(v):
1768 return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1769
1770 def _multiplication_table_from_matrix_basis(basis):
1771 """
1772 At least three of the five simple Euclidean Jordan algebras have the
1773 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1774 multiplication on the right is matrix multiplication. Given a basis
1775 for the underlying matrix space, this function returns a
1776 multiplication table (obtained by looping through the basis
1777 elements) for an algebra of those matrices. A reordered copy
1778 of the basis is also returned to work around the fact that
1779 the ``span()`` in this function will change the order of the basis
1780 from what we think it is, to... something else.
1781 """
1782 # In S^2, for example, we nominally have four coordinates even
1783 # though the space is of dimension three only. The vector space V
1784 # is supposed to hold the entire long vector, and the subspace W
1785 # of V will be spanned by the vectors that arise from symmetric
1786 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1787 field = basis[0].base_ring()
1788 dimension = basis[0].nrows()
1789
1790 V = VectorSpace(field, dimension**2)
1791 W = V.span( _mat2vec(s) for s in basis )
1792
1793 # Taking the span above reorders our basis (thanks, jerk!) so we
1794 # need to put our "matrix basis" in the same order as the
1795 # (reordered) vector basis.
1796 S = tuple( _vec2mat(b) for b in W.basis() )
1797
1798 Qs = []
1799 for s in S:
1800 # Brute force the multiplication-by-s matrix by looping
1801 # through all elements of the basis and doing the computation
1802 # to find out what the corresponding row should be. BEWARE:
1803 # these multiplication tables won't be symmetric! It therefore
1804 # becomes REALLY IMPORTANT that the underlying algebra
1805 # constructor uses ROW vectors and not COLUMN vectors. That's
1806 # why we're computing rows here and not columns.
1807 Q_rows = []
1808 for t in S:
1809 this_row = _mat2vec((s*t + t*s)/2)
1810 Q_rows.append(W.coordinates(this_row))
1811 Q = matrix(field, W.dimension(), Q_rows)
1812 Qs.append(Q)
1813
1814 return (Qs, S)
1815
1816
1817 def _embed_complex_matrix(M):
1818 """
1819 Embed the n-by-n complex matrix ``M`` into the space of real
1820 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1821 bi` to the block matrix ``[[a,b],[-b,a]]``.
1822
1823 SETUP::
1824
1825 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
1826
1827 EXAMPLES::
1828
1829 sage: F = QuadraticField(-1,'i')
1830 sage: x1 = F(4 - 2*i)
1831 sage: x2 = F(1 + 2*i)
1832 sage: x3 = F(-i)
1833 sage: x4 = F(6)
1834 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1835 sage: _embed_complex_matrix(M)
1836 [ 4 -2| 1 2]
1837 [ 2 4|-2 1]
1838 [-----+-----]
1839 [ 0 -1| 6 0]
1840 [ 1 0| 0 6]
1841
1842 TESTS:
1843
1844 Embedding is a homomorphism (isomorphism, in fact)::
1845
1846 sage: set_random_seed()
1847 sage: n = ZZ.random_element(5)
1848 sage: F = QuadraticField(-1, 'i')
1849 sage: X = random_matrix(F, n)
1850 sage: Y = random_matrix(F, n)
1851 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1852 sage: expected = _embed_complex_matrix(X*Y)
1853 sage: actual == expected
1854 True
1855
1856 """
1857 n = M.nrows()
1858 if M.ncols() != n:
1859 raise ValueError("the matrix 'M' must be square")
1860 field = M.base_ring()
1861 blocks = []
1862 for z in M.list():
1863 a = z.real()
1864 b = z.imag()
1865 blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1866
1867 # We can drop the imaginaries here.
1868 return matrix.block(field.base_ring(), n, blocks)
1869
1870
1871 def _unembed_complex_matrix(M):
1872 """
1873 The inverse of _embed_complex_matrix().
1874
1875 SETUP::
1876
1877 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
1878 ....: _unembed_complex_matrix)
1879
1880 EXAMPLES::
1881
1882 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1883 ....: [-2, 1, -4, 3],
1884 ....: [ 9, 10, 11, 12],
1885 ....: [-10, 9, -12, 11] ])
1886 sage: _unembed_complex_matrix(A)
1887 [ 2*i + 1 4*i + 3]
1888 [ 10*i + 9 12*i + 11]
1889
1890 TESTS:
1891
1892 Unembedding is the inverse of embedding::
1893
1894 sage: set_random_seed()
1895 sage: F = QuadraticField(-1, 'i')
1896 sage: M = random_matrix(F, 3)
1897 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1898 True
1899
1900 """
1901 n = ZZ(M.nrows())
1902 if M.ncols() != n:
1903 raise ValueError("the matrix 'M' must be square")
1904 if not n.mod(2).is_zero():
1905 raise ValueError("the matrix 'M' must be a complex embedding")
1906
1907 F = QuadraticField(-1, 'i')
1908 i = F.gen()
1909
1910 # Go top-left to bottom-right (reading order), converting every
1911 # 2-by-2 block we see to a single complex element.
1912 elements = []
1913 for k in xrange(n/2):
1914 for j in xrange(n/2):
1915 submat = M[2*k:2*k+2,2*j:2*j+2]
1916 if submat[0,0] != submat[1,1]:
1917 raise ValueError('bad on-diagonal submatrix')
1918 if submat[0,1] != -submat[1,0]:
1919 raise ValueError('bad off-diagonal submatrix')
1920 z = submat[0,0] + submat[0,1]*i
1921 elements.append(z)
1922
1923 return matrix(F, n/2, elements)
1924
1925
1926 def _embed_quaternion_matrix(M):
1927 """
1928 Embed the n-by-n quaternion matrix ``M`` into the space of real
1929 matrices of size 4n-by-4n by first sending each quaternion entry
1930 `z = a + bi + cj + dk` to the block-complex matrix
1931 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1932 a real matrix.
1933
1934 SETUP::
1935
1936 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
1937
1938 EXAMPLES::
1939
1940 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1941 sage: i,j,k = Q.gens()
1942 sage: x = 1 + 2*i + 3*j + 4*k
1943 sage: M = matrix(Q, 1, [[x]])
1944 sage: _embed_quaternion_matrix(M)
1945 [ 1 2 3 4]
1946 [-2 1 -4 3]
1947 [-3 4 1 -2]
1948 [-4 -3 2 1]
1949
1950 Embedding is a homomorphism (isomorphism, in fact)::
1951
1952 sage: set_random_seed()
1953 sage: n = ZZ.random_element(5)
1954 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1955 sage: X = random_matrix(Q, n)
1956 sage: Y = random_matrix(Q, n)
1957 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1958 sage: expected = _embed_quaternion_matrix(X*Y)
1959 sage: actual == expected
1960 True
1961
1962 """
1963 quaternions = M.base_ring()
1964 n = M.nrows()
1965 if M.ncols() != n:
1966 raise ValueError("the matrix 'M' must be square")
1967
1968 F = QuadraticField(-1, 'i')
1969 i = F.gen()
1970
1971 blocks = []
1972 for z in M.list():
1973 t = z.coefficient_tuple()
1974 a = t[0]
1975 b = t[1]
1976 c = t[2]
1977 d = t[3]
1978 cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1979 [-c + d*i, a - b*i]])
1980 blocks.append(_embed_complex_matrix(cplx_matrix))
1981
1982 # We should have real entries by now, so use the realest field
1983 # we've got for the return value.
1984 return matrix.block(quaternions.base_ring(), n, blocks)
1985
1986
1987 def _unembed_quaternion_matrix(M):
1988 """
1989 The inverse of _embed_quaternion_matrix().
1990
1991 SETUP::
1992
1993 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1994 ....: _unembed_quaternion_matrix)
1995
1996 EXAMPLES::
1997
1998 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1999 ....: [-2, 1, -4, 3],
2000 ....: [-3, 4, 1, -2],
2001 ....: [-4, -3, 2, 1]])
2002 sage: _unembed_quaternion_matrix(M)
2003 [1 + 2*i + 3*j + 4*k]
2004
2005 TESTS:
2006
2007 Unembedding is the inverse of embedding::
2008
2009 sage: set_random_seed()
2010 sage: Q = QuaternionAlgebra(QQ, -1, -1)
2011 sage: M = random_matrix(Q, 3)
2012 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
2013 True
2014
2015 """
2016 n = ZZ(M.nrows())
2017 if M.ncols() != n:
2018 raise ValueError("the matrix 'M' must be square")
2019 if not n.mod(4).is_zero():
2020 raise ValueError("the matrix 'M' must be a complex embedding")
2021
2022 Q = QuaternionAlgebra(QQ,-1,-1)
2023 i,j,k = Q.gens()
2024
2025 # Go top-left to bottom-right (reading order), converting every
2026 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2027 # quaternion block.
2028 elements = []
2029 for l in xrange(n/4):
2030 for m in xrange(n/4):
2031 submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
2032 if submat[0,0] != submat[1,1].conjugate():
2033 raise ValueError('bad on-diagonal submatrix')
2034 if submat[0,1] != -submat[1,0].conjugate():
2035 raise ValueError('bad off-diagonal submatrix')
2036 z = submat[0,0].real() + submat[0,0].imag()*i
2037 z += submat[0,1].real()*j + submat[0,1].imag()*k
2038 elements.append(z)
2039
2040 return matrix(Q, n/4, elements)
2041
2042
2043 # The usual inner product on R^n.
2044 def _usual_ip(x,y):
2045 return x.vector().inner_product(y.vector())
2046
2047 # The inner product used for the real symmetric simple EJA.
2048 # We keep it as a separate function because e.g. the complex
2049 # algebra uses the same inner product, except divided by 2.
2050 def _matrix_ip(X,Y):
2051 X_mat = X.natural_representation()
2052 Y_mat = Y.natural_representation()
2053 return (X_mat*Y_mat).trace()
2054
2055
2056 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
2057 """
2058 The rank-n simple EJA consisting of real symmetric n-by-n
2059 matrices, the usual symmetric Jordan product, and the trace inner
2060 product. It has dimension `(n^2 + n)/2` over the reals.
2061
2062 SETUP::
2063
2064 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
2065
2066 EXAMPLES::
2067
2068 sage: J = RealSymmetricEJA(2)
2069 sage: e0, e1, e2 = J.gens()
2070 sage: e0*e0
2071 e0
2072 sage: e1*e1
2073 e0 + e2
2074 sage: e2*e2
2075 e2
2076
2077 TESTS:
2078
2079 The degree of this algebra is `(n^2 + n) / 2`::
2080
2081 sage: set_random_seed()
2082 sage: n = ZZ.random_element(1,5)
2083 sage: J = RealSymmetricEJA(n)
2084 sage: J.degree() == (n^2 + n)/2
2085 True
2086
2087 The Jordan multiplication is what we think it is::
2088
2089 sage: set_random_seed()
2090 sage: n = ZZ.random_element(1,5)
2091 sage: J = RealSymmetricEJA(n)
2092 sage: x = J.random_element()
2093 sage: y = J.random_element()
2094 sage: actual = (x*y).natural_representation()
2095 sage: X = x.natural_representation()
2096 sage: Y = y.natural_representation()
2097 sage: expected = (X*Y + Y*X)/2
2098 sage: actual == expected
2099 True
2100 sage: J(expected) == x*y
2101 True
2102
2103 """
2104 @staticmethod
2105 def __classcall_private__(cls, n, field=QQ):
2106 S = _real_symmetric_basis(n, field=field)
2107 (Qs, T) = _multiplication_table_from_matrix_basis(S)
2108
2109 fdeja = super(RealSymmetricEJA, cls)
2110 return fdeja.__classcall_private__(cls,
2111 field,
2112 Qs,
2113 rank=n,
2114 natural_basis=T)
2115
2116 def inner_product(self, x, y):
2117 return _matrix_ip(x,y)
2118
2119
2120 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
2121 """
2122 The rank-n simple EJA consisting of complex Hermitian n-by-n
2123 matrices over the real numbers, the usual symmetric Jordan product,
2124 and the real-part-of-trace inner product. It has dimension `n^2` over
2125 the reals.
2126
2127 SETUP::
2128
2129 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2130
2131 TESTS:
2132
2133 The degree of this algebra is `n^2`::
2134
2135 sage: set_random_seed()
2136 sage: n = ZZ.random_element(1,5)
2137 sage: J = ComplexHermitianEJA(n)
2138 sage: J.degree() == n^2
2139 True
2140
2141 The Jordan multiplication is what we think it is::
2142
2143 sage: set_random_seed()
2144 sage: n = ZZ.random_element(1,5)
2145 sage: J = ComplexHermitianEJA(n)
2146 sage: x = J.random_element()
2147 sage: y = J.random_element()
2148 sage: actual = (x*y).natural_representation()
2149 sage: X = x.natural_representation()
2150 sage: Y = y.natural_representation()
2151 sage: expected = (X*Y + Y*X)/2
2152 sage: actual == expected
2153 True
2154 sage: J(expected) == x*y
2155 True
2156
2157 """
2158 @staticmethod
2159 def __classcall_private__(cls, n, field=QQ):
2160 S = _complex_hermitian_basis(n)
2161 (Qs, T) = _multiplication_table_from_matrix_basis(S)
2162
2163 fdeja = super(ComplexHermitianEJA, cls)
2164 return fdeja.__classcall_private__(cls,
2165 field,
2166 Qs,
2167 rank=n,
2168 natural_basis=T)
2169
2170 def inner_product(self, x, y):
2171 # Since a+bi on the diagonal is represented as
2172 #
2173 # a + bi = [ a b ]
2174 # [ -b a ],
2175 #
2176 # we'll double-count the "a" entries if we take the trace of
2177 # the embedding.
2178 return _matrix_ip(x,y)/2
2179
2180
2181 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
2182 """
2183 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2184 matrices, the usual symmetric Jordan product, and the
2185 real-part-of-trace inner product. It has dimension `2n^2 - n` over
2186 the reals.
2187
2188 SETUP::
2189
2190 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2191
2192 TESTS:
2193
2194 The degree of this algebra is `n^2`::
2195
2196 sage: set_random_seed()
2197 sage: n = ZZ.random_element(1,5)
2198 sage: J = QuaternionHermitianEJA(n)
2199 sage: J.degree() == 2*(n^2) - n
2200 True
2201
2202 The Jordan multiplication is what we think it is::
2203
2204 sage: set_random_seed()
2205 sage: n = ZZ.random_element(1,5)
2206 sage: J = QuaternionHermitianEJA(n)
2207 sage: x = J.random_element()
2208 sage: y = J.random_element()
2209 sage: actual = (x*y).natural_representation()
2210 sage: X = x.natural_representation()
2211 sage: Y = y.natural_representation()
2212 sage: expected = (X*Y + Y*X)/2
2213 sage: actual == expected
2214 True
2215 sage: J(expected) == x*y
2216 True
2217
2218 """
2219 @staticmethod
2220 def __classcall_private__(cls, n, field=QQ):
2221 S = _quaternion_hermitian_basis(n)
2222 (Qs, T) = _multiplication_table_from_matrix_basis(S)
2223
2224 fdeja = super(QuaternionHermitianEJA, cls)
2225 return fdeja.__classcall_private__(cls,
2226 field,
2227 Qs,
2228 rank=n,
2229 natural_basis=T)
2230
2231 def inner_product(self, x, y):
2232 # Since a+bi+cj+dk on the diagonal is represented as
2233 #
2234 # a + bi +cj + dk = [ a b c d]
2235 # [ -b a -d c]
2236 # [ -c d a -b]
2237 # [ -d -c b a],
2238 #
2239 # we'll quadruple-count the "a" entries if we take the trace of
2240 # the embedding.
2241 return _matrix_ip(x,y)/4
2242
2243
2244 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
2245 """
2246 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2247 with the usual inner product and jordan product ``x*y =
2248 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2249 the reals.
2250
2251 SETUP::
2252
2253 sage: from mjo.eja.eja_algebra import JordanSpinEJA
2254
2255 EXAMPLES:
2256
2257 This multiplication table can be verified by hand::
2258
2259 sage: J = JordanSpinEJA(4)
2260 sage: e0,e1,e2,e3 = J.gens()
2261 sage: e0*e0
2262 e0
2263 sage: e0*e1
2264 e1
2265 sage: e0*e2
2266 e2
2267 sage: e0*e3
2268 e3
2269 sage: e1*e2
2270 0
2271 sage: e1*e3
2272 0
2273 sage: e2*e3
2274 0
2275
2276 """
2277 @staticmethod
2278 def __classcall_private__(cls, n, field=QQ):
2279 Qs = []
2280 id_matrix = matrix.identity(field, n)
2281 for i in xrange(n):
2282 ei = id_matrix.column(i)
2283 Qi = matrix.zero(field, n)
2284 Qi.set_row(0, ei)
2285 Qi.set_column(0, ei)
2286 Qi += matrix.diagonal(n, [ei[0]]*n)
2287 # The addition of the diagonal matrix adds an extra ei[0] in the
2288 # upper-left corner of the matrix.
2289 Qi[0,0] = Qi[0,0] * ~field(2)
2290 Qs.append(Qi)
2291
2292 # The rank of the spin algebra is two, unless we're in a
2293 # one-dimensional ambient space (because the rank is bounded by
2294 # the ambient dimension).
2295 fdeja = super(JordanSpinEJA, cls)
2296 return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2))
2297
2298 def inner_product(self, x, y):
2299 return _usual_ip(x,y)
My personal library of Sage code.