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         """
 464         if not self.is_invertible():
 465             raise ValueError("element is not invertible")
 466
 467         return (~self.quadratic_representation())(self)
 468
 469
 470     def is_invertible(self):
 471         """
 472         Return whether or not this element is invertible.
 473
 474         ALGORITHM:
 475
 476         The usual way to do this is to check if the determinant is
 477         zero, but we need the characteristic polynomial for the
 478         determinant. The minimal polynomial is a lot easier to get,
 479         so we use Corollary 2 in Chapter V of Koecher to check
 480         whether or not the paren't algebra's zero element is a root
 481         of this element's minimal polynomial.
 482
 483         Beware that we can't use the superclass method, because it
 484         relies on the algebra being associative.
 485
 486         SETUP::
 487
 488             sage: from mjo.eja.eja_algebra import random_eja
 489
 490         TESTS:
 491
 492         The identity element is always invertible::
 493
 494             sage: set_random_seed()
 495             sage: J = random_eja()
 496             sage: J.one().is_invertible()
 497             True
 498
 499         The zero element is never invertible in a non-trivial algebra::
 500
 501             sage: set_random_seed()
 502             sage: J = random_eja()
 503             sage: (not J.is_trivial()) and J.zero().is_invertible()
 504             False
 505
 506         """
 507         if self.is_zero():
 508             if self.parent().is_trivial():
 509                 return True
 510             else:
 511                 return False
 512
 513         # In fact, we only need to know if the constant term is non-zero,
 514         # so we can pass in the field's zero element instead.
 515         zero = self.base_ring().zero()
 516         p = self.minimal_polynomial()
 517         return not (p(zero) == zero)
 518
 519
 520     def is_nilpotent(self):
 521         """
 522         Return whether or not some power of this element is zero.
 523
 524         ALGORITHM:
 525
 526         We use Theorem 5 in Chapter III of Koecher, which says that
 527         an element ``x`` is nilpotent if and only if ``x.operator()``
 528         is nilpotent. And it is a basic fact of linear algebra that
 529         an operator on an `n`-dimensional space is nilpotent if and
 530         only if, when raised to the `n`th power, it equals the zero
 531         operator (for example, see Axler Corollary 8.8).
 532
 533         SETUP::
 534
 535             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 536             ....:                                  random_eja)
 537
 538         EXAMPLES::
 539
 540             sage: J = JordanSpinEJA(3)
 541             sage: x = sum(J.gens())
 542             sage: x.is_nilpotent()
 543             False
 544
 545         TESTS:
 546
 547         The identity element is never nilpotent::
 548
 549             sage: set_random_seed()
 550             sage: random_eja().one().is_nilpotent()
 551             False
 552
 553         The additive identity is always nilpotent::
 554
 555             sage: set_random_seed()
 556             sage: random_eja().zero().is_nilpotent()
 557             True
 558
 559         """
 560         P = self.parent()
 561         zero_operator = P.zero().operator()
 562         return self.operator()**P.dimension() == zero_operator
 563
 564
 565     def is_regular(self):
 566         """
 567         Return whether or not this is a regular element.
 568
 569         SETUP::
 570
 571             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 572             ....:                                  random_eja)
 573
 574         EXAMPLES:
 575
 576         The identity element always has degree one, but any element
 577         linearly-independent from it is regular::
 578
 579             sage: J = JordanSpinEJA(5)
 580             sage: J.one().is_regular()
 581             False
 582             sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 583             sage: for x in J.gens():
 584             ....:     (J.one() + x).is_regular()
 585             False
 586             True
 587             True
 588             True
 589             True
 590
 591         TESTS:
 592
 593         The zero element should never be regular, unless the parent
 594         algebra has dimension one::
 595
 596             sage: set_random_seed()
 597             sage: J = random_eja()
 598             sage: J.dimension() == 1 or not J.zero().is_regular()
 599             True
 600
 601         The unit element isn't regular unless the algebra happens to
 602         consist of only its scalar multiples::
 603
 604             sage: set_random_seed()
 605             sage: J = random_eja()
 606             sage: J.dimension() == 1 or not J.one().is_regular()
 607             True
 608
 609         """
 610         return self.degree() == self.parent().rank()
 611
 612
 613     def degree(self):
 614         """
 615         Return the degree of this element, which is defined to be
 616         the degree of its minimal polynomial.
 617
 618         ALGORITHM:
 619
 620         For now, we skip the messy minimal polynomial computation
 621         and instead return the dimension of the vector space spanned
 622         by the powers of this element. The latter is a bit more
 623         straightforward to compute.
 624
 625         SETUP::
 626
 627             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 628             ....:                                  random_eja)
 629
 630         EXAMPLES::
 631
 632             sage: J = JordanSpinEJA(4)
 633             sage: J.one().degree()
 634             1
 635             sage: e0,e1,e2,e3 = J.gens()
 636             sage: (e0 - e1).degree()
 637             2
 638
 639         In the spin factor algebra (of rank two), all elements that
 640         aren't multiples of the identity are regular::
 641
 642             sage: set_random_seed()
 643             sage: J = JordanSpinEJA.random_instance()
 644             sage: x = J.random_element()
 645             sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 646             True
 647
 648         TESTS:
 649
 650         The zero and unit elements are both of degree one::
 651
 652             sage: set_random_seed()
 653             sage: J = random_eja()
 654             sage: J.zero().degree()
 655             1
 656             sage: J.one().degree()
 657             1
 658
 659         Our implementation agrees with the definition::
 660
 661             sage: set_random_seed()
 662             sage: x = random_eja().random_element()
 663             sage: x.degree() == x.minimal_polynomial().degree()
 664             True
 665
 666         """
 667         if self.is_zero() and not self.parent().is_trivial():
 668             # The minimal polynomial of zero in a nontrivial algebra
 669             # is "t"; in a trivial algebra it's "1" by convention
 670             # (it's an empty product).
 671             return 1
 672         return self.subalgebra_generated_by().dimension()
 673
 674
 675     def left_matrix(self):
 676         """
 677         Our parent class defines ``left_matrix`` and ``matrix``
 678         methods whose names are misleading. We don't want them.
 679         """
 680         raise NotImplementedError("use operator().matrix() instead")
 681
 682     matrix = left_matrix
 683
 684
 685     def minimal_polynomial(self):
 686         """
 687         Return the minimal polynomial of this element,
 688         as a function of the variable `t`.
 689
 690         ALGORITHM:
 691
 692         We restrict ourselves to the associative subalgebra
 693         generated by this element, and then return the minimal
 694         polynomial of this element's operator matrix (in that
 695         subalgebra). This works by Baes Proposition 2.3.16.
 696
 697         SETUP::
 698
 699             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 700             ....:                                  RealSymmetricEJA,
 701             ....:                                  random_eja)
 702
 703         TESTS:
 704
 705         The minimal polynomial of the identity and zero elements are
 706         always the same::
 707
 708             sage: set_random_seed()
 709             sage: J = random_eja()
 710             sage: J.one().minimal_polynomial()
 711             t - 1
 712             sage: J.zero().minimal_polynomial()
 713             t
 714
 715         The degree of an element is (by one definition) the degree
 716         of its minimal polynomial::
 717
 718             sage: set_random_seed()
 719             sage: x = random_eja().random_element()
 720             sage: x.degree() == x.minimal_polynomial().degree()
 721             True
 722
 723         The minimal polynomial and the characteristic polynomial coincide
 724         and are known (see Alizadeh, Example 11.11) for all elements of
 725         the spin factor algebra that aren't scalar multiples of the
 726         identity. We require the dimension of the algebra to be at least
 727         two here so that said elements actually exist::
 728
 729             sage: set_random_seed()
 730             sage: n_max = max(2, JordanSpinEJA._max_test_case_size())
 731             sage: n = ZZ.random_element(2, n_max)
 732             sage: J = JordanSpinEJA(n)
 733             sage: y = J.random_element()
 734             sage: while y == y.coefficient(0)*J.one():
 735             ....:     y = J.random_element()
 736             sage: y0 = y.to_vector()[0]
 737             sage: y_bar = y.to_vector()[1:]
 738             sage: actual = y.minimal_polynomial()
 739             sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
 740             sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
 741             sage: bool(actual == expected)
 742             True
 743
 744         The minimal polynomial should always kill its element::
 745
 746             sage: set_random_seed()
 747             sage: x = random_eja().random_element()
 748             sage: p = x.minimal_polynomial()
 749             sage: x.apply_univariate_polynomial(p)
 750             0
 751
 752         The minimal polynomial is invariant under a change of basis,
 753         and in particular, a re-scaling of the basis::
 754
 755             sage: set_random_seed()
 756             sage: n_max = RealSymmetricEJA._max_test_case_size()
 757             sage: n = ZZ.random_element(1, n_max)
 758             sage: J1 = RealSymmetricEJA(n,QQ)
 759             sage: J2 = RealSymmetricEJA(n,QQ,normalize_basis=False)
 760             sage: X = random_matrix(QQ,n)
 761             sage: X = X*X.transpose()
 762             sage: x1 = J1(X)
 763             sage: x2 = J2(X)
 764             sage: x1.minimal_polynomial() == x2.minimal_polynomial()
 765             True
 766
 767         """
 768         if self.is_zero():
 769             # We would generate a zero-dimensional subalgebra
 770             # where the minimal polynomial would be constant.
 771             # That might be correct, but only if *this* algebra
 772             # is trivial too.
 773             if not self.parent().is_trivial():
 774                 # Pretty sure we know what the minimal polynomial of
 775                 # the zero operator is going to be. This ensures
 776                 # consistency of e.g. the polynomial variable returned
 777                 # in the "normal" case without us having to think about it.
 778                 return self.operator().minimal_polynomial()
 779
 780         A = self.subalgebra_generated_by()
 781         return A(self).operator().minimal_polynomial()
 782
 783
 784
 785     def natural_representation(self):
 786         """
 787         Return a more-natural representation of this element.
 788
 789         Every finite-dimensional Euclidean Jordan Algebra is a
 790         direct sum of five simple algebras, four of which comprise
 791         Hermitian matrices. This method returns the original
 792         "natural" representation of this element as a Hermitian
 793         matrix, if it has one. If not, you get the usual representation.
 794
 795         SETUP::
 796
 797             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 798             ....:                                  QuaternionHermitianEJA)
 799
 800         EXAMPLES::
 801
 802             sage: J = ComplexHermitianEJA(3)
 803             sage: J.one()
 804             e0 + e3 + e8
 805             sage: J.one().natural_representation()
 806             [1 0 0 0 0 0]
 807             [0 1 0 0 0 0]
 808             [0 0 1 0 0 0]
 809             [0 0 0 1 0 0]
 810             [0 0 0 0 1 0]
 811             [0 0 0 0 0 1]
 812
 813         ::
 814
 815             sage: J = QuaternionHermitianEJA(3)
 816             sage: J.one()
 817             e0 + e5 + e14
 818             sage: J.one().natural_representation()
 819             [1 0 0 0 0 0 0 0 0 0 0 0]
 820             [0 1 0 0 0 0 0 0 0 0 0 0]
 821             [0 0 1 0 0 0 0 0 0 0 0 0]
 822             [0 0 0 1 0 0 0 0 0 0 0 0]
 823             [0 0 0 0 1 0 0 0 0 0 0 0]
 824             [0 0 0 0 0 1 0 0 0 0 0 0]
 825             [0 0 0 0 0 0 1 0 0 0 0 0]
 826             [0 0 0 0 0 0 0 1 0 0 0 0]
 827             [0 0 0 0 0 0 0 0 1 0 0 0]
 828             [0 0 0 0 0 0 0 0 0 1 0 0]
 829             [0 0 0 0 0 0 0 0 0 0 1 0]
 830             [0 0 0 0 0 0 0 0 0 0 0 1]
 831
 832         """
 833         B = self.parent().natural_basis()
 834         W = self.parent().natural_basis_space()
 835         return W.linear_combination(izip(B,self.to_vector()))
 836
 837
 838     def norm(self):
 839         """
 840         The norm of this element with respect to :meth:`inner_product`.
 841
 842         SETUP::
 843
 844             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 845             ....:                                  RealCartesianProductEJA)
 846
 847         EXAMPLES::
 848
 849             sage: J = RealCartesianProductEJA(2)
 850             sage: x = sum(J.gens())
 851             sage: x.norm()
 852             sqrt(2)
 853
 854         ::
 855
 856             sage: J = JordanSpinEJA(4)
 857             sage: x = sum(J.gens())
 858             sage: x.norm()
 859             2
 860
 861         """
 862         return self.inner_product(self).sqrt()
 863
 864
 865     def operator(self):
 866         """
 867         Return the left-multiplication-by-this-element
 868         operator on the ambient algebra.
 869
 870         SETUP::
 871
 872             sage: from mjo.eja.eja_algebra import random_eja
 873
 874         TESTS::
 875
 876             sage: set_random_seed()
 877             sage: J = random_eja()
 878             sage: x,y = J.random_elements(2)
 879             sage: x.operator()(y) == x*y
 880             True
 881             sage: y.operator()(x) == x*y
 882             True
 883
 884         """
 885         P = self.parent()
 886         left_mult_by_self = lambda y: self*y
 887         L = P.module_morphism(function=left_mult_by_self, codomain=P)
 888         return FiniteDimensionalEuclideanJordanAlgebraOperator(
 889                  P,
 890                  P,
 891                  L.matrix() )
 892
 893
 894     def quadratic_representation(self, other=None):
 895         """
 896         Return the quadratic representation of this element.
 897
 898         SETUP::
 899
 900             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 901             ....:                                  random_eja)
 902
 903         EXAMPLES:
 904
 905         The explicit form in the spin factor algebra is given by
 906         Alizadeh's Example 11.12::
 907
 908             sage: set_random_seed()
 909             sage: x = JordanSpinEJA.random_instance().random_element()
 910             sage: x_vec = x.to_vector()
 911             sage: n = x_vec.degree()
 912             sage: x0 = x_vec[0]
 913             sage: x_bar = x_vec[1:]
 914             sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 915             sage: B = 2*x0*x_bar.row()
 916             sage: C = 2*x0*x_bar.column()
 917             sage: D = matrix.identity(QQ, n-1)
 918             sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 919             sage: D = D + 2*x_bar.tensor_product(x_bar)
 920             sage: Q = matrix.block(2,2,[A,B,C,D])
 921             sage: Q == x.quadratic_representation().matrix()
 922             True
 923
 924         Test all of the properties from Theorem 11.2 in Alizadeh::
 925
 926             sage: set_random_seed()
 927             sage: J = random_eja()
 928             sage: x,y = J.random_elements(2)
 929             sage: Lx = x.operator()
 930             sage: Lxx = (x*x).operator()
 931             sage: Qx = x.quadratic_representation()
 932             sage: Qy = y.quadratic_representation()
 933             sage: Qxy = x.quadratic_representation(y)
 934             sage: Qex = J.one().quadratic_representation(x)
 935             sage: n = ZZ.random_element(10)
 936             sage: Qxn = (x^n).quadratic_representation()
 937
 938         Property 1:
 939
 940             sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
 941             True
 942
 943         Property 2 (multiply on the right for :trac:`28272`):
 944
 945             sage: alpha = J.base_ring().random_element()
 946             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
 947             True
 948
 949         Property 3:
 950
 951             sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
 952             True
 953
 954             sage: not x.is_invertible() or (
 955             ....:   ~Qx
 956             ....:   ==
 957             ....:   x.inverse().quadratic_representation() )
 958             True
 959
 960             sage: Qxy(J.one()) == x*y
 961             True
 962
 963         Property 4:
 964
 965             sage: not x.is_invertible() or (
 966             ....:   x.quadratic_representation(x.inverse())*Qx
 967             ....:   == Qx*x.quadratic_representation(x.inverse()) )
 968             True
 969
 970             sage: not x.is_invertible() or (
 971             ....:   x.quadratic_representation(x.inverse())*Qx
 972             ....:   ==
 973             ....:   2*Lx*Qex - Qx )
 974             True
 975
 976             sage: 2*Lx*Qex - Qx == Lxx
 977             True
 978
 979         Property 5:
 980
 981             sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
 982             True
 983
 984         Property 6:
 985
 986             sage: Qxn == (Qx)^n
 987             True
 988
 989         Property 7:
 990
 991             sage: not x.is_invertible() or (
 992             ....:   Qx*x.inverse().operator() == Lx )
 993             True
 994
 995         Property 8:
 996
 997             sage: not x.operator_commutes_with(y) or (
 998             ....:   Qx(y)^n == Qxn(y^n) )
 999             True
1000
1001         """
1002         if other is None:
1003             other=self
1004         elif not other in self.parent():
1005             raise TypeError("'other' must live in the same algebra")
1006
1007         L = self.operator()
1008         M = other.operator()
1009         return ( L*M + M*L - (self*other).operator() )
1010
1011
1012
1013
1014     def subalgebra_generated_by(self, orthonormalize_basis=False):
1015         """
1016         Return the associative subalgebra of the parent EJA generated
1017         by this element.
1018
1019         SETUP::
1020
1021             sage: from mjo.eja.eja_algebra import random_eja
1022
1023         TESTS:
1024
1025         This subalgebra, being composed of only powers, is associative::
1026
1027             sage: set_random_seed()
1028             sage: x0 = random_eja().random_element()
1029             sage: A = x0.subalgebra_generated_by()
1030             sage: x,y,z = A.random_elements(3)
1031             sage: (x*y)*z == x*(y*z)
1032             True
1033
1034         Squaring in the subalgebra should work the same as in
1035         the superalgebra::
1036
1037             sage: set_random_seed()
1038             sage: x = random_eja().random_element()
1039             sage: A = x.subalgebra_generated_by()
1040             sage: A(x^2) == A(x)*A(x)
1041             True
1042
1043         The subalgebra generated by the zero element is trivial::
1044
1045             sage: set_random_seed()
1046             sage: A = random_eja().zero().subalgebra_generated_by()
1047             sage: A
1048             Euclidean Jordan algebra of dimension 0 over...
1049             sage: A.one()
1050             0
1051
1052         """
1053         return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
1054
1055
1056     def subalgebra_idempotent(self):
1057         """
1058         Find an idempotent in the associative subalgebra I generate
1059         using Proposition 2.3.5 in Baes.
1060
1061         SETUP::
1062
1063             sage: from mjo.eja.eja_algebra import random_eja
1064
1065         TESTS::
1066
1067             sage: set_random_seed()
1068             sage: J = random_eja()
1069             sage: x = J.random_element()
1070             sage: while x.is_nilpotent():
1071             ....:     x = J.random_element()
1072             sage: c = x.subalgebra_idempotent()
1073             sage: c^2 == c
1074             True
1075
1076         """
1077         if self.is_nilpotent():
1078             raise ValueError("this only works with non-nilpotent elements!")
1079
1080         J = self.subalgebra_generated_by()
1081         u = J(self)
1082
1083         # The image of the matrix of left-u^m-multiplication
1084         # will be minimal for some natural number s...
1085         s = 0
1086         minimal_dim = J.dimension()
1087         for i in xrange(1, minimal_dim):
1088             this_dim = (u**i).operator().matrix().image().dimension()
1089             if this_dim < minimal_dim:
1090                 minimal_dim = this_dim
1091                 s = i
1092
1093         # Now minimal_matrix should correspond to the smallest
1094         # non-zero subspace in Baes's (or really, Koecher's)
1095         # proposition.
1096         #
1097         # However, we need to restrict the matrix to work on the
1098         # subspace... or do we? Can't we just solve, knowing that
1099         # A(c) = u^(s+1) should have a solution in the big space,
1100         # too?
1101         #
1102         # Beware, solve_right() means that we're using COLUMN vectors.
1103         # Our FiniteDimensionalAlgebraElement superclass uses rows.
1104         u_next = u**(s+1)
1105         A = u_next.operator().matrix()
1106         c = J.from_vector(A.solve_right(u_next.to_vector()))
1107
1108         # Now c is the idempotent we want, but it still lives in the subalgebra.
1109         return c.superalgebra_element()
1110
1111
1112     def trace(self):
1113         """
1114         Return my trace, the sum of my eigenvalues.
1115
1116         SETUP::
1117
1118             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1119             ....:                                  RealCartesianProductEJA,
1120             ....:                                  random_eja)
1121
1122         EXAMPLES::
1123
1124             sage: J = JordanSpinEJA(3)
1125             sage: x = sum(J.gens())
1126             sage: x.trace()
1127             2
1128
1129         ::
1130
1131             sage: J = RealCartesianProductEJA(5)
1132             sage: J.one().trace()
1133             5
1134
1135         TESTS:
1136
1137         The trace of an element is a real number::
1138
1139             sage: set_random_seed()
1140             sage: J = random_eja()
1141             sage: J.random_element().trace() in RLF
1142             True
1143
1144         """
1145         P = self.parent()
1146         r = P.rank()
1147         p = P._charpoly_coeff(r-1)
1148         # The _charpoly_coeff function already adds the factor of
1149         # -1 to ensure that _charpoly_coeff(r-1) is really what
1150         # appears in front of t^{r-1} in the charpoly. However,
1151         # we want the negative of THAT for the trace.
1152         return -p(*self.to_vector())
1153
1154
1155     def trace_inner_product(self, other):
1156         """
1157         Return the trace inner product of myself and ``other``.
1158
1159         SETUP::
1160
1161             sage: from mjo.eja.eja_algebra import random_eja
1162
1163         TESTS:
1164
1165         The trace inner product is commutative, bilinear, and associative::
1166
1167             sage: set_random_seed()
1168             sage: J = random_eja()
1169             sage: x,y,z = J.random_elements(3)
1170             sage: # commutative
1171             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1172             True
1173             sage: # bilinear
1174             sage: a = J.base_ring().random_element();
1175             sage: actual = (a*(x+z)).trace_inner_product(y)
1176             sage: expected = ( a*x.trace_inner_product(y) +
1177             ....:              a*z.trace_inner_product(y) )
1178             sage: actual == expected
1179             True
1180             sage: actual = x.trace_inner_product(a*(y+z))
1181             sage: expected = ( a*x.trace_inner_product(y) +
1182             ....:              a*x.trace_inner_product(z) )
1183             sage: actual == expected
1184             True
1185             sage: # associative
1186             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1187             True
1188
1189         """
1190         if not other in self.parent():
1191             raise TypeError("'other' must live in the same algebra")
1192
1193         return (self*other).trace()
1194
1195
1196     def trace_norm(self):
1197         """
1198         The norm of this element with respect to :meth:`trace_inner_product`.
1199
1200         SETUP::
1201
1202             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1203             ....:                                  RealCartesianProductEJA)
1204
1205         EXAMPLES::
1206
1207             sage: J = RealCartesianProductEJA(2)
1208             sage: x = sum(J.gens())
1209             sage: x.trace_norm()
1210             sqrt(2)
1211
1212         ::
1213
1214             sage: J = JordanSpinEJA(4)
1215             sage: x = sum(J.gens())
1216             sage: x.trace_norm()
1217             2*sqrt(2)
1218
1219         """
1220         return self.trace_inner_product(self).sqrt()