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_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 Koranyi, 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.operator()**(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 (RealCartesianProductEJA,
  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 = RealCartesianProductEJA(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(QQ, '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 RealCartesianProductEJA
 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 = RealCartesianProductEJA(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 = RealCartesianProductEJA(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 = RealCartesianProductEJA(3).random_element()
 164             sage: p = x.characteristic_polynomial()
 165             sage: x.apply_univariate_polynomial(p)
 166             0
 167
 168         """
 169         p = self.parent().characteristic_polynomial()
 170         return p(*self.to_vector())
 171
 172
 173     def inner_product(self, other):
 174         """
 175         Return the parent algebra's inner product of myself and ``other``.
 176
 177         SETUP::
 178
 179             sage: from mjo.eja.eja_algebra import (
 180             ....:   ComplexHermitianEJA,
 181             ....:   JordanSpinEJA,
 182             ....:   QuaternionHermitianEJA,
 183             ....:   RealSymmetricEJA,
 184             ....:   random_eja)
 185
 186         EXAMPLES:
 187
 188         The inner product in the Jordan spin algebra is the usual
 189         inner product on `R^n` (this example only works because the
 190         basis for the Jordan algebra is the standard basis in `R^n`)::
 191
 192             sage: J = JordanSpinEJA(3)
 193             sage: x = vector(QQ,[1,2,3])
 194             sage: y = vector(QQ,[4,5,6])
 195             sage: x.inner_product(y)
 196             32
 197             sage: J.from_vector(x).inner_product(J.from_vector(y))
 198             32
 199
 200         The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
 201         multiplication is the usual matrix multiplication in `S^n`,
 202         so the inner product of the identity matrix with itself
 203         should be the `n`::
 204
 205             sage: J = RealSymmetricEJA(3)
 206             sage: J.one().inner_product(J.one())
 207             3
 208
 209         Likewise, the inner product on `C^n` is `<X,Y> =
 210         Re(trace(X*Y))`, where we must necessarily take the real
 211         part because the product of Hermitian matrices may not be
 212         Hermitian::
 213
 214             sage: J = ComplexHermitianEJA(3)
 215             sage: J.one().inner_product(J.one())
 216             3
 217
 218         Ditto for the quaternions::
 219
 220             sage: J = QuaternionHermitianEJA(3)
 221             sage: J.one().inner_product(J.one())
 222             3
 223
 224         TESTS:
 225
 226         Ensure that we can always compute an inner product, and that
 227         it gives us back a real number::
 228
 229             sage: set_random_seed()
 230             sage: J = random_eja()
 231             sage: x = J.random_element()
 232             sage: y = J.random_element()
 233             sage: x.inner_product(y) in RR
 234             True
 235
 236         """
 237         P = self.parent()
 238         if not other in P:
 239             raise TypeError("'other' must live in the same algebra")
 240
 241         return P.inner_product(self, other)
 242
 243
 244     def operator_commutes_with(self, other):
 245         """
 246         Return whether or not this element operator-commutes
 247         with ``other``.
 248
 249         SETUP::
 250
 251             sage: from mjo.eja.eja_algebra import random_eja
 252
 253         EXAMPLES:
 254
 255         The definition of a Jordan algebra says that any element
 256         operator-commutes with its square::
 257
 258             sage: set_random_seed()
 259             sage: x = random_eja().random_element()
 260             sage: x.operator_commutes_with(x^2)
 261             True
 262
 263         TESTS:
 264
 265         Test Lemma 1 from Chapter III of Koecher::
 266
 267             sage: set_random_seed()
 268             sage: J = random_eja()
 269             sage: u = J.random_element()
 270             sage: v = J.random_element()
 271             sage: lhs = u.operator_commutes_with(u*v)
 272             sage: rhs = v.operator_commutes_with(u^2)
 273             sage: lhs == rhs
 274             True
 275
 276         Test the first polarization identity from my notes, Koecher
 277         Chapter III, or from Baes (2.3)::
 278
 279             sage: set_random_seed()
 280             sage: J = random_eja()
 281             sage: x = J.random_element()
 282             sage: y = J.random_element()
 283             sage: Lx = x.operator()
 284             sage: Ly = y.operator()
 285             sage: Lxx = (x*x).operator()
 286             sage: Lxy = (x*y).operator()
 287             sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
 288             True
 289
 290         Test the second polarization identity from my notes or from
 291         Baes (2.4)::
 292
 293             sage: set_random_seed()
 294             sage: J = random_eja()
 295             sage: x = J.random_element()
 296             sage: y = J.random_element()
 297             sage: z = J.random_element()
 298             sage: Lx = x.operator()
 299             sage: Ly = y.operator()
 300             sage: Lz = z.operator()
 301             sage: Lzy = (z*y).operator()
 302             sage: Lxy = (x*y).operator()
 303             sage: Lxz = (x*z).operator()
 304             sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
 305             True
 306
 307         Test the third polarization identity from my notes or from
 308         Baes (2.5)::
 309
 310             sage: set_random_seed()
 311             sage: J = random_eja()
 312             sage: u = J.random_element()
 313             sage: y = J.random_element()
 314             sage: z = J.random_element()
 315             sage: Lu = u.operator()
 316             sage: Ly = y.operator()
 317             sage: Lz = z.operator()
 318             sage: Lzy = (z*y).operator()
 319             sage: Luy = (u*y).operator()
 320             sage: Luz = (u*z).operator()
 321             sage: Luyz = (u*(y*z)).operator()
 322             sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
 323             sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
 324             sage: bool(lhs == rhs)
 325             True
 326
 327         """
 328         if not other in self.parent():
 329             raise TypeError("'other' must live in the same algebra")
 330
 331         A = self.operator()
 332         B = other.operator()
 333         return (A*B == B*A)
 334
 335
 336     def det(self):
 337         """
 338         Return my determinant, the product of my eigenvalues.
 339
 340         SETUP::
 341
 342             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 343             ....:                                  random_eja)
 344
 345         EXAMPLES::
 346
 347             sage: J = JordanSpinEJA(2)
 348             sage: e0,e1 = J.gens()
 349             sage: x = sum( J.gens() )
 350             sage: x.det()
 351             0
 352
 353         ::
 354
 355             sage: J = JordanSpinEJA(3)
 356             sage: e0,e1,e2 = J.gens()
 357             sage: x = sum( J.gens() )
 358             sage: x.det()
 359             -1
 360
 361         TESTS:
 362
 363         An element is invertible if and only if its determinant is
 364         non-zero::
 365
 366             sage: set_random_seed()
 367             sage: x = random_eja().random_element()
 368             sage: x.is_invertible() == (x.det() != 0)
 369             True
 370
 371         """
 372         P = self.parent()
 373         r = P.rank()
 374         p = P._charpoly_coeff(0)
 375         # The _charpoly_coeff function already adds the factor of
 376         # -1 to ensure that _charpoly_coeff(0) is really what
 377         # appears in front of t^{0} in the charpoly. However,
 378         # we want (-1)^r times THAT for the determinant.
 379         return ((-1)**r)*p(*self.to_vector())
 380
 381
 382     def inverse(self):
 383         """
 384         Return the Jordan-multiplicative inverse of this element.
 385
 386         ALGORITHM:
 387
 388         We appeal to the quadratic representation as in Koecher's
 389         Theorem 12 in Chapter III, Section 5.
 390
 391         SETUP::
 392
 393             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 394             ....:                                  random_eja)
 395
 396         EXAMPLES:
 397
 398         The inverse in the spin factor algebra is given in Alizadeh's
 399         Example 11.11::
 400
 401             sage: set_random_seed()
 402             sage: n = ZZ.random_element(1,10)
 403             sage: J = JordanSpinEJA(n)
 404             sage: x = J.random_element()
 405             sage: while not x.is_invertible():
 406             ....:     x = J.random_element()
 407             sage: x_vec = x.to_vector()
 408             sage: x0 = x_vec[0]
 409             sage: x_bar = x_vec[1:]
 410             sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
 411             sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
 412             sage: x_inverse = coeff*inv_vec
 413             sage: x.inverse() == J.from_vector(x_inverse)
 414             True
 415
 416         TESTS:
 417
 418         The identity element is its own inverse::
 419
 420             sage: set_random_seed()
 421             sage: J = random_eja()
 422             sage: J.one().inverse() == J.one()
 423             True
 424
 425         If an element has an inverse, it acts like one::
 426
 427             sage: set_random_seed()
 428             sage: J = random_eja()
 429             sage: x = J.random_element()
 430             sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
 431             True
 432
 433         The inverse of the inverse is what we started with::
 434
 435             sage: set_random_seed()
 436             sage: J = random_eja()
 437             sage: x = J.random_element()
 438             sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
 439             True
 440
 441         The zero element is never invertible::
 442
 443             sage: set_random_seed()
 444             sage: J = random_eja().zero().inverse()
 445             Traceback (most recent call last):
 446             ...
 447             ValueError: element is not invertible
 448
 449         """
 450         if not self.is_invertible():
 451             raise ValueError("element is not invertible")
 452
 453         return (~self.quadratic_representation())(self)
 454
 455
 456     def is_invertible(self):
 457         """
 458         Return whether or not this element is invertible.
 459
 460         ALGORITHM:
 461
 462         The usual way to do this is to check if the determinant is
 463         zero, but we need the characteristic polynomial for the
 464         determinant. The minimal polynomial is a lot easier to get,
 465         so we use Corollary 2 in Chapter V of Koecher to check
 466         whether or not the paren't algebra's zero element is a root
 467         of this element's minimal polynomial.
 468
 469         Beware that we can't use the superclass method, because it
 470         relies on the algebra being associative.
 471
 472         SETUP::
 473
 474             sage: from mjo.eja.eja_algebra import random_eja
 475
 476         TESTS:
 477
 478         The identity element is always invertible::
 479
 480             sage: set_random_seed()
 481             sage: J = random_eja()
 482             sage: J.one().is_invertible()
 483             True
 484
 485         The zero element is never invertible in a non-trivial algebra::
 486
 487             sage: set_random_seed()
 488             sage: J = random_eja()
 489             sage: (not J.is_trivial()) and J.zero().is_invertible()
 490             False
 491
 492         """
 493         if self.is_zero():
 494             if self.parent().is_trivial():
 495                 return True
 496             else:
 497                 return False
 498
 499         # In fact, we only need to know if the constant term is non-zero,
 500         # so we can pass in the field's zero element instead.
 501         zero = self.base_ring().zero()
 502         p = self.minimal_polynomial()
 503         return not (p(zero) == zero)
 504
 505
 506     def is_nilpotent(self):
 507         """
 508         Return whether or not some power of this element is zero.
 509
 510         ALGORITHM:
 511
 512         We use Theorem 5 in Chapter III of Koecher, which says that
 513         an element ``x`` is nilpotent if and only if ``x.operator()``
 514         is nilpotent. And it is a basic fact of linear algebra that
 515         an operator on an `n`-dimensional space is nilpotent if and
 516         only if, when raised to the `n`th power, it equals the zero
 517         operator (for example, see Axler Corollary 8.8).
 518
 519         SETUP::
 520
 521             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 522             ....:                                  random_eja)
 523
 524         EXAMPLES::
 525
 526             sage: J = JordanSpinEJA(3)
 527             sage: x = sum(J.gens())
 528             sage: x.is_nilpotent()
 529             False
 530
 531         TESTS:
 532
 533         The identity element is never nilpotent::
 534
 535             sage: set_random_seed()
 536             sage: random_eja().one().is_nilpotent()
 537             False
 538
 539         The additive identity is always nilpotent::
 540
 541             sage: set_random_seed()
 542             sage: random_eja().zero().is_nilpotent()
 543             True
 544
 545         """
 546         P = self.parent()
 547         zero_operator = P.zero().operator()
 548         return self.operator()**P.dimension() == zero_operator
 549
 550
 551     def is_regular(self):
 552         """
 553         Return whether or not this is a regular element.
 554
 555         SETUP::
 556
 557             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 558             ....:                                  random_eja)
 559
 560         EXAMPLES:
 561
 562         The identity element always has degree one, but any element
 563         linearly-independent from it is regular::
 564
 565             sage: J = JordanSpinEJA(5)
 566             sage: J.one().is_regular()
 567             False
 568             sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 569             sage: for x in J.gens():
 570             ....:     (J.one() + x).is_regular()
 571             False
 572             True
 573             True
 574             True
 575             True
 576
 577         TESTS:
 578
 579         The zero element should never be regular, unless the parent
 580         algebra has dimension one::
 581
 582             sage: set_random_seed()
 583             sage: J = random_eja()
 584             sage: J.dimension() == 1 or not J.zero().is_regular()
 585             True
 586
 587         The unit element isn't regular unless the algebra happens to
 588         consist of only its scalar multiples::
 589
 590             sage: set_random_seed()
 591             sage: J = random_eja()
 592             sage: J.dimension() == 1 or not J.one().is_regular()
 593             True
 594
 595         """
 596         return self.degree() == self.parent().rank()
 597
 598
 599     def degree(self):
 600         """
 601         Return the degree of this element, which is defined to be
 602         the degree of its minimal polynomial.
 603
 604         ALGORITHM:
 605
 606         For now, we skip the messy minimal polynomial computation
 607         and instead return the dimension of the vector space spanned
 608         by the powers of this element. The latter is a bit more
 609         straightforward to compute.
 610
 611         SETUP::
 612
 613             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 614             ....:                                  random_eja)
 615
 616         EXAMPLES::
 617
 618             sage: J = JordanSpinEJA(4)
 619             sage: J.one().degree()
 620             1
 621             sage: e0,e1,e2,e3 = J.gens()
 622             sage: (e0 - e1).degree()
 623             2
 624
 625         In the spin factor algebra (of rank two), all elements that
 626         aren't multiples of the identity are regular::
 627
 628             sage: set_random_seed()
 629             sage: n = ZZ.random_element(1,10)
 630             sage: J = JordanSpinEJA(n)
 631             sage: x = J.random_element()
 632             sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 633             True
 634
 635         TESTS:
 636
 637         The zero and unit elements are both of degree one::
 638
 639             sage: set_random_seed()
 640             sage: J = random_eja()
 641             sage: J.zero().degree()
 642             1
 643             sage: J.one().degree()
 644             1
 645
 646         Our implementation agrees with the definition::
 647
 648             sage: set_random_seed()
 649             sage: x = random_eja().random_element()
 650             sage: x.degree() == x.minimal_polynomial().degree()
 651             True
 652
 653         """
 654         if self.is_zero() and not self.parent().is_trivial():
 655             # The minimal polynomial of zero in a nontrivial algebra
 656             # is "t"; in a trivial algebra it's "1" by convention
 657             # (it's an empty product).
 658             return 1
 659         return self.subalgebra_generated_by().dimension()
 660
 661
 662     def left_matrix(self):
 663         """
 664         Our parent class defines ``left_matrix`` and ``matrix``
 665         methods whose names are misleading. We don't want them.
 666         """
 667         raise NotImplementedError("use operator().matrix() instead")
 668
 669     matrix = left_matrix
 670
 671
 672     def minimal_polynomial(self):
 673         """
 674         Return the minimal polynomial of this element,
 675         as a function of the variable `t`.
 676
 677         ALGORITHM:
 678
 679         We restrict ourselves to the associative subalgebra
 680         generated by this element, and then return the minimal
 681         polynomial of this element's operator matrix (in that
 682         subalgebra). This works by Baes Proposition 2.3.16.
 683
 684         SETUP::
 685
 686             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 687             ....:                                  random_eja)
 688
 689         TESTS:
 690
 691         The minimal polynomial of the identity and zero elements are
 692         always the same::
 693
 694             sage: set_random_seed()
 695             sage: J = random_eja()
 696             sage: J.one().minimal_polynomial()
 697             t - 1
 698             sage: J.zero().minimal_polynomial()
 699             t
 700
 701         The degree of an element is (by one definition) the degree
 702         of its minimal polynomial::
 703
 704             sage: set_random_seed()
 705             sage: x = random_eja().random_element()
 706             sage: x.degree() == x.minimal_polynomial().degree()
 707             True
 708
 709         The minimal polynomial and the characteristic polynomial coincide
 710         and are known (see Alizadeh, Example 11.11) for all elements of
 711         the spin factor algebra that aren't scalar multiples of the
 712         identity::
 713
 714             sage: set_random_seed()
 715             sage: n = ZZ.random_element(2,10)
 716             sage: J = JordanSpinEJA(n)
 717             sage: y = J.random_element()
 718             sage: while y == y.coefficient(0)*J.one():
 719             ....:     y = J.random_element()
 720             sage: y0 = y.to_vector()[0]
 721             sage: y_bar = y.to_vector()[1:]
 722             sage: actual = y.minimal_polynomial()
 723             sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
 724             sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
 725             sage: bool(actual == expected)
 726             True
 727
 728         The minimal polynomial should always kill its element::
 729
 730             sage: set_random_seed()
 731             sage: x = random_eja().random_element()
 732             sage: p = x.minimal_polynomial()
 733             sage: x.apply_univariate_polynomial(p)
 734             0
 735
 736         """
 737         if self.is_zero():
 738             # We would generate a zero-dimensional subalgebra
 739             # where the minimal polynomial would be constant.
 740             # That might be correct, but only if *this* algebra
 741             # is trivial too.
 742             if not self.parent().is_trivial():
 743                 # Pretty sure we know what the minimal polynomial of
 744                 # the zero operator is going to be. This ensures
 745                 # consistency of e.g. the polynomial variable returned
 746                 # in the "normal" case without us having to think about it.
 747                 return self.operator().minimal_polynomial()
 748
 749         A = self.subalgebra_generated_by()
 750         return A(self).operator().minimal_polynomial()
 751
 752
 753
 754     def natural_representation(self):
 755         """
 756         Return a more-natural representation of this element.
 757
 758         Every finite-dimensional Euclidean Jordan Algebra is a
 759         direct sum of five simple algebras, four of which comprise
 760         Hermitian matrices. This method returns the original
 761         "natural" representation of this element as a Hermitian
 762         matrix, if it has one. If not, you get the usual representation.
 763
 764         SETUP::
 765
 766             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 767             ....:                                  QuaternionHermitianEJA)
 768
 769         EXAMPLES::
 770
 771             sage: J = ComplexHermitianEJA(3)
 772             sage: J.one()
 773             e0 + e3 + e8
 774             sage: J.one().natural_representation()
 775             [1 0 0 0 0 0]
 776             [0 1 0 0 0 0]
 777             [0 0 1 0 0 0]
 778             [0 0 0 1 0 0]
 779             [0 0 0 0 1 0]
 780             [0 0 0 0 0 1]
 781
 782         ::
 783
 784             sage: J = QuaternionHermitianEJA(3)
 785             sage: J.one()
 786             e0 + e5 + e14
 787             sage: J.one().natural_representation()
 788             [1 0 0 0 0 0 0 0 0 0 0 0]
 789             [0 1 0 0 0 0 0 0 0 0 0 0]
 790             [0 0 1 0 0 0 0 0 0 0 0 0]
 791             [0 0 0 1 0 0 0 0 0 0 0 0]
 792             [0 0 0 0 1 0 0 0 0 0 0 0]
 793             [0 0 0 0 0 1 0 0 0 0 0 0]
 794             [0 0 0 0 0 0 1 0 0 0 0 0]
 795             [0 0 0 0 0 0 0 1 0 0 0 0]
 796             [0 0 0 0 0 0 0 0 1 0 0 0]
 797             [0 0 0 0 0 0 0 0 0 1 0 0]
 798             [0 0 0 0 0 0 0 0 0 0 1 0]
 799             [0 0 0 0 0 0 0 0 0 0 0 1]
 800
 801         """
 802         B = self.parent().natural_basis()
 803         W = B[0].matrix_space()
 804         return W.linear_combination(zip(B,self.to_vector()))
 805
 806
 807     def operator(self):
 808         """
 809         Return the left-multiplication-by-this-element
 810         operator on the ambient algebra.
 811
 812         SETUP::
 813
 814             sage: from mjo.eja.eja_algebra import random_eja
 815
 816         TESTS::
 817
 818             sage: set_random_seed()
 819             sage: J = random_eja()
 820             sage: x = J.random_element()
 821             sage: y = J.random_element()
 822             sage: x.operator()(y) == x*y
 823             True
 824             sage: y.operator()(x) == x*y
 825             True
 826
 827         """
 828         P = self.parent()
 829         left_mult_by_self = lambda y: self*y
 830         L = P.module_morphism(function=left_mult_by_self, codomain=P)
 831         return FiniteDimensionalEuclideanJordanAlgebraOperator(
 832                  P,
 833                  P,
 834                  L.matrix() )
 835
 836
 837     def quadratic_representation(self, other=None):
 838         """
 839         Return the quadratic representation of this element.
 840
 841         SETUP::
 842
 843             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 844             ....:                                  random_eja)
 845
 846         EXAMPLES:
 847
 848         The explicit form in the spin factor algebra is given by
 849         Alizadeh's Example 11.12::
 850
 851             sage: set_random_seed()
 852             sage: n = ZZ.random_element(1,10)
 853             sage: J = JordanSpinEJA(n)
 854             sage: x = J.random_element()
 855             sage: x_vec = x.to_vector()
 856             sage: x0 = x_vec[0]
 857             sage: x_bar = x_vec[1:]
 858             sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 859             sage: B = 2*x0*x_bar.row()
 860             sage: C = 2*x0*x_bar.column()
 861             sage: D = matrix.identity(QQ, n-1)
 862             sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 863             sage: D = D + 2*x_bar.tensor_product(x_bar)
 864             sage: Q = matrix.block(2,2,[A,B,C,D])
 865             sage: Q == x.quadratic_representation().matrix()
 866             True
 867
 868         Test all of the properties from Theorem 11.2 in Alizadeh::
 869
 870             sage: set_random_seed()
 871             sage: J = random_eja()
 872             sage: x = J.random_element()
 873             sage: y = J.random_element()
 874             sage: Lx = x.operator()
 875             sage: Lxx = (x*x).operator()
 876             sage: Qx = x.quadratic_representation()
 877             sage: Qy = y.quadratic_representation()
 878             sage: Qxy = x.quadratic_representation(y)
 879             sage: Qex = J.one().quadratic_representation(x)
 880             sage: n = ZZ.random_element(10)
 881             sage: Qxn = (x^n).quadratic_representation()
 882
 883         Property 1:
 884
 885             sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
 886             True
 887
 888         Property 2 (multiply on the right for :trac:`28272`):
 889
 890             sage: alpha = QQ.random_element()
 891             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
 892             True
 893
 894         Property 3:
 895
 896             sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
 897             True
 898
 899             sage: not x.is_invertible() or (
 900             ....:   ~Qx
 901             ....:   ==
 902             ....:   x.inverse().quadratic_representation() )
 903             True
 904
 905             sage: Qxy(J.one()) == x*y
 906             True
 907
 908         Property 4:
 909
 910             sage: not x.is_invertible() or (
 911             ....:   x.quadratic_representation(x.inverse())*Qx
 912             ....:   == Qx*x.quadratic_representation(x.inverse()) )
 913             True
 914
 915             sage: not x.is_invertible() or (
 916             ....:   x.quadratic_representation(x.inverse())*Qx
 917             ....:   ==
 918             ....:   2*x.operator()*Qex - Qx )
 919             True
 920
 921             sage: 2*x.operator()*Qex - Qx == Lxx
 922             True
 923
 924         Property 5:
 925
 926             sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
 927             True
 928
 929         Property 6:
 930
 931             sage: Qxn == (Qx)^n
 932             True
 933
 934         Property 7:
 935
 936             sage: not x.is_invertible() or (
 937             ....:   Qx*x.inverse().operator() == Lx )
 938             True
 939
 940         Property 8:
 941
 942             sage: not x.operator_commutes_with(y) or (
 943             ....:   Qx(y)^n == Qxn(y^n) )
 944             True
 945
 946         """
 947         if other is None:
 948             other=self
 949         elif not other in self.parent():
 950             raise TypeError("'other' must live in the same algebra")
 951
 952         L = self.operator()
 953         M = other.operator()
 954         return ( L*M + M*L - (self*other).operator() )
 955
 956
 957
 958
 959     def subalgebra_generated_by(self):
 960         """
 961         Return the associative subalgebra of the parent EJA generated
 962         by this element.
 963
 964         SETUP::
 965
 966             sage: from mjo.eja.eja_algebra import random_eja
 967
 968         TESTS:
 969
 970         This subalgebra, being composed of only powers, is associative::
 971
 972             sage: set_random_seed()
 973             sage: x0 = random_eja().random_element()
 974             sage: A = x0.subalgebra_generated_by()
 975             sage: x = A.random_element()
 976             sage: y = A.random_element()
 977             sage: z = A.random_element()
 978             sage: (x*y)*z == x*(y*z)
 979             True
 980
 981         Squaring in the subalgebra should work the same as in
 982         the superalgebra::
 983
 984             sage: set_random_seed()
 985             sage: x = random_eja().random_element()
 986             sage: A = x.subalgebra_generated_by()
 987             sage: A(x^2) == A(x)*A(x)
 988             True
 989
 990         The subalgebra generated by the zero element is trivial::
 991
 992             sage: set_random_seed()
 993             sage: A = random_eja().zero().subalgebra_generated_by()
 994             sage: A
 995             Euclidean Jordan algebra of dimension 0 over Rational Field
 996             sage: A.one()
 997             0
 998
 999         """
1000         return FiniteDimensionalEuclideanJordanElementSubalgebra(self)
1001
1002
1003     def subalgebra_idempotent(self):
1004         """
1005         Find an idempotent in the associative subalgebra I generate
1006         using Proposition 2.3.5 in Baes.
1007
1008         SETUP::
1009
1010             sage: from mjo.eja.eja_algebra import random_eja
1011
1012         TESTS::
1013
1014             sage: set_random_seed()
1015             sage: J = random_eja()
1016             sage: x = J.random_element()
1017             sage: while x.is_nilpotent():
1018             ....:     x = J.random_element()
1019             sage: c = x.subalgebra_idempotent()
1020             sage: c^2 == c
1021             True
1022
1023         """
1024         if self.is_nilpotent():
1025             raise ValueError("this only works with non-nilpotent elements!")
1026
1027         J = self.subalgebra_generated_by()
1028         u = J(self)
1029
1030         # The image of the matrix of left-u^m-multiplication
1031         # will be minimal for some natural number s...
1032         s = 0
1033         minimal_dim = J.dimension()
1034         for i in xrange(1, minimal_dim):
1035             this_dim = (u**i).operator().matrix().image().dimension()
1036             if this_dim < minimal_dim:
1037                 minimal_dim = this_dim
1038                 s = i
1039
1040         # Now minimal_matrix should correspond to the smallest
1041         # non-zero subspace in Baes's (or really, Koecher's)
1042         # proposition.
1043         #
1044         # However, we need to restrict the matrix to work on the
1045         # subspace... or do we? Can't we just solve, knowing that
1046         # A(c) = u^(s+1) should have a solution in the big space,
1047         # too?
1048         #
1049         # Beware, solve_right() means that we're using COLUMN vectors.
1050         # Our FiniteDimensionalAlgebraElement superclass uses rows.
1051         u_next = u**(s+1)
1052         A = u_next.operator().matrix()
1053         c = J.from_vector(A.solve_right(u_next.to_vector()))
1054
1055         # Now c is the idempotent we want, but it still lives in the subalgebra.
1056         return c.superalgebra_element()
1057
1058
1059     def trace(self):
1060         """
1061         Return my trace, the sum of my eigenvalues.
1062
1063         SETUP::
1064
1065             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1066             ....:                                  RealCartesianProductEJA,
1067             ....:                                  random_eja)
1068
1069         EXAMPLES::
1070
1071             sage: J = JordanSpinEJA(3)
1072             sage: x = sum(J.gens())
1073             sage: x.trace()
1074             2
1075
1076         ::
1077
1078             sage: J = RealCartesianProductEJA(5)
1079             sage: J.one().trace()
1080             5
1081
1082         TESTS:
1083
1084         The trace of an element is a real number::
1085
1086             sage: set_random_seed()
1087             sage: J = random_eja()
1088             sage: J.random_element().trace() in J.base_ring()
1089             True
1090
1091         """
1092         P = self.parent()
1093         r = P.rank()
1094         p = P._charpoly_coeff(r-1)
1095         # The _charpoly_coeff function already adds the factor of
1096         # -1 to ensure that _charpoly_coeff(r-1) is really what
1097         # appears in front of t^{r-1} in the charpoly. However,
1098         # we want the negative of THAT for the trace.
1099         return -p(*self.to_vector())
1100
1101
1102     def trace_inner_product(self, other):
1103         """
1104         Return the trace inner product of myself and ``other``.
1105
1106         SETUP::
1107
1108             sage: from mjo.eja.eja_algebra import random_eja
1109
1110         TESTS:
1111
1112         The trace inner product is commutative::
1113
1114             sage: set_random_seed()
1115             sage: J = random_eja()
1116             sage: x = J.random_element(); y = J.random_element()
1117             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1118             True
1119
1120         The trace inner product is bilinear::
1121
1122             sage: set_random_seed()
1123             sage: J = random_eja()
1124             sage: x = J.random_element()
1125             sage: y = J.random_element()
1126             sage: z = J.random_element()
1127             sage: a = QQ.random_element();
1128             sage: actual = (a*(x+z)).trace_inner_product(y)
1129             sage: expected = ( a*x.trace_inner_product(y) +
1130             ....:              a*z.trace_inner_product(y) )
1131             sage: actual == expected
1132             True
1133             sage: actual = x.trace_inner_product(a*(y+z))
1134             sage: expected = ( a*x.trace_inner_product(y) +
1135             ....:              a*x.trace_inner_product(z) )
1136             sage: actual == expected
1137             True
1138
1139         The trace inner product satisfies the compatibility
1140         condition in the definition of a Euclidean Jordan algebra::
1141
1142             sage: set_random_seed()
1143             sage: J = random_eja()
1144             sage: x = J.random_element()
1145             sage: y = J.random_element()
1146             sage: z = J.random_element()
1147             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1148             True
1149
1150         """
1151         if not other in self.parent():
1152             raise TypeError("'other' must live in the same algebra")
1153
1154         return (self*other).trace()