mjo/eja/eja_element.py

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