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