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         The characteristic polynomials of the zero and unit elements
 169         should be what we think they are in a subalgebra, too::
 170
 171             sage: J = RealCartesianProductEJA(3)
 172             sage: p1 = J.one().characteristic_polynomial()
 173             sage: q1 = J.zero().characteristic_polynomial()
 174             sage: e0,e1,e2 = J.gens()
 175             sage: A = (e0 + 2*e1 + 3*e2).subalgebra_generated_by() # dim 3
 176             sage: p2 = A.one().characteristic_polynomial()
 177             sage: q2 = A.zero().characteristic_polynomial()
 178             sage: p1 == p2
 179             True
 180             sage: q1 == q2
 181             True
 182
 183         """
 184         p = self.parent().characteristic_polynomial()
 185         return p(*self.to_vector())
 186
 187
 188     def inner_product(self, other):
 189         """
 190         Return the parent algebra's inner product of myself and ``other``.
 191
 192         SETUP::
 193
 194             sage: from mjo.eja.eja_algebra import (
 195             ....:   ComplexHermitianEJA,
 196             ....:   JordanSpinEJA,
 197             ....:   QuaternionHermitianEJA,
 198             ....:   RealSymmetricEJA,
 199             ....:   random_eja)
 200
 201         EXAMPLES:
 202
 203         The inner product in the Jordan spin algebra is the usual
 204         inner product on `R^n` (this example only works because the
 205         basis for the Jordan algebra is the standard basis in `R^n`)::
 206
 207             sage: J = JordanSpinEJA(3)
 208             sage: x = vector(QQ,[1,2,3])
 209             sage: y = vector(QQ,[4,5,6])
 210             sage: x.inner_product(y)
 211             32
 212             sage: J.from_vector(x).inner_product(J.from_vector(y))
 213             32
 214
 215         The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
 216         multiplication is the usual matrix multiplication in `S^n`,
 217         so the inner product of the identity matrix with itself
 218         should be the `n`::
 219
 220             sage: J = RealSymmetricEJA(3)
 221             sage: J.one().inner_product(J.one())
 222             3
 223
 224         Likewise, the inner product on `C^n` is `<X,Y> =
 225         Re(trace(X*Y))`, where we must necessarily take the real
 226         part because the product of Hermitian matrices may not be
 227         Hermitian::
 228
 229             sage: J = ComplexHermitianEJA(3)
 230             sage: J.one().inner_product(J.one())
 231             3
 232
 233         Ditto for the quaternions::
 234
 235             sage: J = QuaternionHermitianEJA(3)
 236             sage: J.one().inner_product(J.one())
 237             3
 238
 239         TESTS:
 240
 241         Ensure that we can always compute an inner product, and that
 242         it gives us back a real number::
 243
 244             sage: set_random_seed()
 245             sage: J = random_eja()
 246             sage: x = J.random_element()
 247             sage: y = J.random_element()
 248             sage: x.inner_product(y) in RR
 249             True
 250
 251         """
 252         P = self.parent()
 253         if not other in P:
 254             raise TypeError("'other' must live in the same algebra")
 255
 256         return P.inner_product(self, other)
 257
 258
 259     def operator_commutes_with(self, other):
 260         """
 261         Return whether or not this element operator-commutes
 262         with ``other``.
 263
 264         SETUP::
 265
 266             sage: from mjo.eja.eja_algebra import random_eja
 267
 268         EXAMPLES:
 269
 270         The definition of a Jordan algebra says that any element
 271         operator-commutes with its square::
 272
 273             sage: set_random_seed()
 274             sage: x = random_eja().random_element()
 275             sage: x.operator_commutes_with(x^2)
 276             True
 277
 278         TESTS:
 279
 280         Test Lemma 1 from Chapter III of Koecher::
 281
 282             sage: set_random_seed()
 283             sage: J = random_eja()
 284             sage: u = J.random_element()
 285             sage: v = J.random_element()
 286             sage: lhs = u.operator_commutes_with(u*v)
 287             sage: rhs = v.operator_commutes_with(u^2)
 288             sage: lhs == rhs
 289             True
 290
 291         Test the first polarization identity from my notes, Koecher
 292         Chapter III, or from Baes (2.3)::
 293
 294             sage: set_random_seed()
 295             sage: J = random_eja()
 296             sage: x = J.random_element()
 297             sage: y = J.random_element()
 298             sage: Lx = x.operator()
 299             sage: Ly = y.operator()
 300             sage: Lxx = (x*x).operator()
 301             sage: Lxy = (x*y).operator()
 302             sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
 303             True
 304
 305         Test the second polarization identity from my notes or from
 306         Baes (2.4)::
 307
 308             sage: set_random_seed()
 309             sage: J = random_eja()
 310             sage: x = J.random_element()
 311             sage: y = J.random_element()
 312             sage: z = J.random_element()
 313             sage: Lx = x.operator()
 314             sage: Ly = y.operator()
 315             sage: Lz = z.operator()
 316             sage: Lzy = (z*y).operator()
 317             sage: Lxy = (x*y).operator()
 318             sage: Lxz = (x*z).operator()
 319             sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
 320             True
 321
 322         Test the third polarization identity from my notes or from
 323         Baes (2.5)::
 324
 325             sage: set_random_seed()
 326             sage: J = random_eja()
 327             sage: u = J.random_element()
 328             sage: y = J.random_element()
 329             sage: z = J.random_element()
 330             sage: Lu = u.operator()
 331             sage: Ly = y.operator()
 332             sage: Lz = z.operator()
 333             sage: Lzy = (z*y).operator()
 334             sage: Luy = (u*y).operator()
 335             sage: Luz = (u*z).operator()
 336             sage: Luyz = (u*(y*z)).operator()
 337             sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
 338             sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
 339             sage: bool(lhs == rhs)
 340             True
 341
 342         """
 343         if not other in self.parent():
 344             raise TypeError("'other' must live in the same algebra")
 345
 346         A = self.operator()
 347         B = other.operator()
 348         return (A*B == B*A)
 349
 350
 351     def det(self):
 352         """
 353         Return my determinant, the product of my eigenvalues.
 354
 355         SETUP::
 356
 357             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 358             ....:                                  random_eja)
 359
 360         EXAMPLES::
 361
 362             sage: J = JordanSpinEJA(2)
 363             sage: e0,e1 = J.gens()
 364             sage: x = sum( J.gens() )
 365             sage: x.det()
 366             0
 367
 368         ::
 369
 370             sage: J = JordanSpinEJA(3)
 371             sage: e0,e1,e2 = J.gens()
 372             sage: x = sum( J.gens() )
 373             sage: x.det()
 374             -1
 375
 376         TESTS:
 377
 378         An element is invertible if and only if its determinant is
 379         non-zero::
 380
 381             sage: set_random_seed()
 382             sage: x = random_eja().random_element()
 383             sage: x.is_invertible() == (x.det() != 0)
 384             True
 385
 386         Ensure that the determinant is multiplicative on an associative
 387         subalgebra as in Faraut and Koranyi's Proposition II.2.2::
 388
 389             sage: set_random_seed()
 390             sage: J = random_eja().random_element().subalgebra_generated_by()
 391             sage: x = J.random_element()
 392             sage: y = J.random_element()
 393             sage: (x*y).det() == x.det()*y.det()
 394             True
 395
 396         """
 397         P = self.parent()
 398         r = P.rank()
 399         p = P._charpoly_coeff(0)
 400         # The _charpoly_coeff function already adds the factor of
 401         # -1 to ensure that _charpoly_coeff(0) is really what
 402         # appears in front of t^{0} in the charpoly. However,
 403         # we want (-1)^r times THAT for the determinant.
 404         return ((-1)**r)*p(*self.to_vector())
 405
 406
 407     def inverse(self):
 408         """
 409         Return the Jordan-multiplicative inverse of this element.
 410
 411         ALGORITHM:
 412
 413         We appeal to the quadratic representation as in Koecher's
 414         Theorem 12 in Chapter III, Section 5.
 415
 416         SETUP::
 417
 418             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 419             ....:                                  random_eja)
 420
 421         EXAMPLES:
 422
 423         The inverse in the spin factor algebra is given in Alizadeh's
 424         Example 11.11::
 425
 426             sage: set_random_seed()
 427             sage: n = ZZ.random_element(1,10)
 428             sage: J = JordanSpinEJA(n)
 429             sage: x = J.random_element()
 430             sage: while not x.is_invertible():
 431             ....:     x = J.random_element()
 432             sage: x_vec = x.to_vector()
 433             sage: x0 = x_vec[0]
 434             sage: x_bar = x_vec[1:]
 435             sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
 436             sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
 437             sage: x_inverse = coeff*inv_vec
 438             sage: x.inverse() == J.from_vector(x_inverse)
 439             True
 440
 441         TESTS:
 442
 443         The identity element is its own inverse::
 444
 445             sage: set_random_seed()
 446             sage: J = random_eja()
 447             sage: J.one().inverse() == J.one()
 448             True
 449
 450         If an element has an inverse, it acts like one::
 451
 452             sage: set_random_seed()
 453             sage: J = random_eja()
 454             sage: x = J.random_element()
 455             sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
 456             True
 457
 458         The inverse of the inverse is what we started with::
 459
 460             sage: set_random_seed()
 461             sage: J = random_eja()
 462             sage: x = J.random_element()
 463             sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
 464             True
 465
 466         The zero element is never invertible::
 467
 468             sage: set_random_seed()
 469             sage: J = random_eja().zero().inverse()
 470             Traceback (most recent call last):
 471             ...
 472             ValueError: element is not invertible
 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: n = ZZ.random_element(1,10)
 655             sage: J = JordanSpinEJA(n)
 656             sage: x = J.random_element()
 657             sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 658             True
 659
 660         TESTS:
 661
 662         The zero and unit elements are both of degree one::
 663
 664             sage: set_random_seed()
 665             sage: J = random_eja()
 666             sage: J.zero().degree()
 667             1
 668             sage: J.one().degree()
 669             1
 670
 671         Our implementation agrees with the definition::
 672
 673             sage: set_random_seed()
 674             sage: x = random_eja().random_element()
 675             sage: x.degree() == x.minimal_polynomial().degree()
 676             True
 677
 678         """
 679         if self.is_zero() and not self.parent().is_trivial():
 680             # The minimal polynomial of zero in a nontrivial algebra
 681             # is "t"; in a trivial algebra it's "1" by convention
 682             # (it's an empty product).
 683             return 1
 684         return self.subalgebra_generated_by().dimension()
 685
 686
 687     def left_matrix(self):
 688         """
 689         Our parent class defines ``left_matrix`` and ``matrix``
 690         methods whose names are misleading. We don't want them.
 691         """
 692         raise NotImplementedError("use operator().matrix() instead")
 693
 694     matrix = left_matrix
 695
 696
 697     def minimal_polynomial(self):
 698         """
 699         Return the minimal polynomial of this element,
 700         as a function of the variable `t`.
 701
 702         ALGORITHM:
 703
 704         We restrict ourselves to the associative subalgebra
 705         generated by this element, and then return the minimal
 706         polynomial of this element's operator matrix (in that
 707         subalgebra). This works by Baes Proposition 2.3.16.
 708
 709         SETUP::
 710
 711             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 712             ....:                                  RealSymmetricEJA,
 713             ....:                                  random_eja)
 714
 715         TESTS:
 716
 717         The minimal polynomial of the identity and zero elements are
 718         always the same::
 719
 720             sage: set_random_seed()
 721             sage: J = random_eja()
 722             sage: J.one().minimal_polynomial()
 723             t - 1
 724             sage: J.zero().minimal_polynomial()
 725             t
 726
 727         The degree of an element is (by one definition) the degree
 728         of its minimal polynomial::
 729
 730             sage: set_random_seed()
 731             sage: x = random_eja().random_element()
 732             sage: x.degree() == x.minimal_polynomial().degree()
 733             True
 734
 735         The minimal polynomial and the characteristic polynomial coincide
 736         and are known (see Alizadeh, Example 11.11) for all elements of
 737         the spin factor algebra that aren't scalar multiples of the
 738         identity::
 739
 740             sage: set_random_seed()
 741             sage: n = ZZ.random_element(2,10)
 742             sage: J = JordanSpinEJA(n)
 743             sage: y = J.random_element()
 744             sage: while y == y.coefficient(0)*J.one():
 745             ....:     y = J.random_element()
 746             sage: y0 = y.to_vector()[0]
 747             sage: y_bar = y.to_vector()[1:]
 748             sage: actual = y.minimal_polynomial()
 749             sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
 750             sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
 751             sage: bool(actual == expected)
 752             True
 753
 754         The minimal polynomial should always kill its element::
 755
 756             sage: set_random_seed()
 757             sage: x = random_eja().random_element()
 758             sage: p = x.minimal_polynomial()
 759             sage: x.apply_univariate_polynomial(p)
 760             0
 761
 762         The minimal polynomial is invariant under a change of basis,
 763         and in particular, a re-scaling of the basis::
 764
 765             sage: set_random_seed()
 766             sage: n = ZZ.random_element(1,5)
 767             sage: J1 = RealSymmetricEJA(n)
 768             sage: J2 = RealSymmetricEJA(n,QQ,False)
 769             sage: X = random_matrix(QQ,n)
 770             sage: X = X*X.transpose()
 771             sage: x1 = J1(X)
 772             sage: x2 = J2(X)
 773             sage: x1.minimal_polynomial() == x2.minimal_polynomial()
 774             True
 775
 776         """
 777         if self.is_zero():
 778             # We would generate a zero-dimensional subalgebra
 779             # where the minimal polynomial would be constant.
 780             # That might be correct, but only if *this* algebra
 781             # is trivial too.
 782             if not self.parent().is_trivial():
 783                 # Pretty sure we know what the minimal polynomial of
 784                 # the zero operator is going to be. This ensures
 785                 # consistency of e.g. the polynomial variable returned
 786                 # in the "normal" case without us having to think about it.
 787                 return self.operator().minimal_polynomial()
 788
 789         A = self.subalgebra_generated_by()
 790         return A(self).operator().minimal_polynomial()
 791
 792
 793
 794     def natural_representation(self):
 795         """
 796         Return a more-natural representation of this element.
 797
 798         Every finite-dimensional Euclidean Jordan Algebra is a
 799         direct sum of five simple algebras, four of which comprise
 800         Hermitian matrices. This method returns the original
 801         "natural" representation of this element as a Hermitian
 802         matrix, if it has one. If not, you get the usual representation.
 803
 804         SETUP::
 805
 806             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 807             ....:                                  QuaternionHermitianEJA)
 808
 809         EXAMPLES::
 810
 811             sage: J = ComplexHermitianEJA(3)
 812             sage: J.one()
 813             e0 + e3 + e8
 814             sage: J.one().natural_representation()
 815             [1 0 0 0 0 0]
 816             [0 1 0 0 0 0]
 817             [0 0 1 0 0 0]
 818             [0 0 0 1 0 0]
 819             [0 0 0 0 1 0]
 820             [0 0 0 0 0 1]
 821
 822         ::
 823
 824             sage: J = QuaternionHermitianEJA(3)
 825             sage: J.one()
 826             e0 + e5 + e14
 827             sage: J.one().natural_representation()
 828             [1 0 0 0 0 0 0 0 0 0 0 0]
 829             [0 1 0 0 0 0 0 0 0 0 0 0]
 830             [0 0 1 0 0 0 0 0 0 0 0 0]
 831             [0 0 0 1 0 0 0 0 0 0 0 0]
 832             [0 0 0 0 1 0 0 0 0 0 0 0]
 833             [0 0 0 0 0 1 0 0 0 0 0 0]
 834             [0 0 0 0 0 0 1 0 0 0 0 0]
 835             [0 0 0 0 0 0 0 1 0 0 0 0]
 836             [0 0 0 0 0 0 0 0 1 0 0 0]
 837             [0 0 0 0 0 0 0 0 0 1 0 0]
 838             [0 0 0 0 0 0 0 0 0 0 1 0]
 839             [0 0 0 0 0 0 0 0 0 0 0 1]
 840
 841         """
 842         B = self.parent().natural_basis()
 843         W = self.parent().natural_basis_space()
 844         return W.linear_combination(zip(B,self.to_vector()))
 845
 846
 847     def norm(self):
 848         """
 849         The norm of this element with respect to :meth:`inner_product`.
 850
 851         SETUP::
 852
 853             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 854             ....:                                  RealCartesianProductEJA)
 855
 856         EXAMPLES::
 857
 858             sage: J = RealCartesianProductEJA(2)
 859             sage: x = sum(J.gens())
 860             sage: x.norm()
 861             sqrt(2)
 862
 863         ::
 864
 865             sage: J = JordanSpinEJA(4)
 866             sage: x = sum(J.gens())
 867             sage: x.norm()
 868             2
 869
 870         """
 871         return self.inner_product(self).sqrt()
 872
 873
 874     def operator(self):
 875         """
 876         Return the left-multiplication-by-this-element
 877         operator on the ambient algebra.
 878
 879         SETUP::
 880
 881             sage: from mjo.eja.eja_algebra import random_eja
 882
 883         TESTS::
 884
 885             sage: set_random_seed()
 886             sage: J = random_eja()
 887             sage: x = J.random_element()
 888             sage: y = J.random_element()
 889             sage: x.operator()(y) == x*y
 890             True
 891             sage: y.operator()(x) == x*y
 892             True
 893
 894         """
 895         P = self.parent()
 896         left_mult_by_self = lambda y: self*y
 897         L = P.module_morphism(function=left_mult_by_self, codomain=P)
 898         return FiniteDimensionalEuclideanJordanAlgebraOperator(
 899                  P,
 900                  P,
 901                  L.matrix() )
 902
 903
 904     def quadratic_representation(self, other=None):
 905         """
 906         Return the quadratic representation of this element.
 907
 908         SETUP::
 909
 910             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 911             ....:                                  random_eja)
 912
 913         EXAMPLES:
 914
 915         The explicit form in the spin factor algebra is given by
 916         Alizadeh's Example 11.12::
 917
 918             sage: set_random_seed()
 919             sage: n = ZZ.random_element(1,10)
 920             sage: J = JordanSpinEJA(n)
 921             sage: x = J.random_element()
 922             sage: x_vec = x.to_vector()
 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 = J.random_element()
 940             sage: y = J.random_element()
 941             sage: Lx = x.operator()
 942             sage: Lxx = (x*x).operator()
 943             sage: Qx = x.quadratic_representation()
 944             sage: Qy = y.quadratic_representation()
 945             sage: Qxy = x.quadratic_representation(y)
 946             sage: Qex = J.one().quadratic_representation(x)
 947             sage: n = ZZ.random_element(10)
 948             sage: Qxn = (x^n).quadratic_representation()
 949
 950         Property 1:
 951
 952             sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
 953             True
 954
 955         Property 2 (multiply on the right for :trac:`28272`):
 956
 957             sage: alpha = J.base_ring().random_element()
 958             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
 959             True
 960
 961         Property 3:
 962
 963             sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
 964             True
 965
 966             sage: not x.is_invertible() or (
 967             ....:   ~Qx
 968             ....:   ==
 969             ....:   x.inverse().quadratic_representation() )
 970             True
 971
 972             sage: Qxy(J.one()) == x*y
 973             True
 974
 975         Property 4:
 976
 977             sage: not x.is_invertible() or (
 978             ....:   x.quadratic_representation(x.inverse())*Qx
 979             ....:   == Qx*x.quadratic_representation(x.inverse()) )
 980             True
 981
 982             sage: not x.is_invertible() or (
 983             ....:   x.quadratic_representation(x.inverse())*Qx
 984             ....:   ==
 985             ....:   2*x.operator()*Qex - Qx )
 986             True
 987
 988             sage: 2*x.operator()*Qex - Qx == Lxx
 989             True
 990
 991         Property 5:
 992
 993             sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
 994             True
 995
 996         Property 6:
 997
 998             sage: Qxn == (Qx)^n
 999             True
1000
1001         Property 7:
1002
1003             sage: not x.is_invertible() or (
1004             ....:   Qx*x.inverse().operator() == Lx )
1005             True
1006
1007         Property 8:
1008
1009             sage: not x.operator_commutes_with(y) or (
1010             ....:   Qx(y)^n == Qxn(y^n) )
1011             True
1012
1013         """
1014         if other is None:
1015             other=self
1016         elif not other in self.parent():
1017             raise TypeError("'other' must live in the same algebra")
1018
1019         L = self.operator()
1020         M = other.operator()
1021         return ( L*M + M*L - (self*other).operator() )
1022
1023
1024
1025
1026     def subalgebra_generated_by(self):
1027         """
1028         Return the associative subalgebra of the parent EJA generated
1029         by this element.
1030
1031         SETUP::
1032
1033             sage: from mjo.eja.eja_algebra import random_eja
1034
1035         TESTS:
1036
1037         This subalgebra, being composed of only powers, is associative::
1038
1039             sage: set_random_seed()
1040             sage: x0 = random_eja().random_element()
1041             sage: A = x0.subalgebra_generated_by()
1042             sage: x = A.random_element()
1043             sage: y = A.random_element()
1044             sage: z = A.random_element()
1045             sage: (x*y)*z == x*(y*z)
1046             True
1047
1048         Squaring in the subalgebra should work the same as in
1049         the superalgebra::
1050
1051             sage: set_random_seed()
1052             sage: x = random_eja().random_element()
1053             sage: A = x.subalgebra_generated_by()
1054             sage: A(x^2) == A(x)*A(x)
1055             True
1056
1057         The subalgebra generated by the zero element is trivial::
1058
1059             sage: set_random_seed()
1060             sage: A = random_eja().zero().subalgebra_generated_by()
1061             sage: A
1062             Euclidean Jordan algebra of dimension 0 over...
1063             sage: A.one()
1064             0
1065
1066         """
1067         return FiniteDimensionalEuclideanJordanElementSubalgebra(self)
1068
1069
1070     def subalgebra_idempotent(self):
1071         """
1072         Find an idempotent in the associative subalgebra I generate
1073         using Proposition 2.3.5 in Baes.
1074
1075         SETUP::
1076
1077             sage: from mjo.eja.eja_algebra import random_eja
1078
1079         TESTS::
1080
1081             sage: set_random_seed()
1082             sage: J = random_eja()
1083             sage: x = J.random_element()
1084             sage: while x.is_nilpotent():
1085             ....:     x = J.random_element()
1086             sage: c = x.subalgebra_idempotent()
1087             sage: c^2 == c
1088             True
1089
1090         """
1091         if self.is_nilpotent():
1092             raise ValueError("this only works with non-nilpotent elements!")
1093
1094         J = self.subalgebra_generated_by()
1095         u = J(self)
1096
1097         # The image of the matrix of left-u^m-multiplication
1098         # will be minimal for some natural number s...
1099         s = 0
1100         minimal_dim = J.dimension()
1101         for i in xrange(1, minimal_dim):
1102             this_dim = (u**i).operator().matrix().image().dimension()
1103             if this_dim < minimal_dim:
1104                 minimal_dim = this_dim
1105                 s = i
1106
1107         # Now minimal_matrix should correspond to the smallest
1108         # non-zero subspace in Baes's (or really, Koecher's)
1109         # proposition.
1110         #
1111         # However, we need to restrict the matrix to work on the
1112         # subspace... or do we? Can't we just solve, knowing that
1113         # A(c) = u^(s+1) should have a solution in the big space,
1114         # too?
1115         #
1116         # Beware, solve_right() means that we're using COLUMN vectors.
1117         # Our FiniteDimensionalAlgebraElement superclass uses rows.
1118         u_next = u**(s+1)
1119         A = u_next.operator().matrix()
1120         c = J.from_vector(A.solve_right(u_next.to_vector()))
1121
1122         # Now c is the idempotent we want, but it still lives in the subalgebra.
1123         return c.superalgebra_element()
1124
1125
1126     def trace(self):
1127         """
1128         Return my trace, the sum of my eigenvalues.
1129
1130         SETUP::
1131
1132             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1133             ....:                                  RealCartesianProductEJA,
1134             ....:                                  random_eja)
1135
1136         EXAMPLES::
1137
1138             sage: J = JordanSpinEJA(3)
1139             sage: x = sum(J.gens())
1140             sage: x.trace()
1141             2
1142
1143         ::
1144
1145             sage: J = RealCartesianProductEJA(5)
1146             sage: J.one().trace()
1147             5
1148
1149         TESTS:
1150
1151         The trace of an element is a real number::
1152
1153             sage: set_random_seed()
1154             sage: J = random_eja()
1155             sage: J.random_element().trace() in J.base_ring()
1156             True
1157
1158         """
1159         P = self.parent()
1160         r = P.rank()
1161         p = P._charpoly_coeff(r-1)
1162         # The _charpoly_coeff function already adds the factor of
1163         # -1 to ensure that _charpoly_coeff(r-1) is really what
1164         # appears in front of t^{r-1} in the charpoly. However,
1165         # we want the negative of THAT for the trace.
1166         return -p(*self.to_vector())
1167
1168
1169     def trace_inner_product(self, other):
1170         """
1171         Return the trace inner product of myself and ``other``.
1172
1173         SETUP::
1174
1175             sage: from mjo.eja.eja_algebra import random_eja
1176
1177         TESTS:
1178
1179         The trace inner product is commutative, bilinear, and satisfies
1180         the Jordan axiom:
1181
1182             sage: set_random_seed()
1183             sage: J = random_eja()
1184             sage: x = J.random_element();
1185             sage: y = J.random_element()
1186             sage: z = J.random_element()
1187             sage: # commutative
1188             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1189             True
1190             sage: # bilinear
1191             sage: a = J.base_ring().random_element();
1192             sage: actual = (a*(x+z)).trace_inner_product(y)
1193             sage: expected = ( a*x.trace_inner_product(y) +
1194             ....:              a*z.trace_inner_product(y) )
1195             sage: actual == expected
1196             True
1197             sage: actual = x.trace_inner_product(a*(y+z))
1198             sage: expected = ( a*x.trace_inner_product(y) +
1199             ....:              a*x.trace_inner_product(z) )
1200             sage: actual == expected
1201             True
1202             sage: # jordan axiom
1203             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1204             True
1205
1206         """
1207         if not other in self.parent():
1208             raise TypeError("'other' must live in the same algebra")
1209
1210         return (self*other).trace()
1211
1212
1213     def trace_norm(self):
1214         """
1215         The norm of this element with respect to :meth:`trace_inner_product`.
1216
1217         SETUP::
1218
1219             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1220             ....:                                  RealCartesianProductEJA)
1221
1222         EXAMPLES::
1223
1224             sage: J = RealCartesianProductEJA(2)
1225             sage: x = sum(J.gens())
1226             sage: x.trace_norm()
1227             sqrt(2)
1228
1229         ::
1230
1231             sage: J = JordanSpinEJA(4)
1232             sage: x = sum(J.gens())
1233             sage: x.trace_norm()
1234             2*sqrt(2)
1235
1236         """
1237         return self.trace_inner_product(self).sqrt()