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