mjo/eja/eja_element.py

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