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