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