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 FiniteDimensionalEuclideanJordanAlgebraOperator
   6 from mjo.eja.eja_utils import _mat2vec
   7
   8 class FiniteDimensionalEuclideanJordanAlgebraElement(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 FiniteDimensionalEuclideanJordanAlgebraOperator(
1126                  P,
1127                  P,
1128                  L.matrix() )
1129
1130
1131     def quadratic_representation(self, other=None):
1132         """
1133         Return the quadratic representation of this element.
1134
1135         SETUP::
1136
1137             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1138             ....:                                  random_eja)
1139
1140         EXAMPLES:
1141
1142         The explicit form in the spin factor algebra is given by
1143         Alizadeh's Example 11.12::
1144
1145             sage: set_random_seed()
1146             sage: x = JordanSpinEJA.random_instance().random_element()
1147             sage: x_vec = x.to_vector()
1148             sage: Q = matrix.identity(x.base_ring(), 0)
1149             sage: n = x_vec.degree()
1150             sage: if n > 0:
1151             ....:     x0 = x_vec[0]
1152             ....:     x_bar = x_vec[1:]
1153             ....:     A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1154             ....:     B = 2*x0*x_bar.row()
1155             ....:     C = 2*x0*x_bar.column()
1156             ....:     D = matrix.identity(x.base_ring(), n-1)
1157             ....:     D = (x0^2 - x_bar.inner_product(x_bar))*D
1158             ....:     D = D + 2*x_bar.tensor_product(x_bar)
1159             ....:     Q = matrix.block(2,2,[A,B,C,D])
1160             sage: Q == x.quadratic_representation().matrix()
1161             True
1162
1163         Test all of the properties from Theorem 11.2 in Alizadeh::
1164
1165             sage: set_random_seed()
1166             sage: J = random_eja()
1167             sage: x,y = J.random_elements(2)
1168             sage: Lx = x.operator()
1169             sage: Lxx = (x*x).operator()
1170             sage: Qx = x.quadratic_representation()
1171             sage: Qy = y.quadratic_representation()
1172             sage: Qxy = x.quadratic_representation(y)
1173             sage: Qex = J.one().quadratic_representation(x)
1174             sage: n = ZZ.random_element(10)
1175             sage: Qxn = (x^n).quadratic_representation()
1176
1177         Property 1:
1178
1179             sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1180             True
1181
1182         Property 2 (multiply on the right for :trac:`28272`):
1183
1184             sage: alpha = J.base_ring().random_element()
1185             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1186             True
1187
1188         Property 3:
1189
1190             sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1191             True
1192
1193             sage: not x.is_invertible() or (
1194             ....:   ~Qx
1195             ....:   ==
1196             ....:   x.inverse().quadratic_representation() )
1197             True
1198
1199             sage: Qxy(J.one()) == x*y
1200             True
1201
1202         Property 4:
1203
1204             sage: not x.is_invertible() or (
1205             ....:   x.quadratic_representation(x.inverse())*Qx
1206             ....:   == Qx*x.quadratic_representation(x.inverse()) )
1207             True
1208
1209             sage: not x.is_invertible() or (
1210             ....:   x.quadratic_representation(x.inverse())*Qx
1211             ....:   ==
1212             ....:   2*Lx*Qex - Qx )
1213             True
1214
1215             sage: 2*Lx*Qex - Qx == Lxx
1216             True
1217
1218         Property 5:
1219
1220             sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1221             True
1222
1223         Property 6:
1224
1225             sage: Qxn == (Qx)^n
1226             True
1227
1228         Property 7:
1229
1230             sage: not x.is_invertible() or (
1231             ....:   Qx*x.inverse().operator() == Lx )
1232             True
1233
1234         Property 8:
1235
1236             sage: not x.operator_commutes_with(y) or (
1237             ....:   Qx(y)^n == Qxn(y^n) )
1238             True
1239
1240         """
1241         if other is None:
1242             other=self
1243         elif not other in self.parent():
1244             raise TypeError("'other' must live in the same algebra")
1245
1246         L = self.operator()
1247         M = other.operator()
1248         return ( L*M + M*L - (self*other).operator() )
1249
1250
1251
1252     def spectral_decomposition(self):
1253         """
1254         Return the unique spectral decomposition of this element.
1255
1256         ALGORITHM:
1257
1258         Following Faraut and Korányi's Theorem III.1.1, we restrict this
1259         element's left-multiplication-by operator to the subalgebra it
1260         generates. We then compute the spectral decomposition of that
1261         operator, and the spectral projectors we get back must be the
1262         left-multiplication-by operators for the idempotents we
1263         seek. Thus applying them to the identity element gives us those
1264         idempotents.
1265
1266         Since the eigenvalues are required to be distinct, we take
1267         the spectral decomposition of the zero element to be zero
1268         times the identity element of the algebra (which is idempotent,
1269         obviously).
1270
1271         SETUP::
1272
1273             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1274
1275         EXAMPLES:
1276
1277         The spectral decomposition of the identity is ``1`` times itself,
1278         and the spectral decomposition of zero is ``0`` times the identity::
1279
1280             sage: J = RealSymmetricEJA(3)
1281             sage: J.one()
1282             e0 + e2 + e5
1283             sage: J.one().spectral_decomposition()
1284             [(1, e0 + e2 + e5)]
1285             sage: J.zero().spectral_decomposition()
1286             [(0, e0 + e2 + e5)]
1287
1288         TESTS::
1289
1290             sage: J = RealSymmetricEJA(4)
1291             sage: x = sum(J.gens())
1292             sage: sd = x.spectral_decomposition()
1293             sage: l0 = sd[0][0]
1294             sage: l1 = sd[1][0]
1295             sage: c0 = sd[0][1]
1296             sage: c1 = sd[1][1]
1297             sage: c0.inner_product(c1) == 0
1298             True
1299             sage: c0.is_idempotent()
1300             True
1301             sage: c1.is_idempotent()
1302             True
1303             sage: c0 + c1 == J.one()
1304             True
1305             sage: l0*c0 + l1*c1 == x
1306             True
1307
1308         The spectral decomposition should work in subalgebras, too::
1309
1310             sage: J = RealSymmetricEJA(4)
1311             sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens()
1312             sage: A = 2*e5 - 2*e8
1313             sage: (lambda1, c1) = A.spectral_decomposition()[1]
1314             sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1315             sage: (f0, f1, f2) = J1.gens()
1316             sage: f0.spectral_decomposition()
1317             [(0, f2), (1, f0)]
1318
1319         """
1320         A = self.subalgebra_generated_by(orthonormalize_basis=True)
1321         result = []
1322         for (evalue, proj) in A(self).operator().spectral_decomposition():
1323             result.append( (evalue, proj(A.one()).superalgebra_element()) )
1324         return result
1325
1326     def subalgebra_generated_by(self, orthonormalize_basis=False):
1327         """
1328         Return the associative subalgebra of the parent EJA generated
1329         by this element.
1330
1331         Since our parent algebra is unital, we want "subalgebra" to mean
1332         "unital subalgebra" as well; thus the subalgebra that an element
1333         generates will itself be a Euclidean Jordan algebra after
1334         restricting the algebra operations appropriately. This is the
1335         subalgebra that Faraut and Korányi work with in section II.2, for
1336         example.
1337
1338         SETUP::
1339
1340             sage: from mjo.eja.eja_algebra import random_eja
1341
1342         TESTS:
1343
1344         This subalgebra, being composed of only powers, is associative::
1345
1346             sage: set_random_seed()
1347             sage: x0 = random_eja().random_element()
1348             sage: A = x0.subalgebra_generated_by()
1349             sage: x,y,z = A.random_elements(3)
1350             sage: (x*y)*z == x*(y*z)
1351             True
1352
1353         Squaring in the subalgebra should work the same as in
1354         the superalgebra::
1355
1356             sage: set_random_seed()
1357             sage: x = random_eja().random_element()
1358             sage: A = x.subalgebra_generated_by()
1359             sage: A(x^2) == A(x)*A(x)
1360             True
1361
1362         By definition, the subalgebra generated by the zero element is
1363         the one-dimensional algebra generated by the identity
1364         element... unless the original algebra was trivial, in which
1365         case the subalgebra is trivial too::
1366
1367             sage: set_random_seed()
1368             sage: A = random_eja().zero().subalgebra_generated_by()
1369             sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1370             True
1371
1372         """
1373         from mjo.eja.eja_element_subalgebra import FiniteDimensionalEuclideanJordanElementSubalgebra
1374         return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
1375
1376
1377     def subalgebra_idempotent(self):
1378         """
1379         Find an idempotent in the associative subalgebra I generate
1380         using Proposition 2.3.5 in Baes.
1381
1382         SETUP::
1383
1384             sage: from mjo.eja.eja_algebra import random_eja
1385
1386         TESTS:
1387
1388         Ensure that we can find an idempotent in a non-trivial algebra
1389         where there are non-nilpotent elements, or that we get the dumb
1390         solution in the trivial algebra::
1391
1392             sage: set_random_seed()
1393             sage: J = random_eja()
1394             sage: x = J.random_element()
1395             sage: while x.is_nilpotent() and not J.is_trivial():
1396             ....:     x = J.random_element()
1397             sage: c = x.subalgebra_idempotent()
1398             sage: c^2 == c
1399             True
1400
1401         """
1402         if self.parent().is_trivial():
1403             return self
1404
1405         if self.is_nilpotent():
1406             raise ValueError("this only works with non-nilpotent elements!")
1407
1408         J = self.subalgebra_generated_by()
1409         u = J(self)
1410
1411         # The image of the matrix of left-u^m-multiplication
1412         # will be minimal for some natural number s...
1413         s = 0
1414         minimal_dim = J.dimension()
1415         for i in range(1, minimal_dim):
1416             this_dim = (u**i).operator().matrix().image().dimension()
1417             if this_dim < minimal_dim:
1418                 minimal_dim = this_dim
1419                 s = i
1420
1421         # Now minimal_matrix should correspond to the smallest
1422         # non-zero subspace in Baes's (or really, Koecher's)
1423         # proposition.
1424         #
1425         # However, we need to restrict the matrix to work on the
1426         # subspace... or do we? Can't we just solve, knowing that
1427         # A(c) = u^(s+1) should have a solution in the big space,
1428         # too?
1429         #
1430         # Beware, solve_right() means that we're using COLUMN vectors.
1431         # Our FiniteDimensionalAlgebraElement superclass uses rows.
1432         u_next = u**(s+1)
1433         A = u_next.operator().matrix()
1434         c = J.from_vector(A.solve_right(u_next.to_vector()))
1435
1436         # Now c is the idempotent we want, but it still lives in the subalgebra.
1437         return c.superalgebra_element()
1438
1439
1440     def trace(self):
1441         """
1442         Return my trace, the sum of my eigenvalues.
1443
1444         In a trivial algebra, however you want to look at it, the trace is
1445         an empty sum for which we declare the result to be zero.
1446
1447         SETUP::
1448
1449             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1450             ....:                                  HadamardEJA,
1451             ....:                                  TrivialEJA,
1452             ....:                                  random_eja)
1453
1454         EXAMPLES::
1455
1456             sage: J = TrivialEJA()
1457             sage: J.zero().trace()
1458             0
1459
1460         ::
1461             sage: J = JordanSpinEJA(3)
1462             sage: x = sum(J.gens())
1463             sage: x.trace()
1464             2
1465
1466         ::
1467
1468             sage: J = HadamardEJA(5)
1469             sage: J.one().trace()
1470             5
1471
1472         TESTS:
1473
1474         The trace of an element is a real number::
1475
1476             sage: set_random_seed()
1477             sage: J = random_eja()
1478             sage: J.random_element().trace() in RLF
1479             True
1480
1481         """
1482         P = self.parent()
1483         r = P.rank()
1484
1485         if r == 0:
1486             # Special case for the trivial algebra where
1487             # the trace is an empty sum.
1488             return P.base_ring().zero()
1489
1490         p = P._charpoly_coefficients()[r-1]
1491         # The _charpoly_coeff function already adds the factor of
1492         # -1 to ensure that _charpoly_coeff(r-1) is really what
1493         # appears in front of t^{r-1} in the charpoly. However,
1494         # we want the negative of THAT for the trace.
1495         return -p(*self.to_vector())
1496
1497
1498     def trace_inner_product(self, other):
1499         """
1500         Return the trace inner product of myself and ``other``.
1501
1502         SETUP::
1503
1504             sage: from mjo.eja.eja_algebra import random_eja
1505
1506         TESTS:
1507
1508         The trace inner product is commutative, bilinear, and associative::
1509
1510             sage: set_random_seed()
1511             sage: J = random_eja()
1512             sage: x,y,z = J.random_elements(3)
1513             sage: # commutative
1514             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1515             True
1516             sage: # bilinear
1517             sage: a = J.base_ring().random_element();
1518             sage: actual = (a*(x+z)).trace_inner_product(y)
1519             sage: expected = ( a*x.trace_inner_product(y) +
1520             ....:              a*z.trace_inner_product(y) )
1521             sage: actual == expected
1522             True
1523             sage: actual = x.trace_inner_product(a*(y+z))
1524             sage: expected = ( a*x.trace_inner_product(y) +
1525             ....:              a*x.trace_inner_product(z) )
1526             sage: actual == expected
1527             True
1528             sage: # associative
1529             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1530             True
1531
1532         """
1533         if not other in self.parent():
1534             raise TypeError("'other' must live in the same algebra")
1535
1536         return (self*other).trace()
1537
1538
1539     def trace_norm(self):
1540         """
1541         The norm of this element with respect to :meth:`trace_inner_product`.
1542
1543         SETUP::
1544
1545             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1546             ....:                                  HadamardEJA)
1547
1548         EXAMPLES::
1549
1550             sage: J = HadamardEJA(2)
1551             sage: x = sum(J.gens())
1552             sage: x.trace_norm()
1553             1.414213562373095?
1554
1555         ::
1556
1557             sage: J = JordanSpinEJA(4)
1558             sage: x = sum(J.gens())
1559             sage: x.trace_norm()
1560             2.828427124746190?
1561
1562         """
1563         return self.trace_inner_product(self).sqrt()