mjo/eja/eja_element.py

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