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
   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: e0,e1,e2 = J.gens()
 171             sage: A = (e0 + 2*e1 + 3*e2).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: e0,e1 = J.gens()
 352             sage: x = sum( J.gens() )
 353             sage: x.det()
 354             0
 355
 356         ::
 357
 358             sage: J = JordanSpinEJA(3)
 359             sage: e0,e1,e2 = J.gens()
 360             sage: x = sum( J.gens() )
 361             sage: x.det()
 362             -1
 363
 364         The determinant of the sole element in the rank-zero trivial
 365         algebra is ``1``, by three paths of reasoning. First, its
 366         characteristic polynomial is a constant ``1``, so the constant
 367         term in that polynomial is ``1``. Second, the characteristic
 368         polynomial evaluated at zero is again ``1``. And finally, the
 369         (empty) product of its eigenvalues is likewise just unity::
 370
 371             sage: J = TrivialEJA()
 372             sage: J.zero().det()
 373             1
 374
 375         TESTS:
 376
 377         An element is invertible if and only if its determinant is
 378         non-zero::
 379
 380             sage: set_random_seed()
 381             sage: x = random_eja().random_element()
 382             sage: x.is_invertible() == (x.det() != 0)
 383             True
 384
 385         Ensure that the determinant is multiplicative on an associative
 386         subalgebra as in Faraut and Korányi's Proposition II.2.2::
 387
 388             sage: set_random_seed()
 389             sage: J = random_eja().random_element().subalgebra_generated_by()
 390             sage: x,y = J.random_elements(2)
 391             sage: (x*y).det() == x.det()*y.det()
 392             True
 393
 394         The determinant in matrix algebras is just the usual determinant::
 395
 396             sage: set_random_seed()
 397             sage: X = matrix.random(QQ,3)
 398             sage: X = X + X.T
 399             sage: J1 = RealSymmetricEJA(3)
 400             sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
 401             sage: expected = X.det()
 402             sage: actual1 = J1(X).det()
 403             sage: actual2 = J2(X).det()
 404             sage: actual1 == expected
 405             True
 406             sage: actual2 == expected
 407             True
 408
 409         ::
 410
 411             sage: set_random_seed()
 412             sage: J1 = ComplexHermitianEJA(2)
 413             sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
 414             sage: X = matrix.random(GaussianIntegers(), 2)
 415             sage: X = X + X.H
 416             sage: expected = AA(X.det())
 417             sage: actual1 = J1(J1.real_embed(X)).det()
 418             sage: actual2 = J2(J2.real_embed(X)).det()
 419             sage: expected == actual1
 420             True
 421             sage: expected == actual2
 422             True
 423
 424         """
 425         P = self.parent()
 426         r = P.rank()
 427
 428         if r == 0:
 429             # Special case, since we don't get the a0=1
 430             # coefficient when the rank of the algebra
 431             # is zero.
 432             return P.base_ring().one()
 433
 434         p = P._charpoly_coefficients()[0]
 435         # The _charpoly_coeff function already adds the factor of -1
 436         # to ensure that _charpoly_coefficients()[0] is really what
 437         # appears in front of t^{0} in the charpoly. However, we want
 438         # (-1)^r times THAT for the determinant.
 439         return ((-1)**r)*p(*self.to_vector())
 440
 441
 442     @cached_method
 443     def inverse(self):
 444         """
 445         Return the Jordan-multiplicative inverse of this element.
 446
 447         ALGORITHM:
 448
 449         In general we appeal to the quadratic representation as in
 450         Koecher's Theorem 12 in Chapter III, Section 5. But if the
 451         parent algebra's "characteristic polynomial of" coefficients
 452         happen to be cached, then we use Proposition II.2.4 in Faraut
 453         and Korányi which gives a formula for the inverse based on the
 454         characteristic polynomial and the Cayley-Hamilton theorem for
 455         Euclidean Jordan algebras::
 456
 457         SETUP::
 458
 459             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 460             ....:                                  JordanSpinEJA,
 461             ....:                                  random_eja)
 462
 463         EXAMPLES:
 464
 465         The inverse in the spin factor algebra is given in Alizadeh's
 466         Example 11.11::
 467
 468             sage: set_random_seed()
 469             sage: J = JordanSpinEJA.random_instance()
 470             sage: x = J.random_element()
 471             sage: while not x.is_invertible():
 472             ....:     x = J.random_element()
 473             sage: x_vec = x.to_vector()
 474             sage: x0 = x_vec[:1]
 475             sage: x_bar = x_vec[1:]
 476             sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
 477             sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
 478             sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
 479             sage: x.inverse() == J.from_vector(x_inverse)
 480             True
 481
 482         Trying to invert a non-invertible element throws an error:
 483
 484             sage: JordanSpinEJA(3).zero().inverse()
 485             Traceback (most recent call last):
 486             ...
 487             ZeroDivisionError: element is not invertible
 488
 489         TESTS:
 490
 491         The identity element is its own inverse::
 492
 493             sage: set_random_seed()
 494             sage: J = random_eja()
 495             sage: J.one().inverse() == J.one()
 496             True
 497
 498         If an element has an inverse, it acts like one::
 499
 500             sage: set_random_seed()
 501             sage: J = random_eja()
 502             sage: x = J.random_element()
 503             sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
 504             True
 505
 506         The inverse of the inverse is what we started with::
 507
 508             sage: set_random_seed()
 509             sage: J = random_eja()
 510             sage: x = J.random_element()
 511             sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
 512             True
 513
 514         Proposition II.2.3 in Faraut and Korányi says that the inverse
 515         of an element is the inverse of its left-multiplication operator
 516         applied to the algebra's identity, when that inverse exists::
 517
 518             sage: set_random_seed()
 519             sage: J = random_eja()
 520             sage: x = J.random_element()
 521             sage: (not x.operator().is_invertible()) or (
 522             ....:    x.operator().inverse()(J.one()) == x.inverse() )
 523             True
 524
 525         Check that the fast (cached) and slow algorithms give the same
 526         answer::
 527
 528             sage: set_random_seed()                              # long time
 529             sage: J = random_eja(field=QQ, orthonormalize=False) # long time
 530             sage: x = J.random_element()                         # long time
 531             sage: while not x.is_invertible():                   # long time
 532             ....:     x = J.random_element()                     # long time
 533             sage: slow = x.inverse()                             # long time
 534             sage: _ = J._charpoly_coefficients()                 # long time
 535             sage: fast = x.inverse()                             # long time
 536             sage: slow == fast                                   # long time
 537             True
 538         """
 539         not_invertible_msg = "element is not invertible"
 540         if self.parent()._charpoly_coefficients.is_in_cache():
 541             # We can invert using our charpoly if it will be fast to
 542             # compute. If the coefficients are cached, our rank had
 543             # better be too!
 544             if self.det().is_zero():
 545                 raise ZeroDivisionError(not_invertible_msg)
 546             r = self.parent().rank()
 547             a = self.characteristic_polynomial().coefficients(sparse=False)
 548             return (-1)**(r+1)*sum(a[i+1]*self**i for i in range(r))/self.det()
 549
 550         try:
 551             inv = (~self.quadratic_representation())(self)
 552             self.is_invertible.set_cache(True)
 553             return inv
 554         except ZeroDivisionError:
 555             self.is_invertible.set_cache(False)
 556             raise ZeroDivisionError(not_invertible_msg)
 557
 558
 559     @cached_method
 560     def is_invertible(self):
 561         """
 562         Return whether or not this element is invertible.
 563
 564         ALGORITHM:
 565
 566         If computing my determinant will be fast, we do so and compare
 567         with zero (Proposition II.2.4 in Faraut and
 568         Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
 569         reduces the problem to the invertibility of my quadratic
 570         representation.
 571
 572         SETUP::
 573
 574             sage: from mjo.eja.eja_algebra import random_eja
 575
 576         TESTS:
 577
 578         The identity element is always invertible::
 579
 580             sage: set_random_seed()
 581             sage: J = random_eja()
 582             sage: J.one().is_invertible()
 583             True
 584
 585         The zero element is never invertible in a non-trivial algebra::
 586
 587             sage: set_random_seed()
 588             sage: J = random_eja()
 589             sage: (not J.is_trivial()) and J.zero().is_invertible()
 590             False
 591
 592         Test that the fast (cached) and slow algorithms give the same
 593         answer::
 594
 595             sage: set_random_seed()                              # long time
 596             sage: J = random_eja(field=QQ, orthonormalize=False) # long time
 597             sage: x = J.random_element()                         # long time
 598             sage: slow = x.is_invertible()                       # long time
 599             sage: _ = J._charpoly_coefficients()                 # long time
 600             sage: fast = x.is_invertible()                       # long time
 601             sage: slow == fast                                   # long time
 602             True
 603         """
 604         if self.is_zero():
 605             if self.parent().is_trivial():
 606                 return True
 607             else:
 608                 return False
 609
 610         if self.parent()._charpoly_coefficients.is_in_cache():
 611             # The determinant will be quicker than inverting the
 612             # quadratic representation, most likely.
 613             return (not self.det().is_zero())
 614
 615         # The easiest way to determine if I'm invertible is to try.
 616         try:
 617             inv = (~self.quadratic_representation())(self)
 618             self.inverse.set_cache(inv)
 619             return True
 620         except ZeroDivisionError:
 621             return False
 622
 623
 624     def is_primitive_idempotent(self):
 625         """
 626         Return whether or not this element is a primitive (or minimal)
 627         idempotent.
 628
 629         A primitive idempotent is a non-zero idempotent that is not
 630         the sum of two other non-zero idempotents. Remark 2.7.15 in
 631         Baes shows that this is what he refers to as a "minimal
 632         idempotent."
 633
 634         An element of a Euclidean Jordan algebra is a minimal idempotent
 635         if it :meth:`is_idempotent` and if its Peirce subalgebra
 636         corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
 637         Proposition 2.7.17).
 638
 639         SETUP::
 640
 641             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 642             ....:                                  RealSymmetricEJA,
 643             ....:                                  TrivialEJA,
 644             ....:                                  random_eja)
 645
 646         WARNING::
 647
 648         This method is sloooooow.
 649
 650         EXAMPLES:
 651
 652         The spectral decomposition of a non-regular element should always
 653         contain at least one non-minimal idempotent::
 654
 655             sage: J = RealSymmetricEJA(3)
 656             sage: x = sum(J.gens())
 657             sage: x.is_regular()
 658             False
 659             sage: [ c.is_primitive_idempotent()
 660             ....:   for (l,c) in x.spectral_decomposition() ]
 661             [False, True]
 662
 663         On the other hand, the spectral decomposition of a regular
 664         element should always be in terms of minimal idempotents::
 665
 666             sage: J = JordanSpinEJA(4)
 667             sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
 668             sage: x.is_regular()
 669             True
 670             sage: [ c.is_primitive_idempotent()
 671             ....:   for (l,c) in x.spectral_decomposition() ]
 672             [True, True]
 673
 674         TESTS:
 675
 676         The identity element is minimal only in an EJA of rank one::
 677
 678             sage: set_random_seed()
 679             sage: J = random_eja()
 680             sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
 681             True
 682
 683         A non-idempotent cannot be a minimal idempotent::
 684
 685             sage: set_random_seed()
 686             sage: J = JordanSpinEJA(4)
 687             sage: x = J.random_element()
 688             sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
 689             False
 690
 691         Proposition 2.7.19 in Baes says that an element is a minimal
 692         idempotent if and only if it's idempotent with trace equal to
 693         unity::
 694
 695             sage: set_random_seed()
 696             sage: J = JordanSpinEJA(4)
 697             sage: x = J.random_element()
 698             sage: expected = (x.is_idempotent() and x.trace() == 1)
 699             sage: actual = x.is_primitive_idempotent()
 700             sage: actual == expected
 701             True
 702
 703         Primitive idempotents must be non-zero::
 704
 705             sage: set_random_seed()
 706             sage: J = random_eja()
 707             sage: J.zero().is_idempotent()
 708             True
 709             sage: J.zero().is_primitive_idempotent()
 710             False
 711
 712         As a consequence of the fact that primitive idempotents must
 713         be non-zero, there are no primitive idempotents in a trivial
 714         Euclidean Jordan algebra::
 715
 716             sage: J = TrivialEJA()
 717             sage: J.one().is_idempotent()
 718             True
 719             sage: J.one().is_primitive_idempotent()
 720             False
 721
 722         """
 723         if not self.is_idempotent():
 724             return False
 725
 726         if self.is_zero():
 727             return False
 728
 729         (_,_,J1) = self.parent().peirce_decomposition(self)
 730         return (J1.dimension() == 1)
 731
 732
 733     def is_nilpotent(self):
 734         """
 735         Return whether or not some power of this element is zero.
 736
 737         ALGORITHM:
 738
 739         We use Theorem 5 in Chapter III of Koecher, which says that
 740         an element ``x`` is nilpotent if and only if ``x.operator()``
 741         is nilpotent. And it is a basic fact of linear algebra that
 742         an operator on an `n`-dimensional space is nilpotent if and
 743         only if, when raised to the `n`th power, it equals the zero
 744         operator (for example, see Axler Corollary 8.8).
 745
 746         SETUP::
 747
 748             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 749             ....:                                  random_eja)
 750
 751         EXAMPLES::
 752
 753             sage: J = JordanSpinEJA(3)
 754             sage: x = sum(J.gens())
 755             sage: x.is_nilpotent()
 756             False
 757
 758         TESTS:
 759
 760         The identity element is never nilpotent, except in a trivial EJA::
 761
 762             sage: set_random_seed()
 763             sage: J = random_eja()
 764             sage: J.one().is_nilpotent() and not J.is_trivial()
 765             False
 766
 767         The additive identity is always nilpotent::
 768
 769             sage: set_random_seed()
 770             sage: random_eja().zero().is_nilpotent()
 771             True
 772
 773         """
 774         P = self.parent()
 775         zero_operator = P.zero().operator()
 776         return self.operator()**P.dimension() == zero_operator
 777
 778
 779     def is_regular(self):
 780         """
 781         Return whether or not this is a regular element.
 782
 783         SETUP::
 784
 785             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 786             ....:                                  random_eja)
 787
 788         EXAMPLES:
 789
 790         The identity element always has degree one, but any element
 791         linearly-independent from it is regular::
 792
 793             sage: J = JordanSpinEJA(5)
 794             sage: J.one().is_regular()
 795             False
 796             sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 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: e0,e1,e2,e3 = J.gens()
 847             sage: (e0 - e1).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             # We would generate a zero-dimensional subalgebra
1051             # where the minimal polynomial would be constant.
1052             # That might be correct, but only if *this* algebra
1053             # is trivial too.
1054             if not self.parent().is_trivial():
1055                 # Pretty sure we know what the minimal polynomial of
1056                 # the zero operator is going to be. This ensures
1057                 # consistency of e.g. the polynomial variable returned
1058                 # in the "normal" case without us having to think about it.
1059                 return self.operator().minimal_polynomial()
1060
1061         A = self.subalgebra_generated_by(orthonormalize=False)
1062         return A(self).operator().minimal_polynomial()
1063
1064
1065
1066     def to_matrix(self):
1067         """
1068         Return an (often more natural) representation of this element as a
1069         matrix.
1070
1071         Every finite-dimensional Euclidean Jordan Algebra is a direct
1072         sum of five simple algebras, four of which comprise Hermitian
1073         matrices. This method returns a "natural" matrix
1074         representation of this element as either a Hermitian matrix or
1075         column vector.
1076
1077         SETUP::
1078
1079             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1080             ....:                                  QuaternionHermitianEJA)
1081
1082         EXAMPLES::
1083
1084             sage: J = ComplexHermitianEJA(3)
1085             sage: J.one()
1086             e0 + e3 + e8
1087             sage: J.one().to_matrix()
1088             [1 0 0 0 0 0]
1089             [0 1 0 0 0 0]
1090             [0 0 1 0 0 0]
1091             [0 0 0 1 0 0]
1092             [0 0 0 0 1 0]
1093             [0 0 0 0 0 1]
1094
1095         ::
1096
1097             sage: J = QuaternionHermitianEJA(2)
1098             sage: J.one()
1099             e0 + e5
1100             sage: J.one().to_matrix()
1101             [1 0 0 0 0 0 0 0]
1102             [0 1 0 0 0 0 0 0]
1103             [0 0 1 0 0 0 0 0]
1104             [0 0 0 1 0 0 0 0]
1105             [0 0 0 0 1 0 0 0]
1106             [0 0 0 0 0 1 0 0]
1107             [0 0 0 0 0 0 1 0]
1108             [0 0 0 0 0 0 0 1]
1109
1110         """
1111         B = self.parent().matrix_basis()
1112         W = self.parent().matrix_space()
1113
1114         # This is just a manual "from_vector()", but of course
1115         # matrix spaces aren't vector spaces in sage, so they
1116         # don't have a from_vector() method.
1117         return W.linear_combination( zip(B, self.to_vector()) )
1118
1119
1120
1121     def norm(self):
1122         """
1123         The norm of this element with respect to :meth:`inner_product`.
1124
1125         SETUP::
1126
1127             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1128             ....:                                  HadamardEJA)
1129
1130         EXAMPLES::
1131
1132             sage: J = HadamardEJA(2)
1133             sage: x = sum(J.gens())
1134             sage: x.norm()
1135             1.414213562373095?
1136
1137         ::
1138
1139             sage: J = JordanSpinEJA(4)
1140             sage: x = sum(J.gens())
1141             sage: x.norm()
1142             2
1143
1144         """
1145         return self.inner_product(self).sqrt()
1146
1147
1148     def operator(self):
1149         """
1150         Return the left-multiplication-by-this-element
1151         operator on the ambient algebra.
1152
1153         SETUP::
1154
1155             sage: from mjo.eja.eja_algebra import random_eja
1156
1157         TESTS::
1158
1159             sage: set_random_seed()
1160             sage: J = random_eja()
1161             sage: x,y = J.random_elements(2)
1162             sage: x.operator()(y) == x*y
1163             True
1164             sage: y.operator()(x) == x*y
1165             True
1166
1167         """
1168         P = self.parent()
1169         left_mult_by_self = lambda y: self*y
1170         L = P.module_morphism(function=left_mult_by_self, codomain=P)
1171         return FiniteDimensionalEJAOperator(P, P, L.matrix() )
1172
1173
1174     def quadratic_representation(self, other=None):
1175         """
1176         Return the quadratic representation of this element.
1177
1178         SETUP::
1179
1180             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1181             ....:                                  random_eja)
1182
1183         EXAMPLES:
1184
1185         The explicit form in the spin factor algebra is given by
1186         Alizadeh's Example 11.12::
1187
1188             sage: set_random_seed()
1189             sage: x = JordanSpinEJA.random_instance().random_element()
1190             sage: x_vec = x.to_vector()
1191             sage: Q = matrix.identity(x.base_ring(), 0)
1192             sage: n = x_vec.degree()
1193             sage: if n > 0:
1194             ....:     x0 = x_vec[0]
1195             ....:     x_bar = x_vec[1:]
1196             ....:     A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1197             ....:     B = 2*x0*x_bar.row()
1198             ....:     C = 2*x0*x_bar.column()
1199             ....:     D = matrix.identity(x.base_ring(), n-1)
1200             ....:     D = (x0^2 - x_bar.inner_product(x_bar))*D
1201             ....:     D = D + 2*x_bar.tensor_product(x_bar)
1202             ....:     Q = matrix.block(2,2,[A,B,C,D])
1203             sage: Q == x.quadratic_representation().matrix()
1204             True
1205
1206         Test all of the properties from Theorem 11.2 in Alizadeh::
1207
1208             sage: set_random_seed()
1209             sage: J = random_eja()
1210             sage: x,y = J.random_elements(2)
1211             sage: Lx = x.operator()
1212             sage: Lxx = (x*x).operator()
1213             sage: Qx = x.quadratic_representation()
1214             sage: Qy = y.quadratic_representation()
1215             sage: Qxy = x.quadratic_representation(y)
1216             sage: Qex = J.one().quadratic_representation(x)
1217             sage: n = ZZ.random_element(10)
1218             sage: Qxn = (x^n).quadratic_representation()
1219
1220         Property 1:
1221
1222             sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1223             True
1224
1225         Property 2 (multiply on the right for :trac:`28272`):
1226
1227             sage: alpha = J.base_ring().random_element()
1228             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1229             True
1230
1231         Property 3:
1232
1233             sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1234             True
1235
1236             sage: not x.is_invertible() or (
1237             ....:   ~Qx
1238             ....:   ==
1239             ....:   x.inverse().quadratic_representation() )
1240             True
1241
1242             sage: Qxy(J.one()) == x*y
1243             True
1244
1245         Property 4:
1246
1247             sage: not x.is_invertible() or (
1248             ....:   x.quadratic_representation(x.inverse())*Qx
1249             ....:   == Qx*x.quadratic_representation(x.inverse()) )
1250             True
1251
1252             sage: not x.is_invertible() or (
1253             ....:   x.quadratic_representation(x.inverse())*Qx
1254             ....:   ==
1255             ....:   2*Lx*Qex - Qx )
1256             True
1257
1258             sage: 2*Lx*Qex - Qx == Lxx
1259             True
1260
1261         Property 5:
1262
1263             sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1264             True
1265
1266         Property 6:
1267
1268             sage: Qxn == (Qx)^n
1269             True
1270
1271         Property 7:
1272
1273             sage: not x.is_invertible() or (
1274             ....:   Qx*x.inverse().operator() == Lx )
1275             True
1276
1277         Property 8:
1278
1279             sage: not x.operator_commutes_with(y) or (
1280             ....:   Qx(y)^n == Qxn(y^n) )
1281             True
1282
1283         """
1284         if other is None:
1285             other=self
1286         elif not other in self.parent():
1287             raise TypeError("'other' must live in the same algebra")
1288
1289         L = self.operator()
1290         M = other.operator()
1291         return ( L*M + M*L - (self*other).operator() )
1292
1293
1294
1295     def spectral_decomposition(self):
1296         """
1297         Return the unique spectral decomposition of this element.
1298
1299         ALGORITHM:
1300
1301         Following Faraut and Korányi's Theorem III.1.1, we restrict this
1302         element's left-multiplication-by operator to the subalgebra it
1303         generates. We then compute the spectral decomposition of that
1304         operator, and the spectral projectors we get back must be the
1305         left-multiplication-by operators for the idempotents we
1306         seek. Thus applying them to the identity element gives us those
1307         idempotents.
1308
1309         Since the eigenvalues are required to be distinct, we take
1310         the spectral decomposition of the zero element to be zero
1311         times the identity element of the algebra (which is idempotent,
1312         obviously).
1313
1314         SETUP::
1315
1316             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1317
1318         EXAMPLES:
1319
1320         The spectral decomposition of the identity is ``1`` times itself,
1321         and the spectral decomposition of zero is ``0`` times the identity::
1322
1323             sage: J = RealSymmetricEJA(3)
1324             sage: J.one()
1325             e0 + e2 + e5
1326             sage: J.one().spectral_decomposition()
1327             [(1, e0 + e2 + e5)]
1328             sage: J.zero().spectral_decomposition()
1329             [(0, e0 + e2 + e5)]
1330
1331         TESTS::
1332
1333             sage: J = RealSymmetricEJA(4)
1334             sage: x = sum(J.gens())
1335             sage: sd = x.spectral_decomposition()
1336             sage: l0 = sd[0][0]
1337             sage: l1 = sd[1][0]
1338             sage: c0 = sd[0][1]
1339             sage: c1 = sd[1][1]
1340             sage: c0.inner_product(c1) == 0
1341             True
1342             sage: c0.is_idempotent()
1343             True
1344             sage: c1.is_idempotent()
1345             True
1346             sage: c0 + c1 == J.one()
1347             True
1348             sage: l0*c0 + l1*c1 == x
1349             True
1350
1351         The spectral decomposition should work in subalgebras, too::
1352
1353             sage: J = RealSymmetricEJA(4)
1354             sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens()
1355             sage: A = 2*e5 - 2*e8
1356             sage: (lambda1, c1) = A.spectral_decomposition()[1]
1357             sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1358             sage: (f0, f1, f2) = J1.gens()
1359             sage: f0.spectral_decomposition()
1360             [(0, f2), (1, f0)]
1361
1362         """
1363         A = self.subalgebra_generated_by(orthonormalize=True)
1364         result = []
1365         for (evalue, proj) in A(self).operator().spectral_decomposition():
1366             result.append( (evalue, proj(A.one()).superalgebra_element()) )
1367         return result
1368
1369     def subalgebra_generated_by(self, **kwargs):
1370         """
1371         Return the associative subalgebra of the parent EJA generated
1372         by this element.
1373
1374         Since our parent algebra is unital, we want "subalgebra" to mean
1375         "unital subalgebra" as well; thus the subalgebra that an element
1376         generates will itself be a Euclidean Jordan algebra after
1377         restricting the algebra operations appropriately. This is the
1378         subalgebra that Faraut and Korányi work with in section II.2, for
1379         example.
1380
1381         SETUP::
1382
1383             sage: from mjo.eja.eja_algebra import random_eja
1384
1385         TESTS:
1386
1387         This subalgebra, being composed of only powers, is associative::
1388
1389             sage: set_random_seed()
1390             sage: x0 = random_eja().random_element()
1391             sage: A = x0.subalgebra_generated_by()
1392             sage: x,y,z = A.random_elements(3)
1393             sage: (x*y)*z == x*(y*z)
1394             True
1395
1396         Squaring in the subalgebra should work the same as in
1397         the superalgebra::
1398
1399             sage: set_random_seed()
1400             sage: x = random_eja().random_element()
1401             sage: A = x.subalgebra_generated_by()
1402             sage: A(x^2) == A(x)*A(x)
1403             True
1404
1405         By definition, the subalgebra generated by the zero element is
1406         the one-dimensional algebra generated by the identity
1407         element... unless the original algebra was trivial, in which
1408         case the subalgebra is trivial too::
1409
1410             sage: set_random_seed()
1411             sage: A = random_eja().zero().subalgebra_generated_by()
1412             sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1413             True
1414
1415         """
1416         powers = tuple( self**k for k in range(self.degree()) )
1417         A = self.parent().subalgebra(powers, associative=True, **kwargs)
1418         A.one.set_cache(A(self.parent().one()))
1419         return A
1420
1421
1422     def subalgebra_idempotent(self):
1423         """
1424         Find an idempotent in the associative subalgebra I generate
1425         using Proposition 2.3.5 in Baes.
1426
1427         SETUP::
1428
1429             sage: from mjo.eja.eja_algebra import random_eja
1430
1431         TESTS:
1432
1433         Ensure that we can find an idempotent in a non-trivial algebra
1434         where there are non-nilpotent elements, or that we get the dumb
1435         solution in the trivial algebra::
1436
1437             sage: set_random_seed()
1438             sage: J = random_eja()
1439             sage: x = J.random_element()
1440             sage: while x.is_nilpotent() and not J.is_trivial():
1441             ....:     x = J.random_element()
1442             sage: c = x.subalgebra_idempotent()
1443             sage: c^2 == c
1444             True
1445
1446         """
1447         if self.parent().is_trivial():
1448             return self
1449
1450         if self.is_nilpotent():
1451             raise ValueError("this only works with non-nilpotent elements!")
1452
1453         J = self.subalgebra_generated_by()
1454         u = J(self)
1455
1456         # The image of the matrix of left-u^m-multiplication
1457         # will be minimal for some natural number s...
1458         s = 0
1459         minimal_dim = J.dimension()
1460         for i in range(1, minimal_dim):
1461             this_dim = (u**i).operator().matrix().image().dimension()
1462             if this_dim < minimal_dim:
1463                 minimal_dim = this_dim
1464                 s = i
1465
1466         # Now minimal_matrix should correspond to the smallest
1467         # non-zero subspace in Baes's (or really, Koecher's)
1468         # proposition.
1469         #
1470         # However, we need to restrict the matrix to work on the
1471         # subspace... or do we? Can't we just solve, knowing that
1472         # A(c) = u^(s+1) should have a solution in the big space,
1473         # too?
1474         #
1475         # Beware, solve_right() means that we're using COLUMN vectors.
1476         # Our FiniteDimensionalAlgebraElement superclass uses rows.
1477         u_next = u**(s+1)
1478         A = u_next.operator().matrix()
1479         c = J.from_vector(A.solve_right(u_next.to_vector()))
1480
1481         # Now c is the idempotent we want, but it still lives in the subalgebra.
1482         return c.superalgebra_element()
1483
1484
1485     def trace(self):
1486         """
1487         Return my trace, the sum of my eigenvalues.
1488
1489         In a trivial algebra, however you want to look at it, the trace is
1490         an empty sum for which we declare the result to be zero.
1491
1492         SETUP::
1493
1494             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1495             ....:                                  HadamardEJA,
1496             ....:                                  TrivialEJA,
1497             ....:                                  random_eja)
1498
1499         EXAMPLES::
1500
1501             sage: J = TrivialEJA()
1502             sage: J.zero().trace()
1503             0
1504
1505         ::
1506             sage: J = JordanSpinEJA(3)
1507             sage: x = sum(J.gens())
1508             sage: x.trace()
1509             2
1510
1511         ::
1512
1513             sage: J = HadamardEJA(5)
1514             sage: J.one().trace()
1515             5
1516
1517         TESTS:
1518
1519         The trace of an element is a real number::
1520
1521             sage: set_random_seed()
1522             sage: J = random_eja()
1523             sage: J.random_element().trace() in RLF
1524             True
1525
1526         The trace is linear::
1527
1528             sage: set_random_seed()
1529             sage: J = random_eja()
1530             sage: x,y = J.random_elements(2)
1531             sage: alpha = J.base_ring().random_element()
1532             sage: (alpha*x + y).trace() == alpha*x.trace() + y.trace()
1533             True
1534
1535         """
1536         P = self.parent()
1537         r = P.rank()
1538
1539         if r == 0:
1540             # Special case for the trivial algebra where
1541             # the trace is an empty sum.
1542             return P.base_ring().zero()
1543
1544         p = P._charpoly_coefficients()[r-1]
1545         # The _charpoly_coeff function already adds the factor of
1546         # -1 to ensure that _charpoly_coeff(r-1) is really what
1547         # appears in front of t^{r-1} in the charpoly. However,
1548         # we want the negative of THAT for the trace.
1549         return -p(*self.to_vector())
1550
1551
1552     def trace_inner_product(self, other):
1553         """
1554         Return the trace inner product of myself and ``other``.
1555
1556         SETUP::
1557
1558             sage: from mjo.eja.eja_algebra import random_eja
1559
1560         TESTS:
1561
1562         The trace inner product is commutative, bilinear, and associative::
1563
1564             sage: set_random_seed()
1565             sage: J = random_eja()
1566             sage: x,y,z = J.random_elements(3)
1567             sage: # commutative
1568             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1569             True
1570             sage: # bilinear
1571             sage: a = J.base_ring().random_element();
1572             sage: actual = (a*(x+z)).trace_inner_product(y)
1573             sage: expected = ( a*x.trace_inner_product(y) +
1574             ....:              a*z.trace_inner_product(y) )
1575             sage: actual == expected
1576             True
1577             sage: actual = x.trace_inner_product(a*(y+z))
1578             sage: expected = ( a*x.trace_inner_product(y) +
1579             ....:              a*x.trace_inner_product(z) )
1580             sage: actual == expected
1581             True
1582             sage: # associative
1583             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1584             True
1585
1586         """
1587         if not other in self.parent():
1588             raise TypeError("'other' must live in the same algebra")
1589
1590         return (self*other).trace()
1591
1592
1593     def trace_norm(self):
1594         """
1595         The norm of this element with respect to :meth:`trace_inner_product`.
1596
1597         SETUP::
1598
1599             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1600             ....:                                  HadamardEJA)
1601
1602         EXAMPLES::
1603
1604             sage: J = HadamardEJA(2)
1605             sage: x = sum(J.gens())
1606             sage: x.trace_norm()
1607             1.414213562373095?
1608
1609         ::
1610
1611             sage: J = JordanSpinEJA(4)
1612             sage: x = sum(J.gens())
1613             sage: x.trace_norm()
1614             2.828427124746190?
1615
1616         """
1617         return self.trace_inner_product(self).sqrt()
1618
1619
1620
1621 class CartesianProductEJAElement(FiniteDimensionalEJAElement):
1622
1623     def to_matrix(self):
1624         r"""
1625         SETUP::
1626
1627             sage: from mjo.eja.eja_algebra import (HadamardEJA,
1628             ....:                                  RealSymmetricEJA)
1629
1630         EXAMPLES::
1631
1632             sage: J1 = HadamardEJA(1)
1633             sage: J2 = RealSymmetricEJA(2)
1634             sage: J = cartesian_product([J1,J2])
1635             sage: x = sum(J.gens())
1636             sage: x.to_matrix()[0]
1637             [1]
1638             sage: x.to_matrix()[1]
1639             [                  1 0.7071067811865475?]
1640             [0.7071067811865475?                   1]
1641
1642         """
1643         B = self.parent().matrix_basis()
1644         W = self.parent().matrix_space()
1645
1646         # Aaaaand linear combinations don't work in Cartesian
1647         # product spaces, even though they provide a method
1648         # with that name.
1649         pairs = zip(B, self.to_vector())
1650         return sum( ( W(tuple(alpha*b_i for b_i in b))
1651                       for (b,alpha) in pairs ),
1652                     W.zero())