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         We appeal to the quadratic representation as in Koecher's
 448         Theorem 12 in Chapter III, Section 5.
 449
 450         SETUP::
 451
 452             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 453             ....:                                  JordanSpinEJA,
 454             ....:                                  random_eja)
 455
 456         EXAMPLES:
 457
 458         The inverse in the spin factor algebra is given in Alizadeh's
 459         Example 11.11::
 460
 461             sage: set_random_seed()
 462             sage: J = JordanSpinEJA.random_instance()
 463             sage: x = J.random_element()
 464             sage: while not x.is_invertible():
 465             ....:     x = J.random_element()
 466             sage: x_vec = x.to_vector()
 467             sage: x0 = x_vec[:1]
 468             sage: x_bar = x_vec[1:]
 469             sage: coeff = x0.inner_product(x0) - x_bar.inner_product(x_bar)
 470             sage: x_inverse = x_vec.parent()(x0.list() + (-x_bar).list())
 471             sage: if not coeff.is_zero(): x_inverse = x_inverse/coeff
 472             sage: x.inverse() == J.from_vector(x_inverse)
 473             True
 474
 475         Trying to invert a non-invertible element throws an error:
 476
 477             sage: JordanSpinEJA(3).zero().inverse()
 478             Traceback (most recent call last):
 479             ...
 480             ValueError: element is not invertible
 481
 482         TESTS:
 483
 484         The identity element is its own inverse::
 485
 486             sage: set_random_seed()
 487             sage: J = random_eja()
 488             sage: J.one().inverse() == J.one()
 489             True
 490
 491         If an element has an inverse, it acts like one::
 492
 493             sage: set_random_seed()
 494             sage: J = random_eja()
 495             sage: x = J.random_element()
 496             sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
 497             True
 498
 499         The inverse of the inverse is what we started with::
 500
 501             sage: set_random_seed()
 502             sage: J = random_eja()
 503             sage: x = J.random_element()
 504             sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
 505             True
 506
 507         Proposition II.2.3 in Faraut and Korányi says that the inverse
 508         of an element is the inverse of its left-multiplication operator
 509         applied to the algebra's identity, when that inverse exists::
 510
 511             sage: set_random_seed()
 512             sage: J = random_eja()
 513             sage: x = J.random_element()
 514             sage: (not x.operator().is_invertible()) or (
 515             ....:    x.operator().inverse()(J.one()) == x.inverse() )
 516             True
 517
 518         Proposition II.2.4 in Faraut and Korányi gives a formula for
 519         the inverse based on the characteristic polynomial and the
 520         Cayley-Hamilton theorem for Euclidean Jordan algebras::
 521
 522             sage: set_random_seed()
 523             sage: J = ComplexHermitianEJA(3)
 524             sage: x = J.random_element()
 525             sage: while not x.is_invertible():
 526             ....:     x = J.random_element()
 527             sage: r = J.rank()
 528             sage: a = x.characteristic_polynomial().coefficients(sparse=False)
 529             sage: expected  = (-1)^(r+1)/x.det()
 530             sage: expected *= sum( a[i+1]*x^i for i in range(r) )
 531             sage: x.inverse() == expected
 532             True
 533
 534         """
 535         if not self.is_invertible():
 536             raise ValueError("element is not invertible")
 537
 538         if self.parent()._charpoly_coefficients.is_in_cache():
 539             # We can invert using our charpoly if it will be fast to
 540             # compute. If the coefficients are cached, our rank had
 541             # better be too!
 542             r = self.parent().rank()
 543             a = self.characteristic_polynomial().coefficients(sparse=False)
 544             return (-1)**(r+1)*sum(a[i+1]*self**i for i in range(r))/self.det()
 545
 546         return (~self.quadratic_representation())(self)
 547
 548
 549     def is_invertible(self):
 550         """
 551         Return whether or not this element is invertible.
 552
 553         ALGORITHM:
 554
 555         The usual way to do this is to check if the determinant is
 556         zero, but we need the characteristic polynomial for the
 557         determinant. The minimal polynomial is a lot easier to get,
 558         so we use Corollary 2 in Chapter V of Koecher to check
 559         whether or not the parent algebra's zero element is a root
 560         of this element's minimal polynomial.
 561
 562         That is... unless the coefficients of our algebra's
 563         "characteristic polynomial of" function are already cached!
 564         In that case, we just use the determinant (which will be fast
 565         as a result).
 566
 567         Beware that we can't use the superclass method, because it
 568         relies on the algebra being associative.
 569
 570         SETUP::
 571
 572             sage: from mjo.eja.eja_algebra import random_eja
 573
 574         TESTS:
 575
 576         The identity element is always invertible::
 577
 578             sage: set_random_seed()
 579             sage: J = random_eja()
 580             sage: J.one().is_invertible()
 581             True
 582
 583         The zero element is never invertible in a non-trivial algebra::
 584
 585             sage: set_random_seed()
 586             sage: J = random_eja()
 587             sage: (not J.is_trivial()) and J.zero().is_invertible()
 588             False
 589
 590         """
 591         if self.is_zero():
 592             if self.parent().is_trivial():
 593                 return True
 594             else:
 595                 return False
 596
 597         if self.parent()._charpoly_coefficients.is_in_cache():
 598             # The determinant will be quicker than computing the minimal
 599             # polynomial from scratch, most likely.
 600             return (not self.det().is_zero())
 601
 602         # In fact, we only need to know if the constant term is non-zero,
 603         # so we can pass in the field's zero element instead.
 604         zero = self.base_ring().zero()
 605         p = self.minimal_polynomial()
 606         return not (p(zero) == zero)
 607
 608
 609     def is_primitive_idempotent(self):
 610         """
 611         Return whether or not this element is a primitive (or minimal)
 612         idempotent.
 613
 614         A primitive idempotent is a non-zero idempotent that is not
 615         the sum of two other non-zero idempotents. Remark 2.7.15 in
 616         Baes shows that this is what he refers to as a "minimal
 617         idempotent."
 618
 619         An element of a Euclidean Jordan algebra is a minimal idempotent
 620         if it :meth:`is_idempotent` and if its Peirce subalgebra
 621         corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
 622         Proposition 2.7.17).
 623
 624         SETUP::
 625
 626             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 627             ....:                                  RealSymmetricEJA,
 628             ....:                                  TrivialEJA,
 629             ....:                                  random_eja)
 630
 631         WARNING::
 632
 633         This method is sloooooow.
 634
 635         EXAMPLES:
 636
 637         The spectral decomposition of a non-regular element should always
 638         contain at least one non-minimal idempotent::
 639
 640             sage: J = RealSymmetricEJA(3)
 641             sage: x = sum(J.gens())
 642             sage: x.is_regular()
 643             False
 644             sage: [ c.is_primitive_idempotent()
 645             ....:   for (l,c) in x.spectral_decomposition() ]
 646             [False, True]
 647
 648         On the other hand, the spectral decomposition of a regular
 649         element should always be in terms of minimal idempotents::
 650
 651             sage: J = JordanSpinEJA(4)
 652             sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
 653             sage: x.is_regular()
 654             True
 655             sage: [ c.is_primitive_idempotent()
 656             ....:   for (l,c) in x.spectral_decomposition() ]
 657             [True, True]
 658
 659         TESTS:
 660
 661         The identity element is minimal only in an EJA of rank one::
 662
 663             sage: set_random_seed()
 664             sage: J = random_eja()
 665             sage: J.rank() == 1 or not J.one().is_primitive_idempotent()
 666             True
 667
 668         A non-idempotent cannot be a minimal idempotent::
 669
 670             sage: set_random_seed()
 671             sage: J = JordanSpinEJA(4)
 672             sage: x = J.random_element()
 673             sage: (not x.is_idempotent()) and x.is_primitive_idempotent()
 674             False
 675
 676         Proposition 2.7.19 in Baes says that an element is a minimal
 677         idempotent if and only if it's idempotent with trace equal to
 678         unity::
 679
 680             sage: set_random_seed()
 681             sage: J = JordanSpinEJA(4)
 682             sage: x = J.random_element()
 683             sage: expected = (x.is_idempotent() and x.trace() == 1)
 684             sage: actual = x.is_primitive_idempotent()
 685             sage: actual == expected
 686             True
 687
 688         Primitive idempotents must be non-zero::
 689
 690             sage: set_random_seed()
 691             sage: J = random_eja()
 692             sage: J.zero().is_idempotent()
 693             True
 694             sage: J.zero().is_primitive_idempotent()
 695             False
 696
 697         As a consequence of the fact that primitive idempotents must
 698         be non-zero, there are no primitive idempotents in a trivial
 699         Euclidean Jordan algebra::
 700
 701             sage: J = TrivialEJA()
 702             sage: J.one().is_idempotent()
 703             True
 704             sage: J.one().is_primitive_idempotent()
 705             False
 706
 707         """
 708         if not self.is_idempotent():
 709             return False
 710
 711         if self.is_zero():
 712             return False
 713
 714         (_,_,J1) = self.parent().peirce_decomposition(self)
 715         return (J1.dimension() == 1)
 716
 717
 718     def is_nilpotent(self):
 719         """
 720         Return whether or not some power of this element is zero.
 721
 722         ALGORITHM:
 723
 724         We use Theorem 5 in Chapter III of Koecher, which says that
 725         an element ``x`` is nilpotent if and only if ``x.operator()``
 726         is nilpotent. And it is a basic fact of linear algebra that
 727         an operator on an `n`-dimensional space is nilpotent if and
 728         only if, when raised to the `n`th power, it equals the zero
 729         operator (for example, see Axler Corollary 8.8).
 730
 731         SETUP::
 732
 733             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 734             ....:                                  random_eja)
 735
 736         EXAMPLES::
 737
 738             sage: J = JordanSpinEJA(3)
 739             sage: x = sum(J.gens())
 740             sage: x.is_nilpotent()
 741             False
 742
 743         TESTS:
 744
 745         The identity element is never nilpotent, except in a trivial EJA::
 746
 747             sage: set_random_seed()
 748             sage: J = random_eja()
 749             sage: J.one().is_nilpotent() and not J.is_trivial()
 750             False
 751
 752         The additive identity is always nilpotent::
 753
 754             sage: set_random_seed()
 755             sage: random_eja().zero().is_nilpotent()
 756             True
 757
 758         """
 759         P = self.parent()
 760         zero_operator = P.zero().operator()
 761         return self.operator()**P.dimension() == zero_operator
 762
 763
 764     def is_regular(self):
 765         """
 766         Return whether or not this is a regular element.
 767
 768         SETUP::
 769
 770             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 771             ....:                                  random_eja)
 772
 773         EXAMPLES:
 774
 775         The identity element always has degree one, but any element
 776         linearly-independent from it is regular::
 777
 778             sage: J = JordanSpinEJA(5)
 779             sage: J.one().is_regular()
 780             False
 781             sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 782             sage: for x in J.gens():
 783             ....:     (J.one() + x).is_regular()
 784             False
 785             True
 786             True
 787             True
 788             True
 789
 790         TESTS:
 791
 792         The zero element should never be regular, unless the parent
 793         algebra has dimension less than or equal to one::
 794
 795             sage: set_random_seed()
 796             sage: J = random_eja()
 797             sage: J.dimension() <= 1 or not J.zero().is_regular()
 798             True
 799
 800         The unit element isn't regular unless the algebra happens to
 801         consist of only its scalar multiples::
 802
 803             sage: set_random_seed()
 804             sage: J = random_eja()
 805             sage: J.dimension() <= 1 or not J.one().is_regular()
 806             True
 807
 808         """
 809         return self.degree() == self.parent().rank()
 810
 811
 812     def degree(self):
 813         """
 814         Return the degree of this element, which is defined to be
 815         the degree of its minimal polynomial.
 816
 817         ALGORITHM:
 818
 819         For now, we skip the messy minimal polynomial computation
 820         and instead return the dimension of the vector space spanned
 821         by the powers of this element. The latter is a bit more
 822         straightforward to compute.
 823
 824         SETUP::
 825
 826             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 827             ....:                                  random_eja)
 828
 829         EXAMPLES::
 830
 831             sage: J = JordanSpinEJA(4)
 832             sage: J.one().degree()
 833             1
 834             sage: e0,e1,e2,e3 = J.gens()
 835             sage: (e0 - e1).degree()
 836             2
 837
 838         In the spin factor algebra (of rank two), all elements that
 839         aren't multiples of the identity are regular::
 840
 841             sage: set_random_seed()
 842             sage: J = JordanSpinEJA.random_instance()
 843             sage: n = J.dimension()
 844             sage: x = J.random_element()
 845             sage: x.degree() == min(n,2) or (x == x.coefficient(0)*J.one())
 846             True
 847
 848         TESTS:
 849
 850         The zero and unit elements are both of degree one in nontrivial
 851         algebras::
 852
 853             sage: set_random_seed()
 854             sage: J = random_eja()
 855             sage: d = J.zero().degree()
 856             sage: (J.is_trivial() and d == 0) or d == 1
 857             True
 858             sage: d = J.one().degree()
 859             sage: (J.is_trivial() and d == 0) or d == 1
 860             True
 861
 862         Our implementation agrees with the definition::
 863
 864             sage: set_random_seed()
 865             sage: x = random_eja().random_element()
 866             sage: x.degree() == x.minimal_polynomial().degree()
 867             True
 868
 869         """
 870         if self.is_zero() and not self.parent().is_trivial():
 871             # The minimal polynomial of zero in a nontrivial algebra
 872             # is "t"; in a trivial algebra it's "1" by convention
 873             # (it's an empty product).
 874             return 1
 875         return self.subalgebra_generated_by().dimension()
 876
 877
 878     def left_matrix(self):
 879         """
 880         Our parent class defines ``left_matrix`` and ``matrix``
 881         methods whose names are misleading. We don't want them.
 882         """
 883         raise NotImplementedError("use operator().matrix() instead")
 884
 885     matrix = left_matrix
 886
 887
 888     def minimal_polynomial(self):
 889         """
 890         Return the minimal polynomial of this element,
 891         as a function of the variable `t`.
 892
 893         ALGORITHM:
 894
 895         We restrict ourselves to the associative subalgebra
 896         generated by this element, and then return the minimal
 897         polynomial of this element's operator matrix (in that
 898         subalgebra). This works by Baes Proposition 2.3.16.
 899
 900         SETUP::
 901
 902             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 903             ....:                                  RealSymmetricEJA,
 904             ....:                                  TrivialEJA,
 905             ....:                                  random_eja)
 906
 907         EXAMPLES:
 908
 909         Keeping in mind that the polynomial ``1`` evaluates the identity
 910         element (also the zero element) of the trivial algebra, it is clear
 911         that the polynomial ``1`` is the minimal polynomial of the only
 912         element in a trivial algebra::
 913
 914             sage: J = TrivialEJA()
 915             sage: J.one().minimal_polynomial()
 916             1
 917             sage: J.zero().minimal_polynomial()
 918             1
 919
 920         TESTS:
 921
 922         The minimal polynomial of the identity and zero elements are
 923         always the same, except in trivial algebras where the minimal
 924         polynomial of the unit/zero element is ``1``::
 925
 926             sage: set_random_seed()
 927             sage: J = random_eja()
 928             sage: mu = J.one().minimal_polynomial()
 929             sage: t = mu.parent().gen()
 930             sage: mu + int(J.is_trivial())*(t-2)
 931             t - 1
 932             sage: mu = J.zero().minimal_polynomial()
 933             sage: t = mu.parent().gen()
 934             sage: mu + int(J.is_trivial())*(t-1)
 935             t
 936
 937         The degree of an element is (by one definition) the degree
 938         of its minimal polynomial::
 939
 940             sage: set_random_seed()
 941             sage: x = random_eja().random_element()
 942             sage: x.degree() == x.minimal_polynomial().degree()
 943             True
 944
 945         The minimal polynomial and the characteristic polynomial coincide
 946         and are known (see Alizadeh, Example 11.11) for all elements of
 947         the spin factor algebra that aren't scalar multiples of the
 948         identity. We require the dimension of the algebra to be at least
 949         two here so that said elements actually exist::
 950
 951             sage: set_random_seed()
 952             sage: n_max = max(2, JordanSpinEJA._max_random_instance_size())
 953             sage: n = ZZ.random_element(2, n_max)
 954             sage: J = JordanSpinEJA(n)
 955             sage: y = J.random_element()
 956             sage: while y == y.coefficient(0)*J.one():
 957             ....:     y = J.random_element()
 958             sage: y0 = y.to_vector()[0]
 959             sage: y_bar = y.to_vector()[1:]
 960             sage: actual = y.minimal_polynomial()
 961             sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
 962             sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
 963             sage: bool(actual == expected)
 964             True
 965
 966         The minimal polynomial should always kill its element::
 967
 968             sage: set_random_seed()
 969             sage: x = random_eja().random_element()
 970             sage: p = x.minimal_polynomial()
 971             sage: x.apply_univariate_polynomial(p)
 972             0
 973
 974         The minimal polynomial is invariant under a change of basis,
 975         and in particular, a re-scaling of the basis::
 976
 977             sage: set_random_seed()
 978             sage: n_max = RealSymmetricEJA._max_random_instance_size()
 979             sage: n = ZZ.random_element(1, n_max)
 980             sage: J1 = RealSymmetricEJA(n)
 981             sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
 982             sage: X = random_matrix(AA,n)
 983             sage: X = X*X.transpose()
 984             sage: x1 = J1(X)
 985             sage: x2 = J2(X)
 986             sage: x1.minimal_polynomial() == x2.minimal_polynomial()
 987             True
 988
 989         """
 990         if self.is_zero():
 991             # We would generate a zero-dimensional subalgebra
 992             # where the minimal polynomial would be constant.
 993             # That might be correct, but only if *this* algebra
 994             # is trivial too.
 995             if not self.parent().is_trivial():
 996                 # Pretty sure we know what the minimal polynomial of
 997                 # the zero operator is going to be. This ensures
 998                 # consistency of e.g. the polynomial variable returned
 999                 # in the "normal" case without us having to think about it.
1000                 return self.operator().minimal_polynomial()
1001
1002         A = self.subalgebra_generated_by(orthonormalize_basis=False)
1003         return A(self).operator().minimal_polynomial()
1004
1005
1006
1007     def to_matrix(self):
1008         """
1009         Return an (often more natural) representation of this element as a
1010         matrix.
1011
1012         Every finite-dimensional Euclidean Jordan Algebra is a direct
1013         sum of five simple algebras, four of which comprise Hermitian
1014         matrices. This method returns a "natural" matrix
1015         representation of this element as either a Hermitian matrix or
1016         column vector.
1017
1018         SETUP::
1019
1020             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1021             ....:                                  QuaternionHermitianEJA)
1022
1023         EXAMPLES::
1024
1025             sage: J = ComplexHermitianEJA(3)
1026             sage: J.one()
1027             e0 + e3 + e8
1028             sage: J.one().to_matrix()
1029             [1 0 0 0 0 0]
1030             [0 1 0 0 0 0]
1031             [0 0 1 0 0 0]
1032             [0 0 0 1 0 0]
1033             [0 0 0 0 1 0]
1034             [0 0 0 0 0 1]
1035
1036         ::
1037
1038             sage: J = QuaternionHermitianEJA(2)
1039             sage: J.one()
1040             e0 + e5
1041             sage: J.one().to_matrix()
1042             [1 0 0 0 0 0 0 0]
1043             [0 1 0 0 0 0 0 0]
1044             [0 0 1 0 0 0 0 0]
1045             [0 0 0 1 0 0 0 0]
1046             [0 0 0 0 1 0 0 0]
1047             [0 0 0 0 0 1 0 0]
1048             [0 0 0 0 0 0 1 0]
1049             [0 0 0 0 0 0 0 1]
1050
1051         """
1052         B = self.parent().matrix_basis()
1053         W = self.parent().matrix_space()
1054
1055         # This is just a manual "from_vector()", but of course
1056         # matrix spaces aren't vector spaces in sage, so they
1057         # don't have a from_vector() method.
1058         return W.linear_combination( zip(B, self.to_vector()) )
1059
1060
1061     def norm(self):
1062         """
1063         The norm of this element with respect to :meth:`inner_product`.
1064
1065         SETUP::
1066
1067             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1068             ....:                                  HadamardEJA)
1069
1070         EXAMPLES::
1071
1072             sage: J = HadamardEJA(2)
1073             sage: x = sum(J.gens())
1074             sage: x.norm()
1075             1.414213562373095?
1076
1077         ::
1078
1079             sage: J = JordanSpinEJA(4)
1080             sage: x = sum(J.gens())
1081             sage: x.norm()
1082             2
1083
1084         """
1085         return self.inner_product(self).sqrt()
1086
1087
1088     def operator(self):
1089         """
1090         Return the left-multiplication-by-this-element
1091         operator on the ambient algebra.
1092
1093         SETUP::
1094
1095             sage: from mjo.eja.eja_algebra import random_eja
1096
1097         TESTS::
1098
1099             sage: set_random_seed()
1100             sage: J = random_eja()
1101             sage: x,y = J.random_elements(2)
1102             sage: x.operator()(y) == x*y
1103             True
1104             sage: y.operator()(x) == x*y
1105             True
1106
1107         """
1108         P = self.parent()
1109         left_mult_by_self = lambda y: self*y
1110         L = P.module_morphism(function=left_mult_by_self, codomain=P)
1111         return FiniteDimensionalEuclideanJordanAlgebraOperator(
1112                  P,
1113                  P,
1114                  L.matrix() )
1115
1116
1117     def quadratic_representation(self, other=None):
1118         """
1119         Return the quadratic representation of this element.
1120
1121         SETUP::
1122
1123             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1124             ....:                                  random_eja)
1125
1126         EXAMPLES:
1127
1128         The explicit form in the spin factor algebra is given by
1129         Alizadeh's Example 11.12::
1130
1131             sage: set_random_seed()
1132             sage: x = JordanSpinEJA.random_instance().random_element()
1133             sage: x_vec = x.to_vector()
1134             sage: Q = matrix.identity(x.base_ring(), 0)
1135             sage: n = x_vec.degree()
1136             sage: if n > 0:
1137             ....:     x0 = x_vec[0]
1138             ....:     x_bar = x_vec[1:]
1139             ....:     A = matrix(x.base_ring(), 1, [x_vec.inner_product(x_vec)])
1140             ....:     B = 2*x0*x_bar.row()
1141             ....:     C = 2*x0*x_bar.column()
1142             ....:     D = matrix.identity(x.base_ring(), n-1)
1143             ....:     D = (x0^2 - x_bar.inner_product(x_bar))*D
1144             ....:     D = D + 2*x_bar.tensor_product(x_bar)
1145             ....:     Q = matrix.block(2,2,[A,B,C,D])
1146             sage: Q == x.quadratic_representation().matrix()
1147             True
1148
1149         Test all of the properties from Theorem 11.2 in Alizadeh::
1150
1151             sage: set_random_seed()
1152             sage: J = random_eja()
1153             sage: x,y = J.random_elements(2)
1154             sage: Lx = x.operator()
1155             sage: Lxx = (x*x).operator()
1156             sage: Qx = x.quadratic_representation()
1157             sage: Qy = y.quadratic_representation()
1158             sage: Qxy = x.quadratic_representation(y)
1159             sage: Qex = J.one().quadratic_representation(x)
1160             sage: n = ZZ.random_element(10)
1161             sage: Qxn = (x^n).quadratic_representation()
1162
1163         Property 1:
1164
1165             sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1166             True
1167
1168         Property 2 (multiply on the right for :trac:`28272`):
1169
1170             sage: alpha = J.base_ring().random_element()
1171             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1172             True
1173
1174         Property 3:
1175
1176             sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1177             True
1178
1179             sage: not x.is_invertible() or (
1180             ....:   ~Qx
1181             ....:   ==
1182             ....:   x.inverse().quadratic_representation() )
1183             True
1184
1185             sage: Qxy(J.one()) == x*y
1186             True
1187
1188         Property 4:
1189
1190             sage: not x.is_invertible() or (
1191             ....:   x.quadratic_representation(x.inverse())*Qx
1192             ....:   == Qx*x.quadratic_representation(x.inverse()) )
1193             True
1194
1195             sage: not x.is_invertible() or (
1196             ....:   x.quadratic_representation(x.inverse())*Qx
1197             ....:   ==
1198             ....:   2*Lx*Qex - Qx )
1199             True
1200
1201             sage: 2*Lx*Qex - Qx == Lxx
1202             True
1203
1204         Property 5:
1205
1206             sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1207             True
1208
1209         Property 6:
1210
1211             sage: Qxn == (Qx)^n
1212             True
1213
1214         Property 7:
1215
1216             sage: not x.is_invertible() or (
1217             ....:   Qx*x.inverse().operator() == Lx )
1218             True
1219
1220         Property 8:
1221
1222             sage: not x.operator_commutes_with(y) or (
1223             ....:   Qx(y)^n == Qxn(y^n) )
1224             True
1225
1226         """
1227         if other is None:
1228             other=self
1229         elif not other in self.parent():
1230             raise TypeError("'other' must live in the same algebra")
1231
1232         L = self.operator()
1233         M = other.operator()
1234         return ( L*M + M*L - (self*other).operator() )
1235
1236
1237
1238     def spectral_decomposition(self):
1239         """
1240         Return the unique spectral decomposition of this element.
1241
1242         ALGORITHM:
1243
1244         Following Faraut and Korányi's Theorem III.1.1, we restrict this
1245         element's left-multiplication-by operator to the subalgebra it
1246         generates. We then compute the spectral decomposition of that
1247         operator, and the spectral projectors we get back must be the
1248         left-multiplication-by operators for the idempotents we
1249         seek. Thus applying them to the identity element gives us those
1250         idempotents.
1251
1252         Since the eigenvalues are required to be distinct, we take
1253         the spectral decomposition of the zero element to be zero
1254         times the identity element of the algebra (which is idempotent,
1255         obviously).
1256
1257         SETUP::
1258
1259             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1260
1261         EXAMPLES:
1262
1263         The spectral decomposition of the identity is ``1`` times itself,
1264         and the spectral decomposition of zero is ``0`` times the identity::
1265
1266             sage: J = RealSymmetricEJA(3)
1267             sage: J.one()
1268             e0 + e2 + e5
1269             sage: J.one().spectral_decomposition()
1270             [(1, e0 + e2 + e5)]
1271             sage: J.zero().spectral_decomposition()
1272             [(0, e0 + e2 + e5)]
1273
1274         TESTS::
1275
1276             sage: J = RealSymmetricEJA(4)
1277             sage: x = sum(J.gens())
1278             sage: sd = x.spectral_decomposition()
1279             sage: l0 = sd[0][0]
1280             sage: l1 = sd[1][0]
1281             sage: c0 = sd[0][1]
1282             sage: c1 = sd[1][1]
1283             sage: c0.inner_product(c1) == 0
1284             True
1285             sage: c0.is_idempotent()
1286             True
1287             sage: c1.is_idempotent()
1288             True
1289             sage: c0 + c1 == J.one()
1290             True
1291             sage: l0*c0 + l1*c1 == x
1292             True
1293
1294         The spectral decomposition should work in subalgebras, too::
1295
1296             sage: J = RealSymmetricEJA(4)
1297             sage: (e0, e1, e2, e3, e4, e5, e6, e7, e8, e9) = J.gens()
1298             sage: A = 2*e5 - 2*e8
1299             sage: (lambda1, c1) = A.spectral_decomposition()[1]
1300             sage: (J0, J5, J1) = J.peirce_decomposition(c1)
1301             sage: (f0, f1, f2) = J1.gens()
1302             sage: f0.spectral_decomposition()
1303             [(0, f2), (1, f0)]
1304
1305         """
1306         A = self.subalgebra_generated_by(orthonormalize_basis=True)
1307         result = []
1308         for (evalue, proj) in A(self).operator().spectral_decomposition():
1309             result.append( (evalue, proj(A.one()).superalgebra_element()) )
1310         return result
1311
1312     def subalgebra_generated_by(self, orthonormalize_basis=False):
1313         """
1314         Return the associative subalgebra of the parent EJA generated
1315         by this element.
1316
1317         Since our parent algebra is unital, we want "subalgebra" to mean
1318         "unital subalgebra" as well; thus the subalgebra that an element
1319         generates will itself be a Euclidean Jordan algebra after
1320         restricting the algebra operations appropriately. This is the
1321         subalgebra that Faraut and Korányi work with in section II.2, for
1322         example.
1323
1324         SETUP::
1325
1326             sage: from mjo.eja.eja_algebra import random_eja
1327
1328         TESTS:
1329
1330         This subalgebra, being composed of only powers, is associative::
1331
1332             sage: set_random_seed()
1333             sage: x0 = random_eja().random_element()
1334             sage: A = x0.subalgebra_generated_by()
1335             sage: x,y,z = A.random_elements(3)
1336             sage: (x*y)*z == x*(y*z)
1337             True
1338
1339         Squaring in the subalgebra should work the same as in
1340         the superalgebra::
1341
1342             sage: set_random_seed()
1343             sage: x = random_eja().random_element()
1344             sage: A = x.subalgebra_generated_by()
1345             sage: A(x^2) == A(x)*A(x)
1346             True
1347
1348         By definition, the subalgebra generated by the zero element is
1349         the one-dimensional algebra generated by the identity
1350         element... unless the original algebra was trivial, in which
1351         case the subalgebra is trivial too::
1352
1353             sage: set_random_seed()
1354             sage: A = random_eja().zero().subalgebra_generated_by()
1355             sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1356             True
1357
1358         """
1359         from mjo.eja.eja_element_subalgebra import FiniteDimensionalEuclideanJordanElementSubalgebra
1360         return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
1361
1362
1363     def subalgebra_idempotent(self):
1364         """
1365         Find an idempotent in the associative subalgebra I generate
1366         using Proposition 2.3.5 in Baes.
1367
1368         SETUP::
1369
1370             sage: from mjo.eja.eja_algebra import random_eja
1371
1372         TESTS:
1373
1374         Ensure that we can find an idempotent in a non-trivial algebra
1375         where there are non-nilpotent elements, or that we get the dumb
1376         solution in the trivial algebra::
1377
1378             sage: set_random_seed()
1379             sage: J = random_eja()
1380             sage: x = J.random_element()
1381             sage: while x.is_nilpotent() and not J.is_trivial():
1382             ....:     x = J.random_element()
1383             sage: c = x.subalgebra_idempotent()
1384             sage: c^2 == c
1385             True
1386
1387         """
1388         if self.parent().is_trivial():
1389             return self
1390
1391         if self.is_nilpotent():
1392             raise ValueError("this only works with non-nilpotent elements!")
1393
1394         J = self.subalgebra_generated_by()
1395         u = J(self)
1396
1397         # The image of the matrix of left-u^m-multiplication
1398         # will be minimal for some natural number s...
1399         s = 0
1400         minimal_dim = J.dimension()
1401         for i in range(1, minimal_dim):
1402             this_dim = (u**i).operator().matrix().image().dimension()
1403             if this_dim < minimal_dim:
1404                 minimal_dim = this_dim
1405                 s = i
1406
1407         # Now minimal_matrix should correspond to the smallest
1408         # non-zero subspace in Baes's (or really, Koecher's)
1409         # proposition.
1410         #
1411         # However, we need to restrict the matrix to work on the
1412         # subspace... or do we? Can't we just solve, knowing that
1413         # A(c) = u^(s+1) should have a solution in the big space,
1414         # too?
1415         #
1416         # Beware, solve_right() means that we're using COLUMN vectors.
1417         # Our FiniteDimensionalAlgebraElement superclass uses rows.
1418         u_next = u**(s+1)
1419         A = u_next.operator().matrix()
1420         c = J.from_vector(A.solve_right(u_next.to_vector()))
1421
1422         # Now c is the idempotent we want, but it still lives in the subalgebra.
1423         return c.superalgebra_element()
1424
1425
1426     def trace(self):
1427         """
1428         Return my trace, the sum of my eigenvalues.
1429
1430         In a trivial algebra, however you want to look at it, the trace is
1431         an empty sum for which we declare the result to be zero.
1432
1433         SETUP::
1434
1435             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1436             ....:                                  HadamardEJA,
1437             ....:                                  TrivialEJA,
1438             ....:                                  random_eja)
1439
1440         EXAMPLES::
1441
1442             sage: J = TrivialEJA()
1443             sage: J.zero().trace()
1444             0
1445
1446         ::
1447             sage: J = JordanSpinEJA(3)
1448             sage: x = sum(J.gens())
1449             sage: x.trace()
1450             2
1451
1452         ::
1453
1454             sage: J = HadamardEJA(5)
1455             sage: J.one().trace()
1456             5
1457
1458         TESTS:
1459
1460         The trace of an element is a real number::
1461
1462             sage: set_random_seed()
1463             sage: J = random_eja()
1464             sage: J.random_element().trace() in RLF
1465             True
1466
1467         """
1468         P = self.parent()
1469         r = P.rank()
1470
1471         if r == 0:
1472             # Special case for the trivial algebra where
1473             # the trace is an empty sum.
1474             return P.base_ring().zero()
1475
1476         p = P._charpoly_coefficients()[r-1]
1477         # The _charpoly_coeff function already adds the factor of
1478         # -1 to ensure that _charpoly_coeff(r-1) is really what
1479         # appears in front of t^{r-1} in the charpoly. However,
1480         # we want the negative of THAT for the trace.
1481         return -p(*self.to_vector())
1482
1483
1484     def trace_inner_product(self, other):
1485         """
1486         Return the trace inner product of myself and ``other``.
1487
1488         SETUP::
1489
1490             sage: from mjo.eja.eja_algebra import random_eja
1491
1492         TESTS:
1493
1494         The trace inner product is commutative, bilinear, and associative::
1495
1496             sage: set_random_seed()
1497             sage: J = random_eja()
1498             sage: x,y,z = J.random_elements(3)
1499             sage: # commutative
1500             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1501             True
1502             sage: # bilinear
1503             sage: a = J.base_ring().random_element();
1504             sage: actual = (a*(x+z)).trace_inner_product(y)
1505             sage: expected = ( a*x.trace_inner_product(y) +
1506             ....:              a*z.trace_inner_product(y) )
1507             sage: actual == expected
1508             True
1509             sage: actual = x.trace_inner_product(a*(y+z))
1510             sage: expected = ( a*x.trace_inner_product(y) +
1511             ....:              a*x.trace_inner_product(z) )
1512             sage: actual == expected
1513             True
1514             sage: # associative
1515             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1516             True
1517
1518         """
1519         if not other in self.parent():
1520             raise TypeError("'other' must live in the same algebra")
1521
1522         return (self*other).trace()
1523
1524
1525     def trace_norm(self):
1526         """
1527         The norm of this element with respect to :meth:`trace_inner_product`.
1528
1529         SETUP::
1530
1531             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1532             ....:                                  HadamardEJA)
1533
1534         EXAMPLES::
1535
1536             sage: J = HadamardEJA(2)
1537             sage: x = sum(J.gens())
1538             sage: x.trace_norm()
1539             1.414213562373095?
1540
1541         ::
1542
1543             sage: J = JordanSpinEJA(4)
1544             sage: x = sum(J.gens())
1545             sage: x.trace_norm()
1546             2.828427124746190?
1547
1548         """
1549         return self.trace_inner_product(self).sqrt()