mjo/eja/eja_element.py

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