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(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::
 486
 487             sage: set_random_seed()
 488             sage: J = random_eja()
 489             sage: J.zero().is_invertible()
 490             False
 491
 492         """
 493         zero = self.parent().zero()
 494         p = self.minimal_polynomial()
 495         return not (p(zero) == zero)
 496
 497
 498     def is_nilpotent(self):
 499         """
 500         Return whether or not some power of this element is zero.
 501
 502         ALGORITHM:
 503
 504         We use Theorem 5 in Chapter III of Koecher, which says that
 505         an element ``x`` is nilpotent if and only if ``x.operator()``
 506         is nilpotent. And it is a basic fact of linear algebra that
 507         an operator on an `n`-dimensional space is nilpotent if and
 508         only if, when raised to the `n`th power, it equals the zero
 509         operator (for example, see Axler Corollary 8.8).
 510
 511         SETUP::
 512
 513             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 514             ....:                                  random_eja)
 515
 516         EXAMPLES::
 517
 518             sage: J = JordanSpinEJA(3)
 519             sage: x = sum(J.gens())
 520             sage: x.is_nilpotent()
 521             False
 522
 523         TESTS:
 524
 525         The identity element is never nilpotent::
 526
 527             sage: set_random_seed()
 528             sage: random_eja().one().is_nilpotent()
 529             False
 530
 531         The additive identity is always nilpotent::
 532
 533             sage: set_random_seed()
 534             sage: random_eja().zero().is_nilpotent()
 535             True
 536
 537         """
 538         P = self.parent()
 539         zero_operator = P.zero().operator()
 540         return self.operator()**P.dimension() == zero_operator
 541
 542
 543     def is_regular(self):
 544         """
 545         Return whether or not this is a regular element.
 546
 547         SETUP::
 548
 549             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 550             ....:                                  random_eja)
 551
 552         EXAMPLES:
 553
 554         The identity element always has degree one, but any element
 555         linearly-independent from it is regular::
 556
 557             sage: J = JordanSpinEJA(5)
 558             sage: J.one().is_regular()
 559             False
 560             sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 561             sage: for x in J.gens():
 562             ....:     (J.one() + x).is_regular()
 563             False
 564             True
 565             True
 566             True
 567             True
 568
 569         TESTS:
 570
 571         The zero element should never be regular, unless the parent
 572         algebra has dimension one::
 573
 574             sage: set_random_seed()
 575             sage: J = random_eja()
 576             sage: J.dimension() == 1 or not J.zero().is_regular()
 577             True
 578
 579         The unit element isn't regular unless the algebra happens to
 580         consist of only its scalar multiples::
 581
 582             sage: set_random_seed()
 583             sage: J = random_eja()
 584             sage: J.dimension() == 1 or not J.one().is_regular()
 585             True
 586
 587         """
 588         return self.degree() == self.parent().rank()
 589
 590
 591     def degree(self):
 592         """
 593         Return the degree of this element, which is defined to be
 594         the degree of its minimal polynomial.
 595
 596         ALGORITHM:
 597
 598         For now, we skip the messy minimal polynomial computation
 599         and instead return the dimension of the vector space spanned
 600         by the powers of this element. The latter is a bit more
 601         straightforward to compute.
 602
 603         SETUP::
 604
 605             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 606             ....:                                  random_eja)
 607
 608         EXAMPLES::
 609
 610             sage: J = JordanSpinEJA(4)
 611             sage: J.one().degree()
 612             1
 613             sage: e0,e1,e2,e3 = J.gens()
 614             sage: (e0 - e1).degree()
 615             2
 616
 617         In the spin factor algebra (of rank two), all elements that
 618         aren't multiples of the identity are regular::
 619
 620             sage: set_random_seed()
 621             sage: n = ZZ.random_element(1,10)
 622             sage: J = JordanSpinEJA(n)
 623             sage: x = J.random_element()
 624             sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 625             True
 626
 627         TESTS:
 628
 629         The zero and unit elements are both of degree one::
 630
 631             sage: set_random_seed()
 632             sage: J = random_eja()
 633             sage: J.zero().degree()
 634             1
 635             sage: J.one().degree()
 636             1
 637
 638         Our implementation agrees with the definition::
 639
 640             sage: set_random_seed()
 641             sage: x = random_eja().random_element()
 642             sage: x.degree() == x.minimal_polynomial().degree()
 643             True
 644
 645         """
 646         return self.subalgebra_generated_by().dimension()
 647
 648
 649     def left_matrix(self):
 650         """
 651         Our parent class defines ``left_matrix`` and ``matrix``
 652         methods whose names are misleading. We don't want them.
 653         """
 654         raise NotImplementedError("use operator().matrix() instead")
 655
 656     matrix = left_matrix
 657
 658
 659     def minimal_polynomial(self):
 660         """
 661         Return the minimal polynomial of this element,
 662         as a function of the variable `t`.
 663
 664         ALGORITHM:
 665
 666         We restrict ourselves to the associative subalgebra
 667         generated by this element, and then return the minimal
 668         polynomial of this element's operator matrix (in that
 669         subalgebra). This works by Baes Proposition 2.3.16.
 670
 671         SETUP::
 672
 673             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 674             ....:                                  random_eja)
 675
 676         TESTS:
 677
 678         The minimal polynomial of the identity and zero elements are
 679         always the same::
 680
 681             sage: set_random_seed()
 682             sage: J = random_eja()
 683             sage: J.one().minimal_polynomial()
 684             t - 1
 685             sage: J.zero().minimal_polynomial()
 686             t
 687
 688         The degree of an element is (by one definition) the degree
 689         of its minimal polynomial::
 690
 691             sage: set_random_seed()
 692             sage: x = random_eja().random_element()
 693             sage: x.degree() == x.minimal_polynomial().degree()
 694             True
 695
 696         The minimal polynomial and the characteristic polynomial coincide
 697         and are known (see Alizadeh, Example 11.11) for all elements of
 698         the spin factor algebra that aren't scalar multiples of the
 699         identity::
 700
 701             sage: set_random_seed()
 702             sage: n = ZZ.random_element(2,10)
 703             sage: J = JordanSpinEJA(n)
 704             sage: y = J.random_element()
 705             sage: while y == y.coefficient(0)*J.one():
 706             ....:     y = J.random_element()
 707             sage: y0 = y.to_vector()[0]
 708             sage: y_bar = y.to_vector()[1:]
 709             sage: actual = y.minimal_polynomial()
 710             sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
 711             sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
 712             sage: bool(actual == expected)
 713             True
 714
 715         The minimal polynomial should always kill its element::
 716
 717             sage: set_random_seed()
 718             sage: x = random_eja().random_element()
 719             sage: p = x.minimal_polynomial()
 720             sage: x.apply_univariate_polynomial(p)
 721             0
 722
 723         """
 724         A = self.subalgebra_generated_by()
 725         return A.element_class(A,self).operator().minimal_polynomial()
 726
 727
 728
 729     def natural_representation(self):
 730         """
 731         Return a more-natural representation of this element.
 732
 733         Every finite-dimensional Euclidean Jordan Algebra is a
 734         direct sum of five simple algebras, four of which comprise
 735         Hermitian matrices. This method returns the original
 736         "natural" representation of this element as a Hermitian
 737         matrix, if it has one. If not, you get the usual representation.
 738
 739         SETUP::
 740
 741             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 742             ....:                                  QuaternionHermitianEJA)
 743
 744         EXAMPLES::
 745
 746             sage: J = ComplexHermitianEJA(3)
 747             sage: J.one()
 748             e0 + e3 + e8
 749             sage: J.one().natural_representation()
 750             [1 0 0 0 0 0]
 751             [0 1 0 0 0 0]
 752             [0 0 1 0 0 0]
 753             [0 0 0 1 0 0]
 754             [0 0 0 0 1 0]
 755             [0 0 0 0 0 1]
 756
 757         ::
 758
 759             sage: J = QuaternionHermitianEJA(3)
 760             sage: J.one()
 761             e0 + e5 + e14
 762             sage: J.one().natural_representation()
 763             [1 0 0 0 0 0 0 0 0 0 0 0]
 764             [0 1 0 0 0 0 0 0 0 0 0 0]
 765             [0 0 1 0 0 0 0 0 0 0 0 0]
 766             [0 0 0 1 0 0 0 0 0 0 0 0]
 767             [0 0 0 0 1 0 0 0 0 0 0 0]
 768             [0 0 0 0 0 1 0 0 0 0 0 0]
 769             [0 0 0 0 0 0 1 0 0 0 0 0]
 770             [0 0 0 0 0 0 0 1 0 0 0 0]
 771             [0 0 0 0 0 0 0 0 1 0 0 0]
 772             [0 0 0 0 0 0 0 0 0 1 0 0]
 773             [0 0 0 0 0 0 0 0 0 0 1 0]
 774             [0 0 0 0 0 0 0 0 0 0 0 1]
 775
 776         """
 777         B = self.parent().natural_basis()
 778         W = B[0].matrix_space()
 779         return W.linear_combination(zip(B,self.to_vector()))
 780
 781
 782     def operator(self):
 783         """
 784         Return the left-multiplication-by-this-element
 785         operator on the ambient algebra.
 786
 787         SETUP::
 788
 789             sage: from mjo.eja.eja_algebra import random_eja
 790
 791         TESTS::
 792
 793             sage: set_random_seed()
 794             sage: J = random_eja()
 795             sage: x = J.random_element()
 796             sage: y = J.random_element()
 797             sage: x.operator()(y) == x*y
 798             True
 799             sage: y.operator()(x) == x*y
 800             True
 801
 802         """
 803         P = self.parent()
 804         return FiniteDimensionalEuclideanJordanAlgebraOperator(
 805                  P,
 806                  P,
 807                  self.to_matrix() )
 808
 809
 810     def quadratic_representation(self, other=None):
 811         """
 812         Return the quadratic representation of this element.
 813
 814         SETUP::
 815
 816             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 817             ....:                                  random_eja)
 818
 819         EXAMPLES:
 820
 821         The explicit form in the spin factor algebra is given by
 822         Alizadeh's Example 11.12::
 823
 824             sage: set_random_seed()
 825             sage: n = ZZ.random_element(1,10)
 826             sage: J = JordanSpinEJA(n)
 827             sage: x = J.random_element()
 828             sage: x_vec = x.to_vector()
 829             sage: x0 = x_vec[0]
 830             sage: x_bar = x_vec[1:]
 831             sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 832             sage: B = 2*x0*x_bar.row()
 833             sage: C = 2*x0*x_bar.column()
 834             sage: D = matrix.identity(QQ, n-1)
 835             sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 836             sage: D = D + 2*x_bar.tensor_product(x_bar)
 837             sage: Q = matrix.block(2,2,[A,B,C,D])
 838             sage: Q == x.quadratic_representation().matrix()
 839             True
 840
 841         Test all of the properties from Theorem 11.2 in Alizadeh::
 842
 843             sage: set_random_seed()
 844             sage: J = random_eja()
 845             sage: x = J.random_element()
 846             sage: y = J.random_element()
 847             sage: Lx = x.operator()
 848             sage: Lxx = (x*x).operator()
 849             sage: Qx = x.quadratic_representation()
 850             sage: Qy = y.quadratic_representation()
 851             sage: Qxy = x.quadratic_representation(y)
 852             sage: Qex = J.one().quadratic_representation(x)
 853             sage: n = ZZ.random_element(10)
 854             sage: Qxn = (x^n).quadratic_representation()
 855
 856         Property 1:
 857
 858             sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
 859             True
 860
 861         Property 2 (multiply on the right for :trac:`28272`):
 862
 863             sage: alpha = QQ.random_element()
 864             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
 865             True
 866
 867         Property 3:
 868
 869             sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
 870             True
 871
 872             sage: not x.is_invertible() or (
 873             ....:   ~Qx
 874             ....:   ==
 875             ....:   x.inverse().quadratic_representation() )
 876             True
 877
 878             sage: Qxy(J.one()) == x*y
 879             True
 880
 881         Property 4:
 882
 883             sage: not x.is_invertible() or (
 884             ....:   x.quadratic_representation(x.inverse())*Qx
 885             ....:   == Qx*x.quadratic_representation(x.inverse()) )
 886             True
 887
 888             sage: not x.is_invertible() or (
 889             ....:   x.quadratic_representation(x.inverse())*Qx
 890             ....:   ==
 891             ....:   2*x.operator()*Qex - Qx )
 892             True
 893
 894             sage: 2*x.operator()*Qex - Qx == Lxx
 895             True
 896
 897         Property 5:
 898
 899             sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
 900             True
 901
 902         Property 6:
 903
 904             sage: Qxn == (Qx)^n
 905             True
 906
 907         Property 7:
 908
 909             sage: not x.is_invertible() or (
 910             ....:   Qx*x.inverse().operator() == Lx )
 911             True
 912
 913         Property 8:
 914
 915             sage: not x.operator_commutes_with(y) or (
 916             ....:   Qx(y)^n == Qxn(y^n) )
 917             True
 918
 919         """
 920         if other is None:
 921             other=self
 922         elif not other in self.parent():
 923             raise TypeError("'other' must live in the same algebra")
 924
 925         L = self.operator()
 926         M = other.operator()
 927         return ( L*M + M*L - (self*other).operator() )
 928
 929
 930
 931
 932     def subalgebra_generated_by(self):
 933         """
 934         Return the associative subalgebra of the parent EJA generated
 935         by this element.
 936
 937         SETUP::
 938
 939             sage: from mjo.eja.eja_algebra import random_eja
 940
 941         TESTS::
 942
 943             sage: set_random_seed()
 944             sage: x = random_eja().random_element()
 945             sage: x.subalgebra_generated_by().is_associative()
 946             True
 947
 948         Squaring in the subalgebra should work the same as in
 949         the superalgebra::
 950
 951             sage: set_random_seed()
 952             sage: x = random_eja().random_element()
 953             sage: A = x.subalgebra_generated_by()
 954             sage: A(x^2) == A(x)*A(x)
 955             True
 956
 957         """
 958         return FiniteDimensionalEuclideanJordanElementSubalgebra(self)
 959
 960
 961     def subalgebra_idempotent(self):
 962         """
 963         Find an idempotent in the associative subalgebra I generate
 964         using Proposition 2.3.5 in Baes.
 965
 966         SETUP::
 967
 968             sage: from mjo.eja.eja_algebra import random_eja
 969
 970         TESTS::
 971
 972             sage: set_random_seed()
 973             sage: J = random_eja()
 974             sage: x = J.random_element()
 975             sage: while x.is_nilpotent():
 976             ....:     x = J.random_element()
 977             sage: c = x.subalgebra_idempotent()
 978             sage: c^2 == c
 979             True
 980
 981         """
 982         if self.is_nilpotent():
 983             raise ValueError("this only works with non-nilpotent elements!")
 984
 985         J = self.subalgebra_generated_by()
 986         u = J.from_vector(self.to_vector())
 987
 988         # The image of the matrix of left-u^m-multiplication
 989         # will be minimal for some natural number s...
 990         s = 0
 991         minimal_dim = J.dimension()
 992         for i in xrange(1, minimal_dim):
 993             this_dim = (u**i).operator().matrix().image().dimension()
 994             if this_dim < minimal_dim:
 995                 minimal_dim = this_dim
 996                 s = i
 997
 998         # Now minimal_matrix should correspond to the smallest
 999         # non-zero subspace in Baes's (or really, Koecher's)
1000         # proposition.
1001         #
1002         # However, we need to restrict the matrix to work on the
1003         # subspace... or do we? Can't we just solve, knowing that
1004         # A(c) = u^(s+1) should have a solution in the big space,
1005         # too?
1006         #
1007         # Beware, solve_right() means that we're using COLUMN vectors.
1008         # Our FiniteDimensionalAlgebraElement superclass uses rows.
1009         u_next = u**(s+1)
1010         A = u_next.operator().matrix()
1011         c = J(A.solve_right(u_next.to_vector()))
1012
1013         # Now c is the idempotent we want, but it still lives in the subalgebra.
1014         return c.superalgebra_element()
1015
1016
1017     def trace(self):
1018         """
1019         Return my trace, the sum of my eigenvalues.
1020
1021         SETUP::
1022
1023             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1024             ....:                                  RealCartesianProductEJA,
1025             ....:                                  random_eja)
1026
1027         EXAMPLES::
1028
1029             sage: J = JordanSpinEJA(3)
1030             sage: x = sum(J.gens())
1031             sage: x.trace()
1032             2
1033
1034         ::
1035
1036             sage: J = RealCartesianProductEJA(5)
1037             sage: J.one().trace()
1038             5
1039
1040         TESTS:
1041
1042         The trace of an element is a real number::
1043
1044             sage: set_random_seed()
1045             sage: J = random_eja()
1046             sage: J.random_element().trace() in J.base_ring()
1047             True
1048
1049         """
1050         P = self.parent()
1051         r = P.rank()
1052         p = P._charpoly_coeff(r-1)
1053         # The _charpoly_coeff function already adds the factor of
1054         # -1 to ensure that _charpoly_coeff(r-1) is really what
1055         # appears in front of t^{r-1} in the charpoly. However,
1056         # we want the negative of THAT for the trace.
1057         return -p(*self.to_vector())
1058
1059
1060     def trace_inner_product(self, other):
1061         """
1062         Return the trace inner product of myself and ``other``.
1063
1064         SETUP::
1065
1066             sage: from mjo.eja.eja_algebra import random_eja
1067
1068         TESTS:
1069
1070         The trace inner product is commutative::
1071
1072             sage: set_random_seed()
1073             sage: J = random_eja()
1074             sage: x = J.random_element(); y = J.random_element()
1075             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1076             True
1077
1078         The trace inner product is bilinear::
1079
1080             sage: set_random_seed()
1081             sage: J = random_eja()
1082             sage: x = J.random_element()
1083             sage: y = J.random_element()
1084             sage: z = J.random_element()
1085             sage: a = QQ.random_element();
1086             sage: actual = (a*(x+z)).trace_inner_product(y)
1087             sage: expected = ( a*x.trace_inner_product(y) +
1088             ....:              a*z.trace_inner_product(y) )
1089             sage: actual == expected
1090             True
1091             sage: actual = x.trace_inner_product(a*(y+z))
1092             sage: expected = ( a*x.trace_inner_product(y) +
1093             ....:              a*x.trace_inner_product(z) )
1094             sage: actual == expected
1095             True
1096
1097         The trace inner product satisfies the compatibility
1098         condition in the definition of a Euclidean Jordan algebra::
1099
1100             sage: set_random_seed()
1101             sage: J = random_eja()
1102             sage: x = J.random_element()
1103             sage: y = J.random_element()
1104             sage: z = J.random_element()
1105             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1106             True
1107
1108         """
1109         if not other in self.parent():
1110             raise TypeError("'other' must live in the same algebra")
1111
1112         return (self*other).trace()