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