]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/eja/eja_algebra.py
12166267cd34410e0a03616a4b2ad5c7f5b8c105
[sage.d.git] / mjo / eja / eja_algebra.py
1 """
2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
6 """
7
8 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
9 from sage.categories.magmatic_algebras import MagmaticAlgebras
10 from sage.combinat.free_module import CombinatorialFreeModule
11 from sage.matrix.constructor import matrix
12 from sage.matrix.matrix_space import MatrixSpace
13 from sage.misc.cachefunc import cached_method
14 from sage.misc.prandom import choice
15 from sage.misc.table import table
16 from sage.modules.free_module import FreeModule, VectorSpace
17 from sage.rings.integer_ring import ZZ
18 from sage.rings.number_field.number_field import NumberField
19 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
20 from sage.rings.rational_field import QQ
21 from sage.rings.real_lazy import CLF
22 from sage.structure.element import is_Matrix
23
24 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
25 from mjo.eja.eja_utils import _mat2vec
26
27 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
28 # This is an ugly hack needed to prevent the category framework
29 # from implementing a coercion from our base ring (e.g. the
30 # rationals) into the algebra. First of all -- such a coercion is
31 # nonsense to begin with. But more importantly, it tries to do so
32 # in the category of rings, and since our algebras aren't
33 # associative they generally won't be rings.
34 _no_generic_basering_coercion = True
35
36 def __init__(self,
37 field,
38 mult_table,
39 rank,
40 prefix='e',
41 category=None,
42 natural_basis=None):
43 """
44 SETUP::
45
46 sage: from mjo.eja.eja_algebra import random_eja
47
48 EXAMPLES:
49
50 By definition, Jordan multiplication commutes::
51
52 sage: set_random_seed()
53 sage: J = random_eja()
54 sage: x = J.random_element()
55 sage: y = J.random_element()
56 sage: x*y == y*x
57 True
58
59 """
60 self._rank = rank
61 self._natural_basis = natural_basis
62
63 if category is None:
64 category = MagmaticAlgebras(field).FiniteDimensional()
65 category = category.WithBasis().Unital()
66
67 fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
68 fda.__init__(field,
69 range(len(mult_table)),
70 prefix=prefix,
71 category=category)
72 self.print_options(bracket='')
73
74 # The multiplication table we're given is necessarily in terms
75 # of vectors, because we don't have an algebra yet for
76 # anything to be an element of. However, it's faster in the
77 # long run to have the multiplication table be in terms of
78 # algebra elements. We do this after calling the superclass
79 # constructor so that from_vector() knows what to do.
80 self._multiplication_table = [ map(lambda x: self.from_vector(x), ls)
81 for ls in mult_table ]
82
83
84 def _element_constructor_(self, elt):
85 """
86 Construct an element of this algebra from its natural
87 representation.
88
89 This gets called only after the parent element _call_ method
90 fails to find a coercion for the argument.
91
92 SETUP::
93
94 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
95 ....: RealCartesianProductEJA,
96 ....: RealSymmetricEJA)
97
98 EXAMPLES:
99
100 The identity in `S^n` is converted to the identity in the EJA::
101
102 sage: J = RealSymmetricEJA(3)
103 sage: I = matrix.identity(QQ,3)
104 sage: J(I) == J.one()
105 True
106
107 This skew-symmetric matrix can't be represented in the EJA::
108
109 sage: J = RealSymmetricEJA(3)
110 sage: A = matrix(QQ,3, lambda i,j: i-j)
111 sage: J(A)
112 Traceback (most recent call last):
113 ...
114 ArithmeticError: vector is not in free module
115
116 TESTS:
117
118 Ensure that we can convert any element of the two non-matrix
119 simple algebras (whose natural representations are their usual
120 vector representations) back and forth faithfully::
121
122 sage: set_random_seed()
123 sage: J = RealCartesianProductEJA(5)
124 sage: x = J.random_element()
125 sage: J(x.to_vector().column()) == x
126 True
127 sage: J = JordanSpinEJA(5)
128 sage: x = J.random_element()
129 sage: J(x.to_vector().column()) == x
130 True
131
132 """
133 if elt == 0:
134 # The superclass implementation of random_element()
135 # needs to be able to coerce "0" into the algebra.
136 return self.zero()
137
138 natural_basis = self.natural_basis()
139 basis_space = natural_basis[0].matrix_space()
140 if elt not in basis_space:
141 raise ValueError("not a naturally-represented algebra element")
142
143 # Thanks for nothing! Matrix spaces aren't vector spaces in
144 # Sage, so we have to figure out its natural-basis coordinates
145 # ourselves. We use the basis space's ring instead of the
146 # element's ring because the basis space might be an algebraic
147 # closure whereas the base ring of the 3-by-3 identity matrix
148 # could be QQ instead of QQbar.
149 V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols())
150 W = V.span_of_basis( _mat2vec(s) for s in natural_basis )
151 coords = W.coordinate_vector(_mat2vec(elt))
152 return self.from_vector(coords)
153
154
155 def _repr_(self):
156 """
157 Return a string representation of ``self``.
158
159 SETUP::
160
161 sage: from mjo.eja.eja_algebra import JordanSpinEJA
162
163 TESTS:
164
165 Ensure that it says what we think it says::
166
167 sage: JordanSpinEJA(2, field=QQ)
168 Euclidean Jordan algebra of dimension 2 over Rational Field
169 sage: JordanSpinEJA(3, field=RDF)
170 Euclidean Jordan algebra of dimension 3 over Real Double Field
171
172 """
173 fmt = "Euclidean Jordan algebra of dimension {} over {}"
174 return fmt.format(self.dimension(), self.base_ring())
175
176 def product_on_basis(self, i, j):
177 return self._multiplication_table[i][j]
178
179 def _a_regular_element(self):
180 """
181 Guess a regular element. Needed to compute the basis for our
182 characteristic polynomial coefficients.
183
184 SETUP::
185
186 sage: from mjo.eja.eja_algebra import random_eja
187
188 TESTS:
189
190 Ensure that this hacky method succeeds for every algebra that we
191 know how to construct::
192
193 sage: set_random_seed()
194 sage: J = random_eja()
195 sage: J._a_regular_element().is_regular()
196 True
197
198 """
199 gs = self.gens()
200 z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
201 if not z.is_regular():
202 raise ValueError("don't know a regular element")
203 return z
204
205
206 @cached_method
207 def _charpoly_basis_space(self):
208 """
209 Return the vector space spanned by the basis used in our
210 characteristic polynomial coefficients. This is used not only to
211 compute those coefficients, but also any time we need to
212 evaluate the coefficients (like when we compute the trace or
213 determinant).
214 """
215 z = self._a_regular_element()
216 # Don't use the parent vector space directly here in case this
217 # happens to be a subalgebra. In that case, we would be e.g.
218 # two-dimensional but span_of_basis() would expect three
219 # coordinates.
220 V = VectorSpace(self.base_ring(), self.vector_space().dimension())
221 basis = [ (z**k).to_vector() for k in range(self.rank()) ]
222 V1 = V.span_of_basis( basis )
223 b = (V1.basis() + V1.complement().basis())
224 return V.span_of_basis(b)
225
226
227 @cached_method
228 def _charpoly_coeff(self, i):
229 """
230 Return the coefficient polynomial "a_{i}" of this algebra's
231 general characteristic polynomial.
232
233 Having this be a separate cached method lets us compute and
234 store the trace/determinant (a_{r-1} and a_{0} respectively)
235 separate from the entire characteristic polynomial.
236 """
237 (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
238 R = A_of_x.base_ring()
239 if i >= self.rank():
240 # Guaranteed by theory
241 return R.zero()
242
243 # Danger: the in-place modification is done for performance
244 # reasons (reconstructing a matrix with huge polynomial
245 # entries is slow), but I don't know how cached_method works,
246 # so it's highly possible that we're modifying some global
247 # list variable by reference, here. In other words, you
248 # probably shouldn't call this method twice on the same
249 # algebra, at the same time, in two threads
250 Ai_orig = A_of_x.column(i)
251 A_of_x.set_column(i,xr)
252 numerator = A_of_x.det()
253 A_of_x.set_column(i,Ai_orig)
254
255 # We're relying on the theory here to ensure that each a_i is
256 # indeed back in R, and the added negative signs are to make
257 # the whole charpoly expression sum to zero.
258 return R(-numerator/detA)
259
260
261 @cached_method
262 def _charpoly_matrix_system(self):
263 """
264 Compute the matrix whose entries A_ij are polynomials in
265 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
266 corresponding to `x^r` and the determinent of the matrix A =
267 [A_ij]. In other words, all of the fixed (cachable) data needed
268 to compute the coefficients of the characteristic polynomial.
269 """
270 r = self.rank()
271 n = self.dimension()
272
273 # Turn my vector space into a module so that "vectors" can
274 # have multivatiate polynomial entries.
275 names = tuple('X' + str(i) for i in range(1,n+1))
276 R = PolynomialRing(self.base_ring(), names)
277
278 # Using change_ring() on the parent's vector space doesn't work
279 # here because, in a subalgebra, that vector space has a basis
280 # and change_ring() tries to bring the basis along with it. And
281 # that doesn't work unless the new ring is a PID, which it usually
282 # won't be.
283 V = FreeModule(R,n)
284
285 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
286 # indeterminates...
287 x = V(names)
288
289 # And figure out the "left multiplication by x" matrix in
290 # that setting.
291 lmbx_cols = []
292 monomial_matrices = [ self.monomial(i).operator().matrix()
293 for i in range(n) ] # don't recompute these!
294 for k in range(n):
295 ek = self.monomial(k).to_vector()
296 lmbx_cols.append(
297 sum( x[i]*(monomial_matrices[i]*ek)
298 for i in range(n) ) )
299 Lx = matrix.column(R, lmbx_cols)
300
301 # Now we can compute powers of x "symbolically"
302 x_powers = [self.one().to_vector(), x]
303 for d in range(2, r+1):
304 x_powers.append( Lx*(x_powers[-1]) )
305
306 idmat = matrix.identity(R, n)
307
308 W = self._charpoly_basis_space()
309 W = W.change_ring(R.fraction_field())
310
311 # Starting with the standard coordinates x = (X1,X2,...,Xn)
312 # and then converting the entries to W-coordinates allows us
313 # to pass in the standard coordinates to the charpoly and get
314 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
315 # we have
316 #
317 # W.coordinates(x^2) eval'd at (standard z-coords)
318 # =
319 # W-coords of (z^2)
320 # =
321 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
322 #
323 # We want the middle equivalent thing in our matrix, but use
324 # the first equivalent thing instead so that we can pass in
325 # standard coordinates.
326 x_powers = [ W.coordinate_vector(xp) for xp in x_powers ]
327 l2 = [idmat.column(k-1) for k in range(r+1, n+1)]
328 A_of_x = matrix.column(R, n, (x_powers[:r] + l2))
329 return (A_of_x, x, x_powers[r], A_of_x.det())
330
331
332 @cached_method
333 def characteristic_polynomial(self):
334 """
335 Return a characteristic polynomial that works for all elements
336 of this algebra.
337
338 The resulting polynomial has `n+1` variables, where `n` is the
339 dimension of this algebra. The first `n` variables correspond to
340 the coordinates of an algebra element: when evaluated at the
341 coordinates of an algebra element with respect to a certain
342 basis, the result is a univariate polynomial (in the one
343 remaining variable ``t``), namely the characteristic polynomial
344 of that element.
345
346 SETUP::
347
348 sage: from mjo.eja.eja_algebra import JordanSpinEJA
349
350 EXAMPLES:
351
352 The characteristic polynomial in the spin algebra is given in
353 Alizadeh, Example 11.11::
354
355 sage: J = JordanSpinEJA(3)
356 sage: p = J.characteristic_polynomial(); p
357 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
358 sage: xvec = J.one().to_vector()
359 sage: p(*xvec)
360 t^2 - 2*t + 1
361
362 """
363 r = self.rank()
364 n = self.dimension()
365
366 # The list of coefficient polynomials a_1, a_2, ..., a_n.
367 a = [ self._charpoly_coeff(i) for i in range(n) ]
368
369 # We go to a bit of trouble here to reorder the
370 # indeterminates, so that it's easier to evaluate the
371 # characteristic polynomial at x's coordinates and get back
372 # something in terms of t, which is what we want.
373 R = a[0].parent()
374 S = PolynomialRing(self.base_ring(),'t')
375 t = S.gen(0)
376 S = PolynomialRing(S, R.variable_names())
377 t = S(t)
378
379 # Note: all entries past the rth should be zero. The
380 # coefficient of the highest power (x^r) is 1, but it doesn't
381 # appear in the solution vector which contains coefficients
382 # for the other powers (to make them sum to x^r).
383 if (r < n):
384 a[r] = 1 # corresponds to x^r
385 else:
386 # When the rank is equal to the dimension, trying to
387 # assign a[r] goes out-of-bounds.
388 a.append(1) # corresponds to x^r
389
390 return sum( a[k]*(t**k) for k in range(len(a)) )
391
392
393 def inner_product(self, x, y):
394 """
395 The inner product associated with this Euclidean Jordan algebra.
396
397 Defaults to the trace inner product, but can be overridden by
398 subclasses if they are sure that the necessary properties are
399 satisfied.
400
401 SETUP::
402
403 sage: from mjo.eja.eja_algebra import random_eja
404
405 EXAMPLES:
406
407 The inner product must satisfy its axiom for this algebra to truly
408 be a Euclidean Jordan Algebra::
409
410 sage: set_random_seed()
411 sage: J = random_eja()
412 sage: x = J.random_element()
413 sage: y = J.random_element()
414 sage: z = J.random_element()
415 sage: (x*y).inner_product(z) == y.inner_product(x*z)
416 True
417
418 """
419 if (not x in self) or (not y in self):
420 raise TypeError("arguments must live in this algebra")
421 return x.trace_inner_product(y)
422
423
424 def is_trivial(self):
425 """
426 Return whether or not this algebra is trivial.
427
428 A trivial algebra contains only the zero element.
429
430 SETUP::
431
432 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
433
434 EXAMPLES::
435
436 sage: J = ComplexHermitianEJA(3)
437 sage: J.is_trivial()
438 False
439 sage: A = J.zero().subalgebra_generated_by()
440 sage: A.is_trivial()
441 True
442
443 """
444 return self.dimension() == 0
445
446
447 def multiplication_table(self):
448 """
449 Return a visual representation of this algebra's multiplication
450 table (on basis elements).
451
452 SETUP::
453
454 sage: from mjo.eja.eja_algebra import JordanSpinEJA
455
456 EXAMPLES::
457
458 sage: J = JordanSpinEJA(4)
459 sage: J.multiplication_table()
460 +----++----+----+----+----+
461 | * || e0 | e1 | e2 | e3 |
462 +====++====+====+====+====+
463 | e0 || e0 | e1 | e2 | e3 |
464 +----++----+----+----+----+
465 | e1 || e1 | e0 | 0 | 0 |
466 +----++----+----+----+----+
467 | e2 || e2 | 0 | e0 | 0 |
468 +----++----+----+----+----+
469 | e3 || e3 | 0 | 0 | e0 |
470 +----++----+----+----+----+
471
472 """
473 M = list(self._multiplication_table) # copy
474 for i in range(len(M)):
475 # M had better be "square"
476 M[i] = [self.monomial(i)] + M[i]
477 M = [["*"] + list(self.gens())] + M
478 return table(M, header_row=True, header_column=True, frame=True)
479
480
481 def natural_basis(self):
482 """
483 Return a more-natural representation of this algebra's basis.
484
485 Every finite-dimensional Euclidean Jordan Algebra is a direct
486 sum of five simple algebras, four of which comprise Hermitian
487 matrices. This method returns the original "natural" basis
488 for our underlying vector space. (Typically, the natural basis
489 is used to construct the multiplication table in the first place.)
490
491 Note that this will always return a matrix. The standard basis
492 in `R^n` will be returned as `n`-by-`1` column matrices.
493
494 SETUP::
495
496 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
497 ....: RealSymmetricEJA)
498
499 EXAMPLES::
500
501 sage: J = RealSymmetricEJA(2)
502 sage: J.basis()
503 Finite family {0: e0, 1: e1, 2: e2}
504 sage: J.natural_basis()
505 (
506 [1 0] [ 0 1/2*sqrt2] [0 0]
507 [0 0], [1/2*sqrt2 0], [0 1]
508 )
509
510 ::
511
512 sage: J = JordanSpinEJA(2)
513 sage: J.basis()
514 Finite family {0: e0, 1: e1}
515 sage: J.natural_basis()
516 (
517 [1] [0]
518 [0], [1]
519 )
520
521 """
522 if self._natural_basis is None:
523 M = self.natural_basis_space()
524 return tuple( M(b.to_vector()) for b in self.basis() )
525 else:
526 return self._natural_basis
527
528
529 def natural_basis_space(self):
530 """
531 Return the matrix space in which this algebra's natural basis
532 elements live.
533 """
534 if self._natural_basis is None or len(self._natural_basis) == 0:
535 return MatrixSpace(self.base_ring(), self.dimension(), 1)
536 else:
537 return self._natural_basis[0].matrix_space()
538
539
540 @cached_method
541 def one(self):
542 """
543 Return the unit element of this algebra.
544
545 SETUP::
546
547 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
548 ....: random_eja)
549
550 EXAMPLES::
551
552 sage: J = RealCartesianProductEJA(5)
553 sage: J.one()
554 e0 + e1 + e2 + e3 + e4
555
556 TESTS:
557
558 The identity element acts like the identity::
559
560 sage: set_random_seed()
561 sage: J = random_eja()
562 sage: x = J.random_element()
563 sage: J.one()*x == x and x*J.one() == x
564 True
565
566 The matrix of the unit element's operator is the identity::
567
568 sage: set_random_seed()
569 sage: J = random_eja()
570 sage: actual = J.one().operator().matrix()
571 sage: expected = matrix.identity(J.base_ring(), J.dimension())
572 sage: actual == expected
573 True
574
575 """
576 # We can brute-force compute the matrices of the operators
577 # that correspond to the basis elements of this algebra.
578 # If some linear combination of those basis elements is the
579 # algebra identity, then the same linear combination of
580 # their matrices has to be the identity matrix.
581 #
582 # Of course, matrices aren't vectors in sage, so we have to
583 # appeal to the "long vectors" isometry.
584 oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
585
586 # Now we use basis linear algebra to find the coefficients,
587 # of the matrices-as-vectors-linear-combination, which should
588 # work for the original algebra basis too.
589 A = matrix.column(self.base_ring(), oper_vecs)
590
591 # We used the isometry on the left-hand side already, but we
592 # still need to do it for the right-hand side. Recall that we
593 # wanted something that summed to the identity matrix.
594 b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
595
596 # Now if there's an identity element in the algebra, this should work.
597 coeffs = A.solve_right(b)
598 return self.linear_combination(zip(self.gens(), coeffs))
599
600
601 def random_element(self):
602 # Temporary workaround for https://trac.sagemath.org/ticket/28327
603 if self.is_trivial():
604 return self.zero()
605 else:
606 s = super(FiniteDimensionalEuclideanJordanAlgebra, self)
607 return s.random_element()
608
609
610 def rank(self):
611 """
612 Return the rank of this EJA.
613
614 ALGORITHM:
615
616 The author knows of no algorithm to compute the rank of an EJA
617 where only the multiplication table is known. In lieu of one, we
618 require the rank to be specified when the algebra is created,
619 and simply pass along that number here.
620
621 SETUP::
622
623 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
624 ....: RealSymmetricEJA,
625 ....: ComplexHermitianEJA,
626 ....: QuaternionHermitianEJA,
627 ....: random_eja)
628
629 EXAMPLES:
630
631 The rank of the Jordan spin algebra is always two::
632
633 sage: JordanSpinEJA(2).rank()
634 2
635 sage: JordanSpinEJA(3).rank()
636 2
637 sage: JordanSpinEJA(4).rank()
638 2
639
640 The rank of the `n`-by-`n` Hermitian real, complex, or
641 quaternion matrices is `n`::
642
643 sage: RealSymmetricEJA(2).rank()
644 2
645 sage: ComplexHermitianEJA(2).rank()
646 2
647 sage: QuaternionHermitianEJA(2).rank()
648 2
649 sage: RealSymmetricEJA(5).rank()
650 5
651 sage: ComplexHermitianEJA(5).rank()
652 5
653 sage: QuaternionHermitianEJA(5).rank()
654 5
655
656 TESTS:
657
658 Ensure that every EJA that we know how to construct has a
659 positive integer rank::
660
661 sage: set_random_seed()
662 sage: r = random_eja().rank()
663 sage: r in ZZ and r > 0
664 True
665
666 """
667 return self._rank
668
669
670 def vector_space(self):
671 """
672 Return the vector space that underlies this algebra.
673
674 SETUP::
675
676 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
677
678 EXAMPLES::
679
680 sage: J = RealSymmetricEJA(2)
681 sage: J.vector_space()
682 Vector space of dimension 3 over...
683
684 """
685 return self.zero().to_vector().parent().ambient_vector_space()
686
687
688 Element = FiniteDimensionalEuclideanJordanAlgebraElement
689
690
691 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
692 """
693 Return the Euclidean Jordan Algebra corresponding to the set
694 `R^n` under the Hadamard product.
695
696 Note: this is nothing more than the Cartesian product of ``n``
697 copies of the spin algebra. Once Cartesian product algebras
698 are implemented, this can go.
699
700 SETUP::
701
702 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
703
704 EXAMPLES:
705
706 This multiplication table can be verified by hand::
707
708 sage: J = RealCartesianProductEJA(3)
709 sage: e0,e1,e2 = J.gens()
710 sage: e0*e0
711 e0
712 sage: e0*e1
713 0
714 sage: e0*e2
715 0
716 sage: e1*e1
717 e1
718 sage: e1*e2
719 0
720 sage: e2*e2
721 e2
722
723 TESTS:
724
725 We can change the generator prefix::
726
727 sage: RealCartesianProductEJA(3, prefix='r').gens()
728 (r0, r1, r2)
729
730 Our inner product satisfies the Jordan axiom::
731
732 sage: set_random_seed()
733 sage: n = ZZ.random_element(1,5)
734 sage: J = RealCartesianProductEJA(n)
735 sage: x = J.random_element()
736 sage: y = J.random_element()
737 sage: z = J.random_element()
738 sage: (x*y).inner_product(z) == y.inner_product(x*z)
739 True
740
741 """
742 def __init__(self, n, field=QQ, **kwargs):
743 V = VectorSpace(field, n)
744 mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
745 for i in range(n) ]
746
747 fdeja = super(RealCartesianProductEJA, self)
748 return fdeja.__init__(field, mult_table, rank=n, **kwargs)
749
750 def inner_product(self, x, y):
751 return _usual_ip(x,y)
752
753
754 def random_eja():
755 """
756 Return a "random" finite-dimensional Euclidean Jordan Algebra.
757
758 ALGORITHM:
759
760 For now, we choose a random natural number ``n`` (greater than zero)
761 and then give you back one of the following:
762
763 * The cartesian product of the rational numbers ``n`` times; this is
764 ``QQ^n`` with the Hadamard product.
765
766 * The Jordan spin algebra on ``QQ^n``.
767
768 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
769 product.
770
771 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
772 in the space of ``2n``-by-``2n`` real symmetric matrices.
773
774 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
775 in the space of ``4n``-by-``4n`` real symmetric matrices.
776
777 Later this might be extended to return Cartesian products of the
778 EJAs above.
779
780 SETUP::
781
782 sage: from mjo.eja.eja_algebra import random_eja
783
784 TESTS::
785
786 sage: random_eja()
787 Euclidean Jordan algebra of dimension...
788
789 """
790
791 # The max_n component lets us choose different upper bounds on the
792 # value "n" that gets passed to the constructor. This is needed
793 # because e.g. R^{10} is reasonable to test, while the Hermitian
794 # 10-by-10 quaternion matrices are not.
795 (constructor, max_n) = choice([(RealCartesianProductEJA, 6),
796 (JordanSpinEJA, 6),
797 (RealSymmetricEJA, 5),
798 (ComplexHermitianEJA, 4),
799 (QuaternionHermitianEJA, 3)])
800 n = ZZ.random_element(1, max_n)
801 return constructor(n, field=QQ)
802
803
804
805 def _real_symmetric_basis(n, field):
806 """
807 Return a basis for the space of real symmetric n-by-n matrices.
808
809 SETUP::
810
811 sage: from mjo.eja.eja_algebra import _real_symmetric_basis
812
813 TESTS::
814
815 sage: set_random_seed()
816 sage: n = ZZ.random_element(1,5)
817 sage: B = _real_symmetric_basis(n, QQbar)
818 sage: all( M.is_symmetric() for M in B)
819 True
820
821 """
822 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
823 # coordinates.
824 S = []
825 for i in xrange(n):
826 for j in xrange(i+1):
827 Eij = matrix(field, n, lambda k,l: k==i and l==j)
828 if i == j:
829 Sij = Eij
830 else:
831 Sij = Eij + Eij.transpose()
832 # Now normalize it.
833 Sij = Sij / _real_symmetric_matrix_ip(Sij,Sij).sqrt()
834 S.append(Sij)
835 return tuple(S)
836
837
838 def _complex_hermitian_basis(n, field):
839 """
840 Returns a basis for the space of complex Hermitian n-by-n matrices.
841
842 SETUP::
843
844 sage: from mjo.eja.eja_algebra import _complex_hermitian_basis
845
846 TESTS::
847
848 sage: set_random_seed()
849 sage: n = ZZ.random_element(1,5)
850 sage: B = _complex_hermitian_basis(n, QQ)
851 sage: all( M.is_symmetric() for M in B)
852 True
853
854 """
855 R = PolynomialRing(field, 'z')
856 z = R.gen()
857 F = NumberField(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
858 I = F.gen()
859
860 # This is like the symmetric case, but we need to be careful:
861 #
862 # * We want conjugate-symmetry, not just symmetry.
863 # * The diagonal will (as a result) be real.
864 #
865 S = []
866 for i in xrange(n):
867 for j in xrange(i+1):
868 Eij = matrix(field, n, lambda k,l: k==i and l==j)
869 if i == j:
870 Sij = _embed_complex_matrix(Eij)
871 S.append(Sij)
872 else:
873 # Beware, orthogonal but not normalized! The second one
874 # has a minus because it's conjugated.
875 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
876 S.append(Sij_real)
877 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
878 S.append(Sij_imag)
879 return tuple(S)
880
881
882 def _quaternion_hermitian_basis(n, field):
883 """
884 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
885
886 SETUP::
887
888 sage: from mjo.eja.eja_algebra import _quaternion_hermitian_basis
889
890 TESTS::
891
892 sage: set_random_seed()
893 sage: n = ZZ.random_element(1,5)
894 sage: B = _quaternion_hermitian_basis(n, QQbar)
895 sage: all( M.is_symmetric() for M in B )
896 True
897
898 """
899 Q = QuaternionAlgebra(QQ,-1,-1)
900 I,J,K = Q.gens()
901
902 # This is like the symmetric case, but we need to be careful:
903 #
904 # * We want conjugate-symmetry, not just symmetry.
905 # * The diagonal will (as a result) be real.
906 #
907 S = []
908 for i in xrange(n):
909 for j in xrange(i+1):
910 Eij = matrix(Q, n, lambda k,l: k==i and l==j)
911 if i == j:
912 Sij = _embed_quaternion_matrix(Eij)
913 S.append(Sij)
914 else:
915 # Beware, orthogonal but not normalized! The second,
916 # third, and fourth ones have a minus because they're
917 # conjugated.
918 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
919 S.append(Sij_real)
920 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
921 S.append(Sij_I)
922 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
923 S.append(Sij_J)
924 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
925 S.append(Sij_K)
926 return tuple(S)
927
928
929
930 def _multiplication_table_from_matrix_basis(basis):
931 """
932 At least three of the five simple Euclidean Jordan algebras have the
933 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
934 multiplication on the right is matrix multiplication. Given a basis
935 for the underlying matrix space, this function returns a
936 multiplication table (obtained by looping through the basis
937 elements) for an algebra of those matrices.
938 """
939 # In S^2, for example, we nominally have four coordinates even
940 # though the space is of dimension three only. The vector space V
941 # is supposed to hold the entire long vector, and the subspace W
942 # of V will be spanned by the vectors that arise from symmetric
943 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
944 field = basis[0].base_ring()
945 dimension = basis[0].nrows()
946
947 V = VectorSpace(field, dimension**2)
948 W = V.span_of_basis( _mat2vec(s) for s in basis )
949 n = len(basis)
950 mult_table = [[W.zero() for j in range(n)] for i in range(n)]
951 for i in range(n):
952 for j in range(n):
953 mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
954 mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
955
956 return mult_table
957
958
959 def _embed_complex_matrix(M):
960 """
961 Embed the n-by-n complex matrix ``M`` into the space of real
962 matrices of size 2n-by-2n via the map the sends each entry `z = a +
963 bi` to the block matrix ``[[a,b],[-b,a]]``.
964
965 SETUP::
966
967 sage: from mjo.eja.eja_algebra import _embed_complex_matrix
968
969 EXAMPLES::
970
971 sage: R = PolynomialRing(QQ, 'z')
972 sage: z = R.gen()
973 sage: F = NumberField(z**2 + 1, 'i', embedding=CLF(-1).sqrt())
974 sage: x1 = F(4 - 2*i)
975 sage: x2 = F(1 + 2*i)
976 sage: x3 = F(-i)
977 sage: x4 = F(6)
978 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
979 sage: _embed_complex_matrix(M)
980 [ 4 -2| 1 2]
981 [ 2 4|-2 1]
982 [-----+-----]
983 [ 0 -1| 6 0]
984 [ 1 0| 0 6]
985
986 TESTS:
987
988 Embedding is a homomorphism (isomorphism, in fact)::
989
990 sage: set_random_seed()
991 sage: n = ZZ.random_element(5)
992 sage: R = PolynomialRing(QQ, 'z')
993 sage: z = R.gen()
994 sage: F = NumberField(z**2 + 1, 'i', embedding=CLF(-1).sqrt())
995 sage: X = random_matrix(F, n)
996 sage: Y = random_matrix(F, n)
997 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
998 sage: expected = _embed_complex_matrix(X*Y)
999 sage: actual == expected
1000 True
1001
1002 """
1003 n = M.nrows()
1004 if M.ncols() != n:
1005 raise ValueError("the matrix 'M' must be square")
1006 field = M.base_ring()
1007 blocks = []
1008 for z in M.list():
1009 a = z.real()
1010 b = z.imag()
1011 blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1012
1013 # We can drop the imaginaries here.
1014 return matrix.block(field.base_ring(), n, blocks)
1015
1016
1017 def _unembed_complex_matrix(M):
1018 """
1019 The inverse of _embed_complex_matrix().
1020
1021 SETUP::
1022
1023 sage: from mjo.eja.eja_algebra import (_embed_complex_matrix,
1024 ....: _unembed_complex_matrix)
1025
1026 EXAMPLES::
1027
1028 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1029 ....: [-2, 1, -4, 3],
1030 ....: [ 9, 10, 11, 12],
1031 ....: [-10, 9, -12, 11] ])
1032 sage: _unembed_complex_matrix(A)
1033 [ 2*i + 1 4*i + 3]
1034 [ 10*i + 9 12*i + 11]
1035
1036 TESTS:
1037
1038 Unembedding is the inverse of embedding::
1039
1040 sage: set_random_seed()
1041 sage: R = PolynomialRing(QQ, 'z')
1042 sage: z = R.gen()
1043 sage: F = NumberField(z**2 + 1, 'i', embedding=CLF(-1).sqrt())
1044 sage: M = random_matrix(F, 3)
1045 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1046 True
1047
1048 """
1049 n = ZZ(M.nrows())
1050 if M.ncols() != n:
1051 raise ValueError("the matrix 'M' must be square")
1052 if not n.mod(2).is_zero():
1053 raise ValueError("the matrix 'M' must be a complex embedding")
1054
1055 R = PolynomialRing(QQ, 'z')
1056 z = R.gen()
1057 F = NumberField(z**2 + 1, 'i', embedding=CLF(-1).sqrt())
1058 i = F.gen()
1059
1060 # Go top-left to bottom-right (reading order), converting every
1061 # 2-by-2 block we see to a single complex element.
1062 elements = []
1063 for k in xrange(n/2):
1064 for j in xrange(n/2):
1065 submat = M[2*k:2*k+2,2*j:2*j+2]
1066 if submat[0,0] != submat[1,1]:
1067 raise ValueError('bad on-diagonal submatrix')
1068 if submat[0,1] != -submat[1,0]:
1069 raise ValueError('bad off-diagonal submatrix')
1070 z = submat[0,0] + submat[0,1]*i
1071 elements.append(z)
1072
1073 return matrix(F, n/2, elements)
1074
1075
1076 def _embed_quaternion_matrix(M):
1077 """
1078 Embed the n-by-n quaternion matrix ``M`` into the space of real
1079 matrices of size 4n-by-4n by first sending each quaternion entry
1080 `z = a + bi + cj + dk` to the block-complex matrix
1081 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1082 a real matrix.
1083
1084 SETUP::
1085
1086 sage: from mjo.eja.eja_algebra import _embed_quaternion_matrix
1087
1088 EXAMPLES::
1089
1090 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1091 sage: i,j,k = Q.gens()
1092 sage: x = 1 + 2*i + 3*j + 4*k
1093 sage: M = matrix(Q, 1, [[x]])
1094 sage: _embed_quaternion_matrix(M)
1095 [ 1 2 3 4]
1096 [-2 1 -4 3]
1097 [-3 4 1 -2]
1098 [-4 -3 2 1]
1099
1100 Embedding is a homomorphism (isomorphism, in fact)::
1101
1102 sage: set_random_seed()
1103 sage: n = ZZ.random_element(5)
1104 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1105 sage: X = random_matrix(Q, n)
1106 sage: Y = random_matrix(Q, n)
1107 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1108 sage: expected = _embed_quaternion_matrix(X*Y)
1109 sage: actual == expected
1110 True
1111
1112 """
1113 quaternions = M.base_ring()
1114 n = M.nrows()
1115 if M.ncols() != n:
1116 raise ValueError("the matrix 'M' must be square")
1117
1118 R = PolynomialRing(QQ, 'z')
1119 z = R.gen()
1120 F = NumberField(z**2 + 1, 'i', embedding=CLF(-1).sqrt())
1121 i = F.gen()
1122
1123 blocks = []
1124 for z in M.list():
1125 t = z.coefficient_tuple()
1126 a = t[0]
1127 b = t[1]
1128 c = t[2]
1129 d = t[3]
1130 cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1131 [-c + d*i, a - b*i]])
1132 blocks.append(_embed_complex_matrix(cplx_matrix))
1133
1134 # We should have real entries by now, so use the realest field
1135 # we've got for the return value.
1136 return matrix.block(quaternions.base_ring(), n, blocks)
1137
1138
1139 def _unembed_quaternion_matrix(M):
1140 """
1141 The inverse of _embed_quaternion_matrix().
1142
1143 SETUP::
1144
1145 sage: from mjo.eja.eja_algebra import (_embed_quaternion_matrix,
1146 ....: _unembed_quaternion_matrix)
1147
1148 EXAMPLES::
1149
1150 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1151 ....: [-2, 1, -4, 3],
1152 ....: [-3, 4, 1, -2],
1153 ....: [-4, -3, 2, 1]])
1154 sage: _unembed_quaternion_matrix(M)
1155 [1 + 2*i + 3*j + 4*k]
1156
1157 TESTS:
1158
1159 Unembedding is the inverse of embedding::
1160
1161 sage: set_random_seed()
1162 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1163 sage: M = random_matrix(Q, 3)
1164 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1165 True
1166
1167 """
1168 n = ZZ(M.nrows())
1169 if M.ncols() != n:
1170 raise ValueError("the matrix 'M' must be square")
1171 if not n.mod(4).is_zero():
1172 raise ValueError("the matrix 'M' must be a complex embedding")
1173
1174 Q = QuaternionAlgebra(QQ,-1,-1)
1175 i,j,k = Q.gens()
1176
1177 # Go top-left to bottom-right (reading order), converting every
1178 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1179 # quaternion block.
1180 elements = []
1181 for l in xrange(n/4):
1182 for m in xrange(n/4):
1183 submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1184 if submat[0,0] != submat[1,1].conjugate():
1185 raise ValueError('bad on-diagonal submatrix')
1186 if submat[0,1] != -submat[1,0].conjugate():
1187 raise ValueError('bad off-diagonal submatrix')
1188 z = submat[0,0].real() + submat[0,0].imag()*i
1189 z += submat[0,1].real()*j + submat[0,1].imag()*k
1190 elements.append(z)
1191
1192 return matrix(Q, n/4, elements)
1193
1194
1195 # The usual inner product on R^n.
1196 def _usual_ip(x,y):
1197 return x.to_vector().inner_product(y.to_vector())
1198
1199 # The inner product used for the real symmetric simple EJA.
1200 # We keep it as a separate function because e.g. the complex
1201 # algebra uses the same inner product, except divided by 2.
1202 def _matrix_ip(X,Y):
1203 X_mat = X.natural_representation()
1204 Y_mat = Y.natural_representation()
1205 return (X_mat*Y_mat).trace()
1206
1207 def _real_symmetric_matrix_ip(X,Y):
1208 return (X*Y).trace()
1209
1210
1211 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
1212 """
1213 The rank-n simple EJA consisting of real symmetric n-by-n
1214 matrices, the usual symmetric Jordan product, and the trace inner
1215 product. It has dimension `(n^2 + n)/2` over the reals.
1216
1217 SETUP::
1218
1219 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1220
1221 EXAMPLES::
1222
1223 sage: J = RealSymmetricEJA(2)
1224 sage: e0, e1, e2 = J.gens()
1225 sage: e0*e0
1226 e0
1227 sage: e1*e1
1228 1/2*e0 + 1/2*e2
1229 sage: e2*e2
1230 e2
1231
1232 TESTS:
1233
1234 The dimension of this algebra is `(n^2 + n) / 2`::
1235
1236 sage: set_random_seed()
1237 sage: n = ZZ.random_element(1,5)
1238 sage: J = RealSymmetricEJA(n)
1239 sage: J.dimension() == (n^2 + n)/2
1240 True
1241
1242 The Jordan multiplication is what we think it is::
1243
1244 sage: set_random_seed()
1245 sage: n = ZZ.random_element(1,5)
1246 sage: J = RealSymmetricEJA(n)
1247 sage: x = J.random_element()
1248 sage: y = J.random_element()
1249 sage: actual = (x*y).natural_representation()
1250 sage: X = x.natural_representation()
1251 sage: Y = y.natural_representation()
1252 sage: expected = (X*Y + Y*X)/2
1253 sage: actual == expected
1254 True
1255 sage: J(expected) == x*y
1256 True
1257
1258 We can change the generator prefix::
1259
1260 sage: RealSymmetricEJA(3, prefix='q').gens()
1261 (q0, q1, q2, q3, q4, q5)
1262
1263 Our inner product satisfies the Jordan axiom::
1264
1265 sage: set_random_seed()
1266 sage: n = ZZ.random_element(1,5)
1267 sage: J = RealSymmetricEJA(n)
1268 sage: x = J.random_element()
1269 sage: y = J.random_element()
1270 sage: z = J.random_element()
1271 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1272 True
1273
1274 Our basis is normalized with respect to the natural inner product::
1275
1276 sage: set_random_seed()
1277 sage: n = ZZ.random_element(1,5)
1278 sage: J = RealSymmetricEJA(n)
1279 sage: all( b.norm() == 1 for b in J.gens() )
1280 True
1281
1282 Left-multiplication operators are symmetric because they satisfy
1283 the Jordan axiom::
1284
1285 sage: set_random_seed()
1286 sage: n = ZZ.random_element(1,5)
1287 sage: x = RealSymmetricEJA(n).random_element()
1288 sage: x.operator().matrix().is_symmetric()
1289 True
1290
1291 """
1292 def __init__(self, n, field=QQ, **kwargs):
1293 if n > 1:
1294 # We'll need sqrt(2) to normalize the basis, and this
1295 # winds up in the multiplication table, so the whole
1296 # algebra needs to be over the field extension.
1297 R = PolynomialRing(field, 'z')
1298 z = R.gen()
1299 field = NumberField(z**2 - 2, 'sqrt2')
1300
1301 S = _real_symmetric_basis(n, field)
1302 Qs = _multiplication_table_from_matrix_basis(S)
1303
1304 fdeja = super(RealSymmetricEJA, self)
1305 return fdeja.__init__(field,
1306 Qs,
1307 rank=n,
1308 natural_basis=S,
1309 **kwargs)
1310
1311 def inner_product(self, x, y):
1312 X = x.natural_representation()
1313 Y = y.natural_representation()
1314 return _real_symmetric_matrix_ip(X,Y)
1315
1316
1317 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1318 """
1319 The rank-n simple EJA consisting of complex Hermitian n-by-n
1320 matrices over the real numbers, the usual symmetric Jordan product,
1321 and the real-part-of-trace inner product. It has dimension `n^2` over
1322 the reals.
1323
1324 SETUP::
1325
1326 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1327
1328 TESTS:
1329
1330 The dimension of this algebra is `n^2`::
1331
1332 sage: set_random_seed()
1333 sage: n = ZZ.random_element(1,5)
1334 sage: J = ComplexHermitianEJA(n)
1335 sage: J.dimension() == n^2
1336 True
1337
1338 The Jordan multiplication is what we think it is::
1339
1340 sage: set_random_seed()
1341 sage: n = ZZ.random_element(1,5)
1342 sage: J = ComplexHermitianEJA(n)
1343 sage: x = J.random_element()
1344 sage: y = J.random_element()
1345 sage: actual = (x*y).natural_representation()
1346 sage: X = x.natural_representation()
1347 sage: Y = y.natural_representation()
1348 sage: expected = (X*Y + Y*X)/2
1349 sage: actual == expected
1350 True
1351 sage: J(expected) == x*y
1352 True
1353
1354 We can change the generator prefix::
1355
1356 sage: ComplexHermitianEJA(2, prefix='z').gens()
1357 (z0, z1, z2, z3)
1358
1359 Our inner product satisfies the Jordan axiom::
1360
1361 sage: set_random_seed()
1362 sage: n = ZZ.random_element(1,5)
1363 sage: J = ComplexHermitianEJA(n)
1364 sage: x = J.random_element()
1365 sage: y = J.random_element()
1366 sage: z = J.random_element()
1367 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1368 True
1369
1370 """
1371 def __init__(self, n, field=QQ, **kwargs):
1372 S = _complex_hermitian_basis(n, field)
1373 Qs = _multiplication_table_from_matrix_basis(S)
1374
1375 fdeja = super(ComplexHermitianEJA, self)
1376 return fdeja.__init__(field,
1377 Qs,
1378 rank=n,
1379 natural_basis=S,
1380 **kwargs)
1381
1382
1383 def inner_product(self, x, y):
1384 # Since a+bi on the diagonal is represented as
1385 #
1386 # a + bi = [ a b ]
1387 # [ -b a ],
1388 #
1389 # we'll double-count the "a" entries if we take the trace of
1390 # the embedding.
1391 return _matrix_ip(x,y)/2
1392
1393
1394 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1395 """
1396 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1397 matrices, the usual symmetric Jordan product, and the
1398 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1399 the reals.
1400
1401 SETUP::
1402
1403 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1404
1405 TESTS:
1406
1407 The dimension of this algebra is `n^2`::
1408
1409 sage: set_random_seed()
1410 sage: n = ZZ.random_element(1,5)
1411 sage: J = QuaternionHermitianEJA(n)
1412 sage: J.dimension() == 2*(n^2) - n
1413 True
1414
1415 The Jordan multiplication is what we think it is::
1416
1417 sage: set_random_seed()
1418 sage: n = ZZ.random_element(1,5)
1419 sage: J = QuaternionHermitianEJA(n)
1420 sage: x = J.random_element()
1421 sage: y = J.random_element()
1422 sage: actual = (x*y).natural_representation()
1423 sage: X = x.natural_representation()
1424 sage: Y = y.natural_representation()
1425 sage: expected = (X*Y + Y*X)/2
1426 sage: actual == expected
1427 True
1428 sage: J(expected) == x*y
1429 True
1430
1431 We can change the generator prefix::
1432
1433 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1434 (a0, a1, a2, a3, a4, a5)
1435
1436 Our inner product satisfies the Jordan axiom::
1437
1438 sage: set_random_seed()
1439 sage: n = ZZ.random_element(1,5)
1440 sage: J = QuaternionHermitianEJA(n)
1441 sage: x = J.random_element()
1442 sage: y = J.random_element()
1443 sage: z = J.random_element()
1444 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1445 True
1446
1447 """
1448 def __init__(self, n, field=QQ, **kwargs):
1449 S = _quaternion_hermitian_basis(n, field)
1450 Qs = _multiplication_table_from_matrix_basis(S)
1451
1452 fdeja = super(QuaternionHermitianEJA, self)
1453 return fdeja.__init__(field,
1454 Qs,
1455 rank=n,
1456 natural_basis=S,
1457 **kwargs)
1458
1459 def inner_product(self, x, y):
1460 # Since a+bi+cj+dk on the diagonal is represented as
1461 #
1462 # a + bi +cj + dk = [ a b c d]
1463 # [ -b a -d c]
1464 # [ -c d a -b]
1465 # [ -d -c b a],
1466 #
1467 # we'll quadruple-count the "a" entries if we take the trace of
1468 # the embedding.
1469 return _matrix_ip(x,y)/4
1470
1471
1472 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
1473 """
1474 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1475 with the usual inner product and jordan product ``x*y =
1476 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1477 the reals.
1478
1479 SETUP::
1480
1481 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1482
1483 EXAMPLES:
1484
1485 This multiplication table can be verified by hand::
1486
1487 sage: J = JordanSpinEJA(4)
1488 sage: e0,e1,e2,e3 = J.gens()
1489 sage: e0*e0
1490 e0
1491 sage: e0*e1
1492 e1
1493 sage: e0*e2
1494 e2
1495 sage: e0*e3
1496 e3
1497 sage: e1*e2
1498 0
1499 sage: e1*e3
1500 0
1501 sage: e2*e3
1502 0
1503
1504 We can change the generator prefix::
1505
1506 sage: JordanSpinEJA(2, prefix='B').gens()
1507 (B0, B1)
1508
1509 Our inner product satisfies the Jordan axiom::
1510
1511 sage: set_random_seed()
1512 sage: n = ZZ.random_element(1,5)
1513 sage: J = JordanSpinEJA(n)
1514 sage: x = J.random_element()
1515 sage: y = J.random_element()
1516 sage: z = J.random_element()
1517 sage: (x*y).inner_product(z) == y.inner_product(x*z)
1518 True
1519
1520 """
1521 def __init__(self, n, field=QQ, **kwargs):
1522 V = VectorSpace(field, n)
1523 mult_table = [[V.zero() for j in range(n)] for i in range(n)]
1524 for i in range(n):
1525 for j in range(n):
1526 x = V.gen(i)
1527 y = V.gen(j)
1528 x0 = x[0]
1529 xbar = x[1:]
1530 y0 = y[0]
1531 ybar = y[1:]
1532 # z = x*y
1533 z0 = x.inner_product(y)
1534 zbar = y0*xbar + x0*ybar
1535 z = V([z0] + zbar.list())
1536 mult_table[i][j] = z
1537
1538 # The rank of the spin algebra is two, unless we're in a
1539 # one-dimensional ambient space (because the rank is bounded by
1540 # the ambient dimension).
1541 fdeja = super(JordanSpinEJA, self)
1542 return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
1543
1544 def inner_product(self, x, y):
1545 return _usual_ip(x,y)