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, unless the parent
 631         algebra has dimension one::
 632
 633             sage: set_random_seed()
 634             sage: J = random_eja()
 635             sage: J.dimension() == 1 or not J.zero().is_regular()
 636             True
 637
 638         The unit element isn't regular unless the algebra happens to
 639         consist of only its scalar multiples::
 640
 641             sage: set_random_seed()
 642             sage: J = random_eja()
 643             sage: J.dimension() == 1 or not J.one().is_regular()
 644             True
 645
 646         """
 647         return self.degree() == self.parent().rank()
 648
 649
 650     def degree(self):
 651         """
 652         Return the degree of this element, which is defined to be
 653         the degree of its minimal polynomial.
 654
 655         ALGORITHM:
 656
 657         For now, we skip the messy minimal polynomial computation
 658         and instead return the dimension of the vector space spanned
 659         by the powers of this element. The latter is a bit more
 660         straightforward to compute.
 661
 662         SETUP::
 663
 664             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 665             ....:                                  random_eja)
 666
 667         EXAMPLES::
 668
 669             sage: J = JordanSpinEJA(4)
 670             sage: J.one().degree()
 671             1
 672             sage: e0,e1,e2,e3 = J.gens()
 673             sage: (e0 - e1).degree()
 674             2
 675
 676         In the spin factor algebra (of rank two), all elements that
 677         aren't multiples of the identity are regular::
 678
 679             sage: set_random_seed()
 680             sage: n = ZZ.random_element(1,10)
 681             sage: J = JordanSpinEJA(n)
 682             sage: x = J.random_element()
 683             sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 684             True
 685
 686         TESTS:
 687
 688         The zero and unit elements are both of degree one::
 689
 690             sage: set_random_seed()
 691             sage: J = random_eja()
 692             sage: J.zero().degree()
 693             1
 694             sage: J.one().degree()
 695             1
 696
 697         Our implementation agrees with the definition::
 698
 699             sage: set_random_seed()
 700             sage: x = random_eja().random_element()
 701             sage: x.degree() == x.minimal_polynomial().degree()
 702             True
 703
 704         """
 705         return self.span_of_powers().dimension()
 706
 707
 708     def left_matrix(self):
 709         """
 710         Our parent class defines ``left_matrix`` and ``matrix``
 711         methods whose names are misleading. We don't want them.
 712         """
 713         raise NotImplementedError("use operator().matrix() instead")
 714
 715     matrix = left_matrix
 716
 717
 718     def minimal_polynomial(self):
 719         """
 720         Return the minimal polynomial of this element,
 721         as a function of the variable `t`.
 722
 723         ALGORITHM:
 724
 725         We restrict ourselves to the associative subalgebra
 726         generated by this element, and then return the minimal
 727         polynomial of this element's operator matrix (in that
 728         subalgebra). This works by Baes Proposition 2.3.16.
 729
 730         SETUP::
 731
 732             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 733             ....:                                  random_eja)
 734
 735         TESTS:
 736
 737         The minimal polynomial of the identity and zero elements are
 738         always the same::
 739
 740             sage: set_random_seed()
 741             sage: J = random_eja()
 742             sage: J.one().minimal_polynomial()
 743             t - 1
 744             sage: J.zero().minimal_polynomial()
 745             t
 746
 747         The degree of an element is (by one definition) the degree
 748         of its minimal polynomial::
 749
 750             sage: set_random_seed()
 751             sage: x = random_eja().random_element()
 752             sage: x.degree() == x.minimal_polynomial().degree()
 753             True
 754
 755         The minimal polynomial and the characteristic polynomial coincide
 756         and are known (see Alizadeh, Example 11.11) for all elements of
 757         the spin factor algebra that aren't scalar multiples of the
 758         identity::
 759
 760             sage: set_random_seed()
 761             sage: n = ZZ.random_element(2,10)
 762             sage: J = JordanSpinEJA(n)
 763             sage: y = J.random_element()
 764             sage: while y == y.coefficient(0)*J.one():
 765             ....:     y = J.random_element()
 766             sage: y0 = y.vector()[0]
 767             sage: y_bar = y.vector()[1:]
 768             sage: actual = y.minimal_polynomial()
 769             sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
 770             sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
 771             sage: bool(actual == expected)
 772             True
 773
 774         The minimal polynomial should always kill its element::
 775
 776             sage: set_random_seed()
 777             sage: x = random_eja().random_element()
 778             sage: p = x.minimal_polynomial()
 779             sage: x.apply_univariate_polynomial(p)
 780             0
 781
 782         """
 783         V = self.span_of_powers()
 784         assoc_subalg = self.subalgebra_generated_by()
 785         # Mis-design warning: the basis used for span_of_powers()
 786         # and subalgebra_generated_by() must be the same, and in
 787         # the same order!
 788         elt = assoc_subalg(V.coordinates(self.vector()))
 789         return elt.operator().minimal_polynomial()
 790
 791
 792
 793     def natural_representation(self):
 794         """
 795         Return a more-natural representation of this element.
 796
 797         Every finite-dimensional Euclidean Jordan Algebra is a
 798         direct sum of five simple algebras, four of which comprise
 799         Hermitian matrices. This method returns the original
 800         "natural" representation of this element as a Hermitian
 801         matrix, if it has one. If not, you get the usual representation.
 802
 803         SETUP::
 804
 805             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 806             ....:                                  QuaternionHermitianEJA)
 807
 808         EXAMPLES::
 809
 810             sage: J = ComplexHermitianEJA(3)
 811             sage: J.one()
 812             e0 + e5 + e8
 813             sage: J.one().natural_representation()
 814             [1 0 0 0 0 0]
 815             [0 1 0 0 0 0]
 816             [0 0 1 0 0 0]
 817             [0 0 0 1 0 0]
 818             [0 0 0 0 1 0]
 819             [0 0 0 0 0 1]
 820
 821         ::
 822
 823             sage: J = QuaternionHermitianEJA(3)
 824             sage: J.one()
 825             e0 + e9 + e14
 826             sage: J.one().natural_representation()
 827             [1 0 0 0 0 0 0 0 0 0 0 0]
 828             [0 1 0 0 0 0 0 0 0 0 0 0]
 829             [0 0 1 0 0 0 0 0 0 0 0 0]
 830             [0 0 0 1 0 0 0 0 0 0 0 0]
 831             [0 0 0 0 1 0 0 0 0 0 0 0]
 832             [0 0 0 0 0 1 0 0 0 0 0 0]
 833             [0 0 0 0 0 0 1 0 0 0 0 0]
 834             [0 0 0 0 0 0 0 1 0 0 0 0]
 835             [0 0 0 0 0 0 0 0 1 0 0 0]
 836             [0 0 0 0 0 0 0 0 0 1 0 0]
 837             [0 0 0 0 0 0 0 0 0 0 1 0]
 838             [0 0 0 0 0 0 0 0 0 0 0 1]
 839
 840         """
 841         B = self.parent().natural_basis()
 842         W = B[0].matrix_space()
 843         return W.linear_combination(zip(self.vector(), B))
 844
 845
 846     def operator(self):
 847         """
 848         Return the left-multiplication-by-this-element
 849         operator on the ambient algebra.
 850
 851         SETUP::
 852
 853             sage: from mjo.eja.eja_algebra import random_eja
 854
 855         TESTS::
 856
 857             sage: set_random_seed()
 858             sage: J = random_eja()
 859             sage: x = J.random_element()
 860             sage: y = J.random_element()
 861             sage: x.operator()(y) == x*y
 862             True
 863             sage: y.operator()(x) == x*y
 864             True
 865
 866         """
 867         P = self.parent()
 868         fda_elt = FiniteDimensionalAlgebraElement(P, self)
 869         return FiniteDimensionalEuclideanJordanAlgebraOperator(
 870                  P,
 871                  P,
 872                  fda_elt.matrix().transpose() )
 873
 874
 875     def quadratic_representation(self, other=None):
 876         """
 877         Return the quadratic representation of this element.
 878
 879         SETUP::
 880
 881             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 882             ....:                                  random_eja)
 883
 884         EXAMPLES:
 885
 886         The explicit form in the spin factor algebra is given by
 887         Alizadeh's Example 11.12::
 888
 889             sage: set_random_seed()
 890             sage: n = ZZ.random_element(1,10)
 891             sage: J = JordanSpinEJA(n)
 892             sage: x = J.random_element()
 893             sage: x_vec = x.vector()
 894             sage: x0 = x_vec[0]
 895             sage: x_bar = x_vec[1:]
 896             sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 897             sage: B = 2*x0*x_bar.row()
 898             sage: C = 2*x0*x_bar.column()
 899             sage: D = matrix.identity(QQ, n-1)
 900             sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 901             sage: D = D + 2*x_bar.tensor_product(x_bar)
 902             sage: Q = matrix.block(2,2,[A,B,C,D])
 903             sage: Q == x.quadratic_representation().matrix()
 904             True
 905
 906         Test all of the properties from Theorem 11.2 in Alizadeh::
 907
 908             sage: set_random_seed()
 909             sage: J = random_eja()
 910             sage: x = J.random_element()
 911             sage: y = J.random_element()
 912             sage: Lx = x.operator()
 913             sage: Lxx = (x*x).operator()
 914             sage: Qx = x.quadratic_representation()
 915             sage: Qy = y.quadratic_representation()
 916             sage: Qxy = x.quadratic_representation(y)
 917             sage: Qex = J.one().quadratic_representation(x)
 918             sage: n = ZZ.random_element(10)
 919             sage: Qxn = (x^n).quadratic_representation()
 920
 921         Property 1:
 922
 923             sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
 924             True
 925
 926         Property 2 (multiply on the right for :trac:`28272`):
 927
 928             sage: alpha = QQ.random_element()
 929             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
 930             True
 931
 932         Property 3:
 933
 934             sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
 935             True
 936
 937             sage: not x.is_invertible() or (
 938             ....:   ~Qx
 939             ....:   ==
 940             ....:   x.inverse().quadratic_representation() )
 941             True
 942
 943             sage: Qxy(J.one()) == x*y
 944             True
 945
 946         Property 4:
 947
 948             sage: not x.is_invertible() or (
 949             ....:   x.quadratic_representation(x.inverse())*Qx
 950             ....:   == Qx*x.quadratic_representation(x.inverse()) )
 951             True
 952
 953             sage: not x.is_invertible() or (
 954             ....:   x.quadratic_representation(x.inverse())*Qx
 955             ....:   ==
 956             ....:   2*x.operator()*Qex - Qx )
 957             True
 958
 959             sage: 2*x.operator()*Qex - Qx == Lxx
 960             True
 961
 962         Property 5:
 963
 964             sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
 965             True
 966
 967         Property 6:
 968
 969             sage: Qxn == (Qx)^n
 970             True
 971
 972         Property 7:
 973
 974             sage: not x.is_invertible() or (
 975             ....:   Qx*x.inverse().operator() == Lx )
 976             True
 977
 978         Property 8:
 979
 980             sage: not x.operator_commutes_with(y) or (
 981             ....:   Qx(y)^n == Qxn(y^n) )
 982             True
 983
 984         """
 985         if other is None:
 986             other=self
 987         elif not other in self.parent():
 988             raise TypeError("'other' must live in the same algebra")
 989
 990         L = self.operator()
 991         M = other.operator()
 992         return ( L*M + M*L - (self*other).operator() )
 993
 994
 995     def span_of_powers(self):
 996         """
 997         Return the vector space spanned by successive powers of
 998         this element.
 999         """
1000         # The dimension of the subalgebra can't be greater than
1001         # the big algebra, so just put everything into a list
1002         # and let span() get rid of the excess.
1003         #
1004         # We do the extra ambient_vector_space() in case we're messing
1005         # with polynomials and the direct parent is a module.
1006         V = self.parent().vector_space()
1007         return V.span( (self**d).vector() for d in xrange(V.dimension()) )
1008
1009
1010     def subalgebra_generated_by(self):
1011         """
1012         Return the associative subalgebra of the parent EJA generated
1013         by this element.
1014
1015         SETUP::
1016
1017             sage: from mjo.eja.eja_algebra import random_eja
1018
1019         TESTS::
1020
1021             sage: set_random_seed()
1022             sage: x = random_eja().random_element()
1023             sage: x.subalgebra_generated_by().is_associative()
1024             True
1025
1026         Squaring in the subalgebra should work the same as in
1027         the superalgebra::
1028
1029             sage: set_random_seed()
1030             sage: x = random_eja().random_element()
1031             sage: u = x.subalgebra_generated_by().random_element()
1032             sage: u.operator()(u) == u^2
1033             True
1034
1035         """
1036         # First get the subspace spanned by the powers of myself...
1037         V = self.span_of_powers()
1038         F = self.base_ring()
1039
1040         # Now figure out the entries of the right-multiplication
1041         # matrix for the successive basis elements b0, b1,... of
1042         # that subspace.
1043         mats = []
1044         for b_right in V.basis():
1045             eja_b_right = self.parent()(b_right)
1046             b_right_rows = []
1047             # The first row of the right-multiplication matrix by
1048             # b1 is what we get if we apply that matrix to b1. The
1049             # second row of the right multiplication matrix by b1
1050             # is what we get when we apply that matrix to b2...
1051             #
1052             # IMPORTANT: this assumes that all vectors are COLUMN
1053             # vectors, unlike our superclass (which uses row vectors).
1054             for b_left in V.basis():
1055                 eja_b_left = self.parent()(b_left)
1056                 # Multiply in the original EJA, but then get the
1057                 # coordinates from the subalgebra in terms of its
1058                 # basis.
1059                 this_row = V.coordinates((eja_b_left*eja_b_right).vector())
1060                 b_right_rows.append(this_row)
1061             b_right_matrix = matrix(F, b_right_rows)
1062             mats.append(b_right_matrix)
1063
1064         # It's an algebra of polynomials in one element, and EJAs
1065         # are power-associative.
1066         #
1067         # TODO: choose generator names intelligently.
1068         #
1069         # The rank is the highest possible degree of a minimal polynomial,
1070         # and is bounded above by the dimension. We know in this case that
1071         # there's an element whose minimal polynomial has the same degree
1072         # as the space's dimension, so that must be its rank too.
1073         return FiniteDimensionalEuclideanJordanAlgebra(
1074                  F,
1075                  mats,
1076                  V.dimension(),
1077                  assume_associative=True,
1078                  names='f')
1079
1080
1081     def subalgebra_idempotent(self):
1082         """
1083         Find an idempotent in the associative subalgebra I generate
1084         using Proposition 2.3.5 in Baes.
1085
1086         SETUP::
1087
1088             sage: from mjo.eja.eja_algebra import random_eja
1089
1090         TESTS::
1091
1092             sage: set_random_seed()
1093             sage: J = random_eja()
1094             sage: x = J.random_element()
1095             sage: while x.is_nilpotent():
1096             ....:     x = J.random_element()
1097             sage: c = x.subalgebra_idempotent()
1098             sage: c^2 == c
1099             True
1100
1101         """
1102         if self.is_nilpotent():
1103             raise ValueError("this only works with non-nilpotent elements!")
1104
1105         V = self.span_of_powers()
1106         J = self.subalgebra_generated_by()
1107         # Mis-design warning: the basis used for span_of_powers()
1108         # and subalgebra_generated_by() must be the same, and in
1109         # the same order!
1110         u = J(V.coordinates(self.vector()))
1111
1112         # The image of the matrix of left-u^m-multiplication
1113         # will be minimal for some natural number s...
1114         s = 0
1115         minimal_dim = V.dimension()
1116         for i in xrange(1, V.dimension()):
1117             this_dim = (u**i).operator().matrix().image().dimension()
1118             if this_dim < minimal_dim:
1119                 minimal_dim = this_dim
1120                 s = i
1121
1122         # Now minimal_matrix should correspond to the smallest
1123         # non-zero subspace in Baes's (or really, Koecher's)
1124         # proposition.
1125         #
1126         # However, we need to restrict the matrix to work on the
1127         # subspace... or do we? Can't we just solve, knowing that
1128         # A(c) = u^(s+1) should have a solution in the big space,
1129         # too?
1130         #
1131         # Beware, solve_right() means that we're using COLUMN vectors.
1132         # Our FiniteDimensionalAlgebraElement superclass uses rows.
1133         u_next = u**(s+1)
1134         A = u_next.operator().matrix()
1135         c_coordinates = A.solve_right(u_next.vector())
1136
1137         # Now c_coordinates is the idempotent we want, but it's in
1138         # the coordinate system of the subalgebra.
1139         #
1140         # We need the basis for J, but as elements of the parent algebra.
1141         #
1142         basis = [self.parent(v) for v in V.basis()]
1143         return self.parent().linear_combination(zip(c_coordinates, basis))
1144
1145
1146     def trace(self):
1147         """
1148         Return my trace, the sum of my eigenvalues.
1149
1150         SETUP::
1151
1152             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1153             ....:                                  RealCartesianProductEJA,
1154             ....:                                  random_eja)
1155
1156         EXAMPLES::
1157
1158             sage: J = JordanSpinEJA(3)
1159             sage: x = sum(J.gens())
1160             sage: x.trace()
1161             2
1162
1163         ::
1164
1165             sage: J = RealCartesianProductEJA(5)
1166             sage: J.one().trace()
1167             5
1168
1169         TESTS:
1170
1171         The trace of an element is a real number::
1172
1173             sage: set_random_seed()
1174             sage: J = random_eja()
1175             sage: J.random_element().trace() in J.base_ring()
1176             True
1177
1178         """
1179         P = self.parent()
1180         r = P.rank()
1181         p = P._charpoly_coeff(r-1)
1182         # The _charpoly_coeff function already adds the factor of
1183         # -1 to ensure that _charpoly_coeff(r-1) is really what
1184         # appears in front of t^{r-1} in the charpoly. However,
1185         # we want the negative of THAT for the trace.
1186         return -p(*self.vector())
1187
1188
1189     def trace_inner_product(self, other):
1190         """
1191         Return the trace inner product of myself and ``other``.
1192
1193         SETUP::
1194
1195             sage: from mjo.eja.eja_algebra import random_eja
1196
1197         TESTS:
1198
1199         The trace inner product is commutative::
1200
1201             sage: set_random_seed()
1202             sage: J = random_eja()
1203             sage: x = J.random_element(); y = J.random_element()
1204             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1205             True
1206
1207         The trace inner product is bilinear::
1208
1209             sage: set_random_seed()
1210             sage: J = random_eja()
1211             sage: x = J.random_element()
1212             sage: y = J.random_element()
1213             sage: z = J.random_element()
1214             sage: a = QQ.random_element();
1215             sage: actual = (a*(x+z)).trace_inner_product(y)
1216             sage: expected = ( a*x.trace_inner_product(y) +
1217             ....:              a*z.trace_inner_product(y) )
1218             sage: actual == expected
1219             True
1220             sage: actual = x.trace_inner_product(a*(y+z))
1221             sage: expected = ( a*x.trace_inner_product(y) +
1222             ....:              a*x.trace_inner_product(z) )
1223             sage: actual == expected
1224             True
1225
1226         The trace inner product satisfies the compatibility
1227         condition in the definition of a Euclidean Jordan algebra::
1228
1229             sage: set_random_seed()
1230             sage: J = random_eja()
1231             sage: x = J.random_element()
1232             sage: y = J.random_element()
1233             sage: z = J.random_element()
1234             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1235             True
1236
1237         """
1238         if not other in self.parent():
1239             raise TypeError("'other' must live in the same algebra")
1240
1241         return (self*other).trace()