]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/eja/eja_element.py
projects / sage.d.git / blob
? search:


6a48309499bce637cd039d4c8631ce451f890ff1
[sage.d.git] / mjo / eja / eja_element.py
1 from sage.matrix.constructor import matrix
2 from sage.misc.cachefunc import cached_method
3 from sage.modules.free_module import VectorSpace
4 from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
5
6 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
7 from mjo.eja.eja_utils import _mat2vec, _scale
8
9 class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
10 """
11 An element of a Euclidean Jordan algebra.
12 """
13
14 def __dir__(self):
15 """
16 Oh man, I should not be doing this. This hides the "disabled"
17 methods ``left_matrix`` and ``matrix`` from introspection;
18 in particular it removes them from tab-completion.
19 """
20 return filter(lambda s: s not in ['left_matrix', 'matrix'],
21 dir(self.__class__) )
22
23
24
25
26 def __pow__(self, n):
27 """
28 Return ``self`` raised to the power ``n``.
29
30 Jordan algebras are always power-associative; see for
31 example Faraut and Korányi, Proposition II.1.2 (ii).
32
33 We have to override this because our superclass uses row
34 vectors instead of column vectors! We, on the other hand,
35 assume column vectors everywhere.
36
37 SETUP::
38
39 sage: from mjo.eja.eja_algebra import random_eja
40
41 TESTS:
42
43 The definition of `x^2` is the unambiguous `x*x`::
44
45 sage: set_random_seed()
46 sage: x = random_eja().random_element()
47 sage: x*x == (x^2)
48 True
49
50 A few examples of power-associativity::
51
52 sage: set_random_seed()
53 sage: x = random_eja().random_element()
54 sage: x*(x*x)*(x*x) == x^5
55 True
56 sage: (x*x)*(x*x*x) == x^5
57 True
58
59 We also know that powers operator-commute (Koecher, Chapter
60 III, Corollary 1)::
61
62 sage: set_random_seed()
63 sage: x = random_eja().random_element()
64 sage: m = ZZ.random_element(0,10)
65 sage: n = ZZ.random_element(0,10)
66 sage: Lxm = (x^m).operator()
67 sage: Lxn = (x^n).operator()
68 sage: Lxm*Lxn == Lxn*Lxm
69 True
70
71 """
72 if n == 0:
73 return self.parent().one()
74 elif n == 1:
75 return self
76 else:
77 return (self**(n-1))*self
78
79
80 def apply_univariate_polynomial(self, p):
81 """
82 Apply the univariate polynomial ``p`` to this element.
83
84 A priori, SageMath won't allow us to apply a univariate
85 polynomial to an element of an EJA, because we don't know
86 that EJAs are rings (they are usually not associative). Of
87 course, we know that EJAs are power-associative, so the
88 operation is ultimately kosher. This function sidesteps
89 the CAS to get the answer we want and expect.
90
91 SETUP::
92
93 sage: from mjo.eja.eja_algebra import (HadamardEJA,
94 ....: random_eja)
95
96 EXAMPLES::
97
98 sage: R = PolynomialRing(QQ, 't')
99 sage: t = R.gen(0)
100 sage: p = t^4 - t^3 + 5*t - 2
101 sage: J = HadamardEJA(5)
102 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
103 True
104
105 TESTS:
106
107 We should always get back an element of the algebra::
108
109 sage: set_random_seed()
110 sage: p = PolynomialRing(AA, 't').random_element()
111 sage: J = random_eja()
112 sage: x = J.random_element()
113 sage: x.apply_univariate_polynomial(p) in J
114 True
115
116 """
117 if len(p.variables()) > 1:
118 raise ValueError("not a univariate polynomial")
119 P = self.parent()
120 R = P.base_ring()
121 # Convert the coeficcients to the parent's base ring,
122 # because a priori they might live in an (unnecessarily)
123 # larger ring for which P.sum() would fail below.
124 cs = [ R(c) for c in p.coefficients(sparse=False) ]
125 return P.sum( cs[k]*(self**k) for k in range(len(cs)) )
126
127
128 def characteristic_polynomial(self):
129 """
130 Return the characteristic polynomial of this element.
131
132 SETUP::
133
134 sage: from mjo.eja.eja_algebra import HadamardEJA
135
136 EXAMPLES:
137
138 The rank of `R^3` is three, and the minimal polynomial of
139 the identity element is `(t-1)` from which it follows that
140 the characteristic polynomial should be `(t-1)^3`::
141
142 sage: J = HadamardEJA(3)
143 sage: J.one().characteristic_polynomial()
144 t^3 - 3*t^2 + 3*t - 1
145
146 Likewise, the characteristic of the zero element in the
147 rank-three algebra `R^{n}` should be `t^{3}`::
148
149 sage: J = HadamardEJA(3)
150 sage: J.zero().characteristic_polynomial()
151 t^3
152
153 TESTS:
154
155 The characteristic polynomial of an element should evaluate
156 to zero on that element::
157
158 sage: set_random_seed()
159 sage: x = HadamardEJA(3).random_element()
160 sage: p = x.characteristic_polynomial()
161 sage: x.apply_univariate_polynomial(p)
162 0
163
164 The characteristic polynomials of the zero and unit elements
165 should be what we think they are in a subalgebra, too::
166
167 sage: J = HadamardEJA(3)
168 sage: p1 = J.one().characteristic_polynomial()
169 sage: q1 = J.zero().characteristic_polynomial()
170 sage: b0,b1,b2 = J.gens()
171 sage: A = (b0 + 2*b1 + 3*b2).subalgebra_generated_by() # dim 3
172 sage: p2 = A.one().characteristic_polynomial()
173 sage: q2 = A.zero().characteristic_polynomial()
174 sage: p1 == p2
175 True
176 sage: q1 == q2
177 True
178
179 """
180 p = self.parent().characteristic_polynomial_of()
181 return p(*self.to_vector())
182
183
184 def inner_product(self, other):
185 """
186 Return the parent algebra's inner product of myself and ``other``.
187
188 SETUP::
189
190 sage: from mjo.eja.eja_algebra import (
191 ....: ComplexHermitianEJA,
192 ....: JordanSpinEJA,
193 ....: QuaternionHermitianEJA,
194 ....: RealSymmetricEJA,
195 ....: random_eja)
196
197 EXAMPLES:
198
199 The inner product in the Jordan spin algebra is the usual
200 inner product on `R^n` (this example only works because the
201 basis for the Jordan algebra is the standard basis in `R^n`)::
202
203 sage: J = JordanSpinEJA(3)
204 sage: x = vector(QQ,[1,2,3])
205 sage: y = vector(QQ,[4,5,6])
206 sage: x.inner_product(y)
207 32
208 sage: J.from_vector(x).inner_product(J.from_vector(y))
209 32
210
211 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
212 multiplication is the usual matrix multiplication in `S^n`,
213 so the inner product of the identity matrix with itself
214 should be the `n`::
215
216 sage: J = RealSymmetricEJA(3)
217 sage: J.one().inner_product(J.one())
218 3
219
220 Likewise, the inner product on `C^n` is `<X,Y> =
221 Re(trace(X*Y))`, where we must necessarily take the real
222 part because the product of Hermitian matrices may not be
223 Hermitian::
224
225 sage: J = ComplexHermitianEJA(3)
226 sage: J.one().inner_product(J.one())
227 3
228
229 Ditto for the quaternions::
230
231 sage: J = QuaternionHermitianEJA(2)
232 sage: J.one().inner_product(J.one())
233 2
234
235 TESTS:
236
237 Ensure that we can always compute an inner product, and that
238 it gives us back a real number::
239
240 sage: set_random_seed()
241 sage: J = random_eja()
242 sage: x,y = J.random_elements(2)
243 sage: x.inner_product(y) in RLF
244 True
245
246 """
247 P = self.parent()
248 if not other in P:
249 raise TypeError("'other' must live in the same algebra")
250
251 return P.inner_product(self, other)
252
253
254 def operator_commutes_with(self, other):
255 """
256 Return whether or not this element operator-commutes
257 with ``other``.
258
259 SETUP::
260
261 sage: from mjo.eja.eja_algebra import random_eja
262
263 EXAMPLES:
264
265 The definition of a Jordan algebra says that any element
266 operator-commutes with its square::
267
268 sage: set_random_seed()
269 sage: x = random_eja().random_element()
270 sage: x.operator_commutes_with(x^2)
271 True
272
273 TESTS:
274
275 Test Lemma 1 from Chapter III of Koecher::
276
277 sage: set_random_seed()
278 sage: u,v = random_eja().random_elements(2)
279 sage: lhs = u.operator_commutes_with(u*v)
280 sage: rhs = v.operator_commutes_with(u^2)
281 sage: lhs == rhs
282 True
283
284 Test the first polarization identity from my notes, Koecher
285 Chapter III, or from Baes (2.3)::
286
287 sage: set_random_seed()
288 sage: x,y = random_eja().random_elements(2)
289 sage: Lx = x.operator()
290 sage: Ly = y.operator()
291 sage: Lxx = (x*x).operator()
292 sage: Lxy = (x*y).operator()
293 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
294 True
295
296 Test the second polarization identity from my notes or from
297 Baes (2.4)::
298
299 sage: set_random_seed()
300 sage: x,y,z = random_eja().random_elements(3)
301 sage: Lx = x.operator()
302 sage: Ly = y.operator()
303 sage: Lz = z.operator()
304 sage: Lzy = (z*y).operator()
305 sage: Lxy = (x*y).operator()
306 sage: Lxz = (x*z).operator()
307 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
308 True
309
310 Test the third polarization identity from my notes or from
311 Baes (2.5)::
312
313 sage: set_random_seed()
314 sage: u,y,z = random_eja().random_elements(3)
315 sage: Lu = u.operator()
316 sage: Ly = y.operator()
317 sage: Lz = z.operator()
318 sage: Lzy = (z*y).operator()
319 sage: Luy = (u*y).operator()
320 sage: Luz = (u*z).operator()
321 sage: Luyz = (u*(y*z)).operator()
322 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
323 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
324 sage: bool(lhs == rhs)
325 True
326
327 """
328 if not other in self.parent():
329 raise TypeError("'other' must live in the same algebra")
330
331 A = self.operator()
332 B = other.operator()
333 return (A*B == B*A)
334
335
336 def det(self):
337 """
338 Return my determinant, the product of my eigenvalues.
339
340 SETUP::
341
342 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
343 ....: TrivialEJA,
344 ....: RealSymmetricEJA,
345 ....: ComplexHermitianEJA,
346 ....: random_eja)
347
348 EXAMPLES::
349
350 sage: J = JordanSpinEJA(2)
351 sage: x = sum( J.gens() )
352 sage: x.det()
353 0
354
355 ::
356
357 sage: J = JordanSpinEJA(3)
358 sage: x = sum( J.gens() )
359 sage: x.det()
360 -1
361
362 The determinant of the sole element in the rank-zero trivial
363 algebra is ``1``, by three paths of reasoning. First, its
364 characteristic polynomial is a constant ``1``, so the constant
365 term in that polynomial is ``1``. Second, the characteristic
366 polynomial evaluated at zero is again ``1``. And finally, the
367 (empty) product of its eigenvalues is likewise just unity::
368
369 sage: J = TrivialEJA()
370 sage: J.zero().det()
371 1
372
373 TESTS:
374
375 An element is invertible if and only if its determinant is
376 non-zero::
377
378 sage: set_random_seed()
379 sage: x = random_eja().random_element()
380 sage: x.is_invertible() == (x.det() != 0)
381 True
382
383 Ensure that the determinant is multiplicative on an associative
384 subalgebra as in Faraut and Korányi's Proposition II.2.2::
385
386 sage: set_random_seed()
387 sage: J = random_eja().random_element().subalgebra_generated_by()
388 sage: x,y = J.random_elements(2)
389 sage: (x*y).det() == x.det()*y.det()
390 True
391
392 The determinant in matrix algebras is just the usual determinant::
393
394 sage: set_random_seed()
395 sage: X = matrix.random(QQ,3)
396 sage: X = X + X.T
397 sage: J1 = RealSymmetricEJA(3)
398 sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
399 sage: expected = X.det()
400 sage: actual1 = J1(X).det()
401 sage: actual2 = J2(X).det()
402 sage: actual1 == expected
403 True
404 sage: actual2 == expected
405 True
406
407 ::
408
409 sage: set_random_seed()
410 sage: J1 = ComplexHermitianEJA(2)
411 sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
412 sage: X = matrix.random(GaussianIntegers(), 2)
413 sage: X = X + X.H
414 sage: expected = AA(X.det())
415 sage: actual1 = J1(J1.real_embed(X)).det()
416 sage: actual2 = J2(J2.real_embed(X)).det()
417 sage: expected == actual1
418 True
419 sage: expected == actual2
420 True
421
422 """
423 P = self.parent()
424 r = P.rank()
425
426 if r == 0:
427 # Special case, since we don't get the a0=1
428 # coefficient when the rank of the algebra
429 # is zero.
430 return P.base_ring().one()
431
432 p = P._charpoly_coefficients()[0]
433 # The _charpoly_coeff function already adds the factor of -1
434 # to ensure that _charpoly_coefficients()[0] is really what
435 # appears in front of t^{0} in the charpoly. However, we want
436 # (-1)^r times THAT for the determinant.
437 return ((-1)**r)*p(*self.to_vector())
438
439
440 @cached_method
441 def inverse(self):
442 """
443 Return the Jordan-multiplicative inverse of this element.
444
445 ALGORITHM:
446
447 In general we appeal to the quadratic representation as in
448 Koecher's Theorem 12 in Chapter III, Section 5. But if the
449 parent algebra's "characteristic polynomial of" coefficients
450 happen to be cached, then we use Proposition II.2.4 in Faraut
451 and Korányi which gives a formula for the inverse based on the
452 characteristic polynomial and the Cayley-Hamilton theorem for
453 Euclidean Jordan algebras::
454
455 SETUP::
456
457 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
458 ....: JordanSpinEJA,
459 ....: random_eja)
460
461 EXAMPLES:
462
463 The inverse in the spin factor algebra is given in Alizadeh's
464 Example 11.11::
465
466 sage: set_random_seed()
467 sage: J = JordanSpinEJA.random_instance()
468 sage: x = J.random_element()
469 sage: while not x.is_invertible():
470 ....: x = J.random_element()
471 sage: x_vec = x.to_vector()
472 sage: x0 = x_vec[:1]
473 sage: x_bar = x_vec[1:]
474 sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
475 sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
476 sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
477 sage: x.inverse() == J.from_vector(x_inverse)
478 True
479
480 Trying to invert a non-invertible element throws an error:
481
482 sage: JordanSpinEJA(3).zero().inverse()
483 Traceback (most recent call last):
484 ...
485 ZeroDivisionError: element is not invertible
486
487 TESTS:
488
489 The identity element is its own inverse::
490
491 sage: set_random_seed()
492 sage: J = random_eja()
493 sage: J.one().inverse() == J.one()
494 True
495
496 If an element has an inverse, it acts like one::
497
498 sage: set_random_seed()
499 sage: J = random_eja()
500 sage: x = J.random_element()
501 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
502 True
503
504 The inverse of the inverse is what we started with::
505
506 sage: set_random_seed()
507 sage: J = random_eja()
508 sage: x = J.random_element()
509 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
510 True
511
512 Proposition II.2.3 in Faraut and Korányi says that the inverse
513 of an element is the inverse of its left-multiplication operator
514 applied to the algebra's identity, when that inverse exists::
515
516 sage: set_random_seed()
517 sage: J = random_eja()
518 sage: x = J.random_element()
519 sage: (not x.operator().is_invertible()) or (
520 ....: x.operator().inverse()(J.one()) == x.inverse() )
521 True
522
523 Check that the fast (cached) and slow algorithms give the same
524 answer::
525
526 sage: set_random_seed() # long time
527 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
528 sage: x = J.random_element() # long time
529 sage: while not x.is_invertible(): # long time
530 ....: x = J.random_element() # long time
531 sage: slow = x.inverse() # long time
532 sage: _ = J._charpoly_coefficients() # long time
533 sage: fast = x.inverse() # long time
534 sage: slow == fast # long time
535 True
536 """
537 not_invertible_msg = "element is not invertible"
538 if self.parent()._charpoly_coefficients.is_in_cache():
539 # We can invert using our charpoly if it will be fast to
540 # compute. If the coefficients are cached, our rank had
541 # better be too!
542 if self.det().is_zero():
543 raise ZeroDivisionError(not_invertible_msg)
544 r = self.parent().rank()
545 a = self.characteristic_polynomial().coefficients(sparse=False)
546 return (-1)**(r+1)*sum(a[i+1]*self**i for i in range(r))/self.det()
547
548 try:
549 inv = (~self.quadratic_representation())(self)
550 self.is_invertible.set_cache(True)
551 return inv
552 except ZeroDivisionError:
553 self.is_invertible.set_cache(False)
554 raise ZeroDivisionError(not_invertible_msg)
555
556
557 @cached_method
558 def is_invertible(self):
559 """
560 Return whether or not this element is invertible.
561
562 ALGORITHM:
563
564 If computing my determinant will be fast, we do so and compare
565 with zero (Proposition II.2.4 in Faraut and
566 Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
567 reduces the problem to the invertibility of my quadratic
568 representation.
569
570 SETUP::
571
572 sage: from mjo.eja.eja_algebra import random_eja
573
574 TESTS:
575
576 The identity element is always invertible::
577
578 sage: set_random_seed()
579 sage: J = random_eja()
580 sage: J.one().is_invertible()
581 True
582
583 The zero element is never invertible in a non-trivial algebra::
584
585 sage: set_random_seed()
586 sage: J = random_eja()
587 sage: (not J.is_trivial()) and J.zero().is_invertible()
588 False
589
590 Test that the fast (cached) and slow algorithms give the same
591 answer::
592
593 sage: set_random_seed() # long time
594 sage: J = random_eja(field=QQ, orthonormalize=False) # long time
595 sage: x = J.random_element() # long time
596 sage: slow = x.is_invertible() # long time
597 sage: _ = J._charpoly_coefficients() # long time
598 sage: fast = x.is_invertible() # long time
599 sage: slow == fast # long time
600 True
601 """
602 if self.is_zero():
603 if self.parent().is_trivial():
604 return True
605 else:
606 return False
607
608 if self.parent()._charpoly_coefficients.is_in_cache():
609 # The determinant will be quicker than inverting the
610 # quadratic representation, most likely.
611 return (not self.det().is_zero())
612
613 # The easiest way to determine if I'm invertible is to try.
614 try:
615 inv = (~self.quadratic_representation())(self)
616 self.inverse.set_cache(inv)
617 return True
618 except ZeroDivisionError:
619 return False
620
621
622 def is_primitive_idempotent(self):
623 """
624 Return whether or not this element is a primitive (or minimal)
625 idempotent.
626
627 A primitive idempotent is a non-zero idempotent that is not
628 the sum of two other non-zero idempotents. Remark 2.7.15 in
629 Baes shows that this is what he refers to as a "minimal
630 idempotent."
631
632 An element of a Euclidean Jordan algebra is a minimal idempotent
633 if it :meth:`is_idempotent` and if its Peirce subalgebra
634 corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
635 Proposition 2.7.17).
636
637 SETUP::
638
639 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
640 ....: RealSymmetricEJA,
641 ....: TrivialEJA,
642 ....: random_eja)
643
644 WARNING::
645
646 This method is sloooooow.
647
648 EXAMPLES:
649
650 The spectral decomposition of a non-regular element should always
651 contain at least one non-minimal idempotent::
652
653 sage: J = RealSymmetricEJA(3)
654 sage: x = sum(J.gens())
655 sage: x.is_regular()
656 False
657 sage: [ c.is_primitive_idempotent()
658 ....: for (l,c) in x.spectral_decomposition() ]
659 [False, True]
660
661 On the other hand, the spectral decomposition of a regular
662 element should always be in terms of minimal idempotents::
663
664 sage: J = JordanSpinEJA(4)
665 sage: x = sum( i*J.monomial(i) for i in range(len(J.gens())) )
666 sage: x.is_regular()
667 True
668 sage: [ c.is_primitive_idempotent()
669 ....: for (l,c) in x.spectral_decomposition() ]
670 [True, True]
671
672 TESTS:
673
674 The identity element is minimal only in an EJA of rank one::
675
676 sage: set_random_seed()
677 sage: J = random_eja()
678 sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
679 True
680
681 A non-idempotent cannot be a minimal idempotent::
682
683 sage: set_random_seed()
684 sage: J = JordanSpinEJA(4)
685 sage: x = J.random_element()
686 sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
687 False
688
689 Proposition 2.7.19 in Baes says that an element is a minimal
690 idempotent if and only if it's idempotent with trace equal to
691 unity::
692
693 sage: set_random_seed()
694 sage: J = JordanSpinEJA(4)
695 sage: x = J.random_element()
696 sage: expected = (x.is_idempotent() and x.trace() == 1)
697 sage: actual = x.is_primitive_idempotent()
698 sage: actual == expected
699 True
700
701 Primitive idempotents must be non-zero::
702
703 sage: set_random_seed()
704 sage: J = random_eja()
705 sage: J.zero().is_idempotent()
706 True
707 sage: J.zero().is_primitive_idempotent()
708 False
709
710 As a consequence of the fact that primitive idempotents must
711 be non-zero, there are no primitive idempotents in a trivial
712 Euclidean Jordan algebra::
713
714 sage: J = TrivialEJA()
715 sage: J.one().is_idempotent()
716 True
717 sage: J.one().is_primitive_idempotent()
718 False
719
720 """
721 if not self.is_idempotent():
722 return False
723
724 if self.is_zero():
725 return False
726
727 (_,_,J1) = self.parent().peirce_decomposition(self)
728 return (J1.dimension() == 1)
729
730
731 def is_nilpotent(self):
732 """
733 Return whether or not some power of this element is zero.
734
735 ALGORITHM:
736
737 We use Theorem 5 in Chapter III of Koecher, which says that
738 an element ``x`` is nilpotent if and only if ``x.operator()``
739 is nilpotent. And it is a basic fact of linear algebra that
740 an operator on an `n`-dimensional space is nilpotent if and
741 only if, when raised to the `n`th power, it equals the zero
742 operator (for example, see Axler Corollary 8.8).
743
744 SETUP::
745
746 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
747 ....: random_eja)
748
749 EXAMPLES::
750
751 sage: J = JordanSpinEJA(3)
752 sage: x = sum(J.gens())
753 sage: x.is_nilpotent()
754 False
755
756 TESTS:
757
758 The identity element is never nilpotent, except in a trivial EJA::
759
760 sage: set_random_seed()
761 sage: J = random_eja()
762 sage: J.one().is_nilpotent() and not J.is_trivial()
763 False
764
765 The additive identity is always nilpotent::
766
767 sage: set_random_seed()
768 sage: random_eja().zero().is_nilpotent()
769 True
770
771 """
772 P = self.parent()
773 zero_operator = P.zero().operator()
774 return self.operator()**P.dimension() == zero_operator
775
776
777 def is_regular(self):
778 """
779 Return whether or not this is a regular element.
780
781 SETUP::
782
783 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
784 ....: random_eja)
785
786 EXAMPLES:
787
788 The identity element always has degree one, but any element
789 linearly-independent from it is regular::
790
791 sage: J = JordanSpinEJA(5)
792 sage: J.one().is_regular()
793 False
794 sage: b0, b1, b2, b3, b4 = J.gens()
795 sage: b0 == J.one()
796 True
797 sage: for x in J.gens():
798 ....: (J.one() + x).is_regular()
799 False
800 True
801 True
802 True
803 True
804
805 TESTS:
806
807 The zero element should never be regular, unless the parent
808 algebra has dimension less than or equal to one::
809
810 sage: set_random_seed()
811 sage: J = random_eja()
812 sage: J.dimension() <= 1 or not J.zero().is_regular()
813 True
814
815 The unit element isn't regular unless the algebra happens to
816 consist of only its scalar multiples::
817
818 sage: set_random_seed()
819 sage: J = random_eja()
820 sage: J.dimension() <= 1 or not J.one().is_regular()
821 True
822
823 """
824 return self.degree() == self.parent().rank()
825
826
827 def degree(self):
828 """
829 Return the degree of this element, which is defined to be
830 the degree of its minimal polynomial.
831
832 ALGORITHM:
833
834 .........
835
836 SETUP::
837
838 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
839 ....: random_eja)
840
841 EXAMPLES::
842
843 sage: J = JordanSpinEJA(4)
844 sage: J.one().degree()
845 1
846 sage: b0,b1,b2,b3 = J.gens()
847 sage: (b0 - b1).degree()
848 2
849
850 In the spin factor algebra (of rank two), all elements that
851 aren't multiples of the identity are regular::
852
853 sage: set_random_seed()
854 sage: J = JordanSpinEJA.random_instance()
855 sage: n = J.dimension()
856 sage: x = J.random_element()
857 sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
858 True
859
860 TESTS:
861
862 The zero and unit elements are both of degree one in nontrivial
863 algebras::
864
865 sage: set_random_seed()
866 sage: J = random_eja()
867 sage: d = J.zero().degree()
868 sage: (J.is_trivial() and d == 0) or d == 1
869 True
870 sage: d = J.one().degree()
871 sage: (J.is_trivial() and d == 0) or d == 1
872 True
873
874 Our implementation agrees with the definition::
875
876 sage: set_random_seed()
877 sage: x = random_eja().random_element()
878 sage: x.degree() == x.minimal_polynomial().degree()
879 True
880
881 """
882 n = self.parent().dimension()
883
884 if n == 0:
885 # The minimal polynomial is an empty product, i.e. the
886 # constant polynomial "1" having degree zero.
887 return 0
888 elif self.is_zero():
889 # The minimal polynomial of zero in a nontrivial algebra
890 # is "t", and is of degree one.
891 return 1
892 elif n == 1:
893 # If this is a nonzero element of a nontrivial algebra, it
894 # has degree at least one. It follows that, in an algebra
895 # of dimension one, the degree must be actually one.
896 return 1
897
898 # BEWARE: The subalgebra_generated_by() method uses the result
899 # of this method to construct a basis for the subalgebra. That
900 # means, in particular, that we cannot implement this method
901 # as ``self.subalgebra_generated_by().dimension()``.
902
903 # Algorithm: keep appending (vector representations of) powers
904 # self as rows to a matrix and echelonizing it. When its rank
905 # stops increasing, we've reached a redundancy.
906
907 # Given the special cases above, we can assume that "self" is
908 # nonzero, the algebra is nontrivial, and that its dimension
909 # is at least two.
910 M = matrix([(self.parent().one()).to_vector()])
911 old_rank = 1
912
913 # Specifying the row-reduction algorithm can e.g. help over
914 # AA because it avoids the RecursionError that gets thrown
915 # when we have to look too hard for a root.
916 #
917 # Beware: QQ supports an entirely different set of "algorithm"
918 # keywords than do AA and RR.
919 algo = None
920 from sage.rings.all import QQ
921 if self.parent().base_ring() is not QQ:
922 algo = "scaled_partial_pivoting"
923
924 for d in range(1,n):
925 M = matrix(M.rows() + [(self**d).to_vector()])
926 M.echelonize(algo)
927 new_rank = M.rank()
928 if new_rank == old_rank:
929 return new_rank
930 else:
931 old_rank = new_rank
932
933 return n
934
935
936
937 def left_matrix(self):
938 """
939 Our parent class defines ``left_matrix`` and ``matrix``
940 methods whose names are misleading. We don't want them.
941 """
942 raise NotImplementedError("use operator().matrix() instead")
943
944 matrix = left_matrix
945
946
947 def minimal_polynomial(self):
948 """
949 Return the minimal polynomial of this element,
950 as a function of the variable `t`.
951
952 ALGORITHM:
953
954 We restrict ourselves to the associative subalgebra
955 generated by this element, and then return the minimal
956 polynomial of this element's operator matrix (in that
957 subalgebra). This works by Baes Proposition 2.3.16.
958
959 SETUP::
960
961 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
962 ....: RealSymmetricEJA,
963 ....: TrivialEJA,
964 ....: random_eja)
965
966 EXAMPLES:
967
968 Keeping in mind that the polynomial ``1`` evaluates the identity
969 element (also the zero element) of the trivial algebra, it is clear
970 that the polynomial ``1`` is the minimal polynomial of the only
971 element in a trivial algebra::
972
973 sage: J = TrivialEJA()
974 sage: J.one().minimal_polynomial()
975 1
976 sage: J.zero().minimal_polynomial()
977 1
978
979 TESTS:
980
981 The minimal polynomial of the identity and zero elements are
982 always the same, except in trivial algebras where the minimal
983 polynomial of the unit/zero element is ``1``::
984
985 sage: set_random_seed()
986 sage: J = random_eja()
987 sage: mu = J.one().minimal_polynomial()
988 sage: t = mu.parent().gen()
989 sage: mu + int(J.is_trivial())*(t-2)
990 t - 1
991 sage: mu = J.zero().minimal_polynomial()
992 sage: t = mu.parent().gen()
993 sage: mu + int(J.is_trivial())*(t-1)
994 t
995
996 The degree of an element is (by one definition) the degree
997 of its minimal polynomial::
998
999 sage: set_random_seed()
1000 sage: x = random_eja().random_element()
1001 sage: x.degree() == x.minimal_polynomial().degree()
1002 True
1003
1004 The minimal polynomial and the characteristic polynomial coincide
1005 and are known (see Alizadeh, Example 11.11) for all elements of
1006 the spin factor algebra that aren't scalar multiples of the
1007 identity. We require the dimension of the algebra to be at least
1008 two here so that said elements actually exist::
1009
1010 sage: set_random_seed()
1011 sage: n_max = max(2, JordanSpinEJA._max_random_instance_size())
1012 sage: n = ZZ.random_element(2, n_max)
1013 sage: J = JordanSpinEJA(n)
1014 sage: y = J.random_element()
1015 sage: while y == y.coefficient(0)*J.one():
1016 ....: y = J.random_element()
1017 sage: y0 = y.to_vector()[0]
1018 sage: y_bar = y.to_vector()[1:]
1019 sage: actual = y.minimal_polynomial()
1020 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1021 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1022 sage: bool(actual == expected)
1023 True
1024
1025 The minimal polynomial should always kill its element::
1026
1027 sage: set_random_seed()
1028 sage: x = random_eja().random_element()
1029 sage: p = x.minimal_polynomial()
1030 sage: x.apply_univariate_polynomial(p)
1031 0
1032
1033 The minimal polynomial is invariant under a change of basis,
1034 and in particular, a re-scaling of the basis::
1035
1036 sage: set_random_seed()
1037 sage: n_max = RealSymmetricEJA._max_random_instance_size()
1038 sage: n = ZZ.random_element(1, n_max)
1039 sage: J1 = RealSymmetricEJA(n)
1040 sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
1041 sage: X = random_matrix(AA,n)
1042 sage: X = X*X.transpose()
1043 sage: x1 = J1(X)
1044 sage: x2 = J2(X)
1045 sage: x1.minimal_polynomial() == x2.minimal_polynomial()
1046 True
1047
1048 """
1049 if self.is_zero():
1050 # Pretty sure we know what the minimal polynomial of
1051 # the zero operator is going to be. This ensures
1052 # consistency of e.g. the polynomial variable returned
1053 # in the "normal" case without us having to think about it.
1054 return self.operator().minimal_polynomial()
1055
1056 # If we don't orthonormalize the subalgebra's basis, then the
1057 # first two monomials in the subalgebra will be self^0 and
1058 # self^1... assuming that self^1 is not a scalar multiple of
1059 # self^0 (the unit element). We special case these to avoid
1060 # having to solve a system to coerce self into the subalgebra.
1061 A = self.subalgebra_generated_by(orthonormalize=False)
1062
1063 if A.dimension() == 1:
1064 # Does a solve to find the scalar multiple alpha such that
1065 # alpha*unit = self. We have to do this because the basis
1066 # for the subalgebra will be [ self^0 ], and not [ self^1 ]!
1067 unit = self.parent().one()
1068 alpha = self.to_vector() / unit.to_vector()
1069 return (unit.operator()*alpha).minimal_polynomial()
1070 else:
1071 # If the dimension of the subalgebra is >= 2, then we just
1072 # use the second basis element.
1073 return A.monomial(1).operator().minimal_polynomial()
1074
1075
1076
1077 def to_matrix(self):
1078 """
1079 Return an (often more natural) representation of this element as a
1080 matrix.
1081
1082 Every finite-dimensional Euclidean Jordan Algebra is a direct
1083 sum of five simple algebras, four of which comprise Hermitian
1084 matrices. This method returns a "natural" matrix
1085 representation of this element as either a Hermitian matrix or
1086 column vector.
1087
1088 SETUP::
1089
1090 sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1091 ....: HadamardEJA,
1092 ....: QuaternionHermitianEJA,
1093 ....: RealSymmetricEJA)
1094
1095 EXAMPLES::
1096
1097 sage: J = ComplexHermitianEJA(3)
1098 sage: J.one()
1099 b0 + b3 + b8
1100 sage: J.one().to_matrix()
1101 [1 0 0 0 0 0]
1102 [0 1 0 0 0 0]
1103 [0 0 1 0 0 0]
1104 [0 0 0 1 0 0]
1105 [0 0 0 0 1 0]
1106 [0 0 0 0 0 1]
1107
1108 ::
1109
1110 sage: J = QuaternionHermitianEJA(2)
1111 sage: J.one()
1112 b0 + b5
1113 sage: J.one().to_matrix()
1114 +---+---+
1115 | 1 | 0 |
1116 +---+---+
1117 | 0 | 1 |
1118 +---+---+
1119
1120 This also works in Cartesian product algebras::
1121
1122 sage: J1 = HadamardEJA(1)
1123 sage: J2 = RealSymmetricEJA(2)
1124 sage: J = cartesian_product([J1,J2])
1125 sage: x = sum(J.gens())
1126 sage: x.to_matrix()[0]
1127 [1]
1128 sage: x.to_matrix()[1]
1129 [ 1 0.7071067811865475?]
1130 [0.7071067811865475? 1]
1131
1132 """
1133 B = self.parent().matrix_basis()
1134 W = self.parent().matrix_space()
1135
1136 if hasattr(W, 'cartesian_factors'):
1137 # Aaaaand linear combinations don't work in Cartesian
1138 # product spaces, even though they provide a method with
1139 # that name. This is hidden behind an "if" because the
1140 # _scale() function is slow.
1141 pairs = zip(B, self.to_vector())
1142 return W.sum( _scale(b, alpha) for (b,alpha) in pairs )
1143 else:
1144 # This is just a manual "from_vector()", but of course
1145 # matrix spaces aren't vector spaces in sage, so they
1146 # don't have a from_vector() method.
1147 return W.linear_combination( zip(B, self.to_vector()) )
1148
1149
1150
1151 def norm(self):
1152 """
1153 The norm of this element with respect to :meth:`inner_product`.
1154
1155 SETUP::
1156
1157 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1158 ....: HadamardEJA)
1159
1160 EXAMPLES::
1161
1162 sage: J = HadamardEJA(2)
1163 sage: x = sum(J.gens())
1164 sage: x.norm()
1165 1.414213562373095?
1166
1167 ::
1168
1169 sage: J = JordanSpinEJA(4)
1170 sage: x = sum(J.gens())
1171 sage: x.norm()
1172 2
1173
1174 """
1175 return self.inner_product(self).sqrt()
1176
1177
1178 def operator(self):
1179 """
1180 Return the left-multiplication-by-this-element
1181 operator on the ambient algebra.
1182
1183 SETUP::
1184
1185 sage: from mjo.eja.eja_algebra import random_eja
1186
1187 TESTS::
1188
1189 sage: set_random_seed()
1190 sage: J = random_eja()
1191 sage: x,y = J.random_elements(2)
1192 sage: x.operator()(y) == x*y
1193 True
1194 sage: y.operator()(x) == x*y
1195 True
1196
1197 """
1198 P = self.parent()
1199 left_mult_by_self = lambda y: self*y
1200 L = P.module_morphism(function=left_mult_by_self, codomain=P)
1201 return FiniteDimensionalEJAOperator(P, P, L.matrix() )
1202
1203
1204 def quadratic_representation(self, other=None):
1205 """
1206 Return the quadratic representation of this element.
1207
1208 SETUP::
1209
1210 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1211 ....: random_eja)
1212
1213 EXAMPLES:
1214
1215 The explicit form in the spin factor algebra is given by
1216 Alizadeh's Example 11.12::
1217
1218 sage: set_random_seed()
1219 sage: x = JordanSpinEJA.random_instance().random_element()
1220 sage: x_vec = x.to_vector()
1221 sage: Q = matrix.identity(x.base_ring(), 0)
1222 sage: n = x_vec.degree()
1223 sage: if n > 0:
1224 ....: x0 = x_vec[0]
1225 ....: x_bar = x_vec[1:]
1226 ....: A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1227 ....: B = 2*x0*x_bar.row()
1228 ....: C = 2*x0*x_bar.column()
1229 ....: D = matrix.identity(x.base_ring(), n-1)
1230 ....: D = (x0^2 - x_bar.inner_product(x_bar))*D
1231 ....: D = D + 2*x_bar.tensor_product(x_bar)
1232 ....: Q = matrix.block(2,2,[A,B,C,D])
1233 sage: Q == x.quadratic_representation().matrix()
1234 True
1235
1236 Test all of the properties from Theorem 11.2 in Alizadeh::
1237
1238 sage: set_random_seed()
1239 sage: J = random_eja()
1240 sage: x,y = J.random_elements(2)
1241 sage: Lx = x.operator()
1242 sage: Lxx = (x*x).operator()
1243 sage: Qx = x.quadratic_representation()
1244 sage: Qy = y.quadratic_representation()
1245 sage: Qxy = x.quadratic_representation(y)
1246 sage: Qex = J.one().quadratic_representation(x)
1247 sage: n = ZZ.random_element(10)
1248 sage: Qxn = (x^n).quadratic_representation()
1249
1250 Property 1:
1251
1252 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1253 True
1254
1255 Property 2 (multiply on the right for :trac:`28272`):
1256
1257 sage: alpha = J.base_ring().random_element()
1258 sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1259 True
1260
1261 Property 3:
1262
1263 sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1264 True
1265
1266 sage: not x.is_invertible() or (
1267 ....: ~Qx
1268 ....: ==
1269 ....: x.inverse().quadratic_representation() )
1270 True
1271
1272 sage: Qxy(J.one()) == x*y
1273 True
1274
1275 Property 4:
1276
1277 sage: not x.is_invertible() or (
1278 ....: x.quadratic_representation(x.inverse())*Qx
1279 ....: == Qx*x.quadratic_representation(x.inverse()) )
1280 True
1281
1282 sage: not x.is_invertible() or (
1283 ....: x.quadratic_representation(x.inverse())*Qx
1284 ....: ==
1285 ....: 2*Lx*Qex - Qx )
1286 True
1287
1288 sage: 2*Lx*Qex - Qx == Lxx
1289 True
1290
1291 Property 5:
1292
1293 sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1294 True
1295
1296 Property 6:
1297
1298 sage: Qxn == (Qx)^n
1299 True
1300
1301 Property 7:
1302
1303 sage: not x.is_invertible() or (
1304 ....: Qx*x.inverse().operator() == Lx )
1305 True
1306
1307 Property 8:
1308
1309 sage: not x.operator_commutes_with(y) or (
1310 ....: Qx(y)^n == Qxn(y^n) )
1311 True
1312
1313 """
1314 if other is None:
1315 other=self
1316 elif not other in self.parent():
1317 raise TypeError("'other' must live in the same algebra")
1318
1319 L = self.operator()
1320 M = other.operator()
1321 return ( L*M + M*L - (self*other).operator() )
1322
1323
1324
1325 def spectral_decomposition(self):
1326 """
1327 Return the unique spectral decomposition of this element.
1328
1329 ALGORITHM:
1330
1331 Following Faraut and Korányi's Theorem III.1.1, we restrict this
1332 element's left-multiplication-by operator to the subalgebra it
1333 generates. We then compute the spectral decomposition of that
1334 operator, and the spectral projectors we get back must be the
1335 left-multiplication-by operators for the idempotents we
1336 seek. Thus applying them to the identity element gives us those
1337 idempotents.
1338
1339 Since the eigenvalues are required to be distinct, we take
1340 the spectral decomposition of the zero element to be zero
1341 times the identity element of the algebra (which is idempotent,
1342 obviously).
1343
1344 SETUP::
1345
1346 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1347
1348 EXAMPLES:
1349
1350 The spectral decomposition of the identity is ``1`` times itself,
1351 and the spectral decomposition of zero is ``0`` times the identity::
1352
1353 sage: J = RealSymmetricEJA(3)
1354 sage: J.one()
1355 b0 + b2 + b5
1356 sage: J.one().spectral_decomposition()
1357 [(1, b0 + b2 + b5)]
1358 sage: J.zero().spectral_decomposition()
1359 [(0, b0 + b2 + b5)]
1360
1361 TESTS::
1362
1363 sage: J = RealSymmetricEJA(4)
1364 sage: x = sum(J.gens())
1365 sage: sd = x.spectral_decomposition()
1366 sage: l0 = sd[0][0]
1367 sage: l1 = sd[1][0]
1368 sage: c0 = sd[0][1]
1369 sage: c1 = sd[1][1]
1370 sage: c0.inner_product(c1) == 0
1371 True
1372 sage: c0.is_idempotent()
1373 True
1374 sage: c1.is_idempotent()
1375 True
1376 sage: c0 + c1 == J.one()
1377 True
1378 sage: l0*c0 + l1*c1 == x
1379 True
1380
1381 The spectral decomposition should work in subalgebras, too::
1382
1383 sage: J = RealSymmetricEJA(4)
1384 sage: (b0, b1, b2, b3, b4, b5, b6, b7, b8, b9) = J.gens()
1385 sage: A = 2*b5 - 2*b8
1386 sage: (lambda1, c1) = A.spectral_decomposition()[1]
1387 sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1388 sage: (f0, f1, f2) = J1.gens()
1389 sage: f0.spectral_decomposition()
1390 [(0, c2), (1, c0)]
1391
1392 """
1393 A = self.subalgebra_generated_by(orthonormalize=True)
1394 result = []
1395 for (evalue, proj) in A(self).operator().spectral_decomposition():
1396 result.append( (evalue, proj(A.one()).superalgebra_element()) )
1397 return result
1398
1399 def subalgebra_generated_by(self, **kwargs):
1400 """
1401 Return the associative subalgebra of the parent EJA generated
1402 by this element.
1403
1404 Since our parent algebra is unital, we want "subalgebra" to mean
1405 "unital subalgebra" as well; thus the subalgebra that an element
1406 generates will itself be a Euclidean Jordan algebra after
1407 restricting the algebra operations appropriately. This is the
1408 subalgebra that Faraut and Korányi work with in section II.2, for
1409 example.
1410
1411 SETUP::
1412
1413 sage: from mjo.eja.eja_algebra import (random_eja,
1414 ....: HadamardEJA,
1415 ....: RealSymmetricEJA)
1416
1417 EXAMPLES:
1418
1419 We can create subalgebras of Cartesian product EJAs that are not
1420 themselves Cartesian product EJAs (they're just "regular" EJAs)::
1421
1422 sage: J1 = HadamardEJA(3)
1423 sage: J2 = RealSymmetricEJA(2)
1424 sage: J = cartesian_product([J1,J2])
1425 sage: J.one().subalgebra_generated_by()
1426 Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
1427
1428 TESTS:
1429
1430 This subalgebra, being composed of only powers, is associative::
1431
1432 sage: set_random_seed()
1433 sage: x0 = random_eja().random_element()
1434 sage: A = x0.subalgebra_generated_by()
1435 sage: x,y,z = A.random_elements(3)
1436 sage: (x*y)*z == x*(y*z)
1437 True
1438
1439 Squaring in the subalgebra should work the same as in
1440 the superalgebra::
1441
1442 sage: set_random_seed()
1443 sage: x = random_eja().random_element()
1444 sage: A = x.subalgebra_generated_by()
1445 sage: A(x^2) == A(x)*A(x)
1446 True
1447
1448 By definition, the subalgebra generated by the zero element is
1449 the one-dimensional algebra generated by the identity
1450 element... unless the original algebra was trivial, in which
1451 case the subalgebra is trivial too::
1452
1453 sage: set_random_seed()
1454 sage: A = random_eja().zero().subalgebra_generated_by()
1455 sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1456 True
1457
1458 """
1459 powers = tuple( self**k for k in range(self.degree()) )
1460 A = self.parent().subalgebra(powers,
1461 associative=True,
1462 check_field=False,
1463 check_axioms=False,
1464 **kwargs)
1465 A.one.set_cache(A(self.parent().one()))
1466 return A
1467
1468
1469 def subalgebra_idempotent(self):
1470 """
1471 Find an idempotent in the associative subalgebra I generate
1472 using Proposition 2.3.5 in Baes.
1473
1474 SETUP::
1475
1476 sage: from mjo.eja.eja_algebra import random_eja
1477
1478 TESTS:
1479
1480 Ensure that we can find an idempotent in a non-trivial algebra
1481 where there are non-nilpotent elements, or that we get the dumb
1482 solution in the trivial algebra::
1483
1484 sage: set_random_seed()
1485 sage: J = random_eja()
1486 sage: x = J.random_element()
1487 sage: while x.is_nilpotent() and not J.is_trivial():
1488 ....: x = J.random_element()
1489 sage: c = x.subalgebra_idempotent()
1490 sage: c^2 == c
1491 True
1492
1493 """
1494 if self.parent().is_trivial():
1495 return self
1496
1497 if self.is_nilpotent():
1498 raise ValueError("this only works with non-nilpotent elements!")
1499
1500 J = self.subalgebra_generated_by()
1501 u = J(self)
1502
1503 # The image of the matrix of left-u^m-multiplication
1504 # will be minimal for some natural number s...
1505 s = 0
1506 minimal_dim = J.dimension()
1507 for i in range(1, minimal_dim):
1508 this_dim = (u**i).operator().matrix().image().dimension()
1509 if this_dim < minimal_dim:
1510 minimal_dim = this_dim
1511 s = i
1512
1513 # Now minimal_matrix should correspond to the smallest
1514 # non-zero subspace in Baes's (or really, Koecher's)
1515 # proposition.
1516 #
1517 # However, we need to restrict the matrix to work on the
1518 # subspace... or do we? Can't we just solve, knowing that
1519 # A(c) = u^(s+1) should have a solution in the big space,
1520 # too?
1521 #
1522 # Beware, solve_right() means that we're using COLUMN vectors.
1523 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1524 u_next = u**(s+1)
1525 A = u_next.operator().matrix()
1526 c = J.from_vector(A.solve_right(u_next.to_vector()))
1527
1528 # Now c is the idempotent we want, but it still lives in the subalgebra.
1529 return c.superalgebra_element()
1530
1531
1532 def trace(self):
1533 """
1534 Return my trace, the sum of my eigenvalues.
1535
1536 In a trivial algebra, however you want to look at it, the trace is
1537 an empty sum for which we declare the result to be zero.
1538
1539 SETUP::
1540
1541 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1542 ....: HadamardEJA,
1543 ....: TrivialEJA,
1544 ....: random_eja)
1545
1546 EXAMPLES::
1547
1548 sage: J = TrivialEJA()
1549 sage: J.zero().trace()
1550 0
1551
1552 ::
1553 sage: J = JordanSpinEJA(3)
1554 sage: x = sum(J.gens())
1555 sage: x.trace()
1556 2
1557
1558 ::
1559
1560 sage: J = HadamardEJA(5)
1561 sage: J.one().trace()
1562 5
1563
1564 TESTS:
1565
1566 The trace of an element is a real number::
1567
1568 sage: set_random_seed()
1569 sage: J = random_eja()
1570 sage: J.random_element().trace() in RLF
1571 True
1572
1573 The trace is linear::
1574
1575 sage: set_random_seed()
1576 sage: J = random_eja()
1577 sage: x,y = J.random_elements(2)
1578 sage: alpha = J.base_ring().random_element()
1579 sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1580 True
1581
1582 """
1583 P = self.parent()
1584 r = P.rank()
1585
1586 if r == 0:
1587 # Special case for the trivial algebra where
1588 # the trace is an empty sum.
1589 return P.base_ring().zero()
1590
1591 p = P._charpoly_coefficients()[r-1]
1592 # The _charpoly_coeff function already adds the factor of
1593 # -1 to ensure that _charpoly_coeff(r-1) is really what
1594 # appears in front of t^{r-1} in the charpoly. However,
1595 # we want the negative of THAT for the trace.
1596 return -p(*self.to_vector())
1597
1598
1599 def trace_inner_product(self, other):
1600 """
1601 Return the trace inner product of myself and ``other``.
1602
1603 SETUP::
1604
1605 sage: from mjo.eja.eja_algebra import random_eja
1606
1607 TESTS:
1608
1609 The trace inner product is commutative, bilinear, and associative::
1610
1611 sage: set_random_seed()
1612 sage: J = random_eja()
1613 sage: x,y,z = J.random_elements(3)
1614 sage: # commutative
1615 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1616 True
1617 sage: # bilinear
1618 sage: a = J.base_ring().random_element();
1619 sage: actual = (a*(x+z)).trace_inner_product(y)
1620 sage: expected = ( a*x.trace_inner_product(y) +
1621 ....: a*z.trace_inner_product(y) )
1622 sage: actual == expected
1623 True
1624 sage: actual = x.trace_inner_product(a*(y+z))
1625 sage: expected = ( a*x.trace_inner_product(y) +
1626 ....: a*x.trace_inner_product(z) )
1627 sage: actual == expected
1628 True
1629 sage: # associative
1630 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1631 True
1632
1633 """
1634 if not other in self.parent():
1635 raise TypeError("'other' must live in the same algebra")
1636
1637 return (self*other).trace()
1638
1639
1640 def trace_norm(self):
1641 """
1642 The norm of this element with respect to :meth:`trace_inner_product`.
1643
1644 SETUP::
1645
1646 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1647 ....: HadamardEJA)
1648
1649 EXAMPLES::
1650
1651 sage: J = HadamardEJA(2)
1652 sage: x = sum(J.gens())
1653 sage: x.trace_norm()
1654 1.414213562373095?
1655
1656 ::
1657
1658 sage: J = JordanSpinEJA(4)
1659 sage: x = sum(J.gens())
1660 sage: x.trace_norm()
1661 2.828427124746190?
1662
1663 """
1664 return self.trace_inner_product(self).sqrt()
My personal library of Sage code.