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: J = JordanSpinEJA.random_instance()
 428             sage: x = J.random_element()
 429             sage: while not x.is_invertible():
 430             ....:     x = J.random_element()
 431             sage: x_vec = x.to_vector()
 432             sage: x0 = x_vec[0]
 433             sage: x_bar = x_vec[1:]
 434             sage: coeff = ~(x0^2 - x_bar.inner_product(x_bar))
 435             sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
 436             sage: x_inverse = coeff*inv_vec
 437             sage: x.inverse() == J.from_vector(x_inverse)
 438             True
 439
 440         TESTS:
 441
 442         The identity element is its own inverse::
 443
 444             sage: set_random_seed()
 445             sage: J = random_eja()
 446             sage: J.one().inverse() == J.one()
 447             True
 448
 449         If an element has an inverse, it acts like one::
 450
 451             sage: set_random_seed()
 452             sage: J = random_eja()
 453             sage: x = J.random_element()
 454             sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
 455             True
 456
 457         The inverse of the inverse is what we started with::
 458
 459             sage: set_random_seed()
 460             sage: J = random_eja()
 461             sage: x = J.random_element()
 462             sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
 463             True
 464
 465         The zero element is never invertible::
 466
 467             sage: set_random_seed()
 468             sage: J = random_eja().zero().inverse()
 469             Traceback (most recent call last):
 470             ...
 471             ValueError: element is not invertible
 472
 473         """
 474         if not self.is_invertible():
 475             raise ValueError("element is not invertible")
 476
 477         return (~self.quadratic_representation())(self)
 478
 479
 480     def is_invertible(self):
 481         """
 482         Return whether or not this element is invertible.
 483
 484         ALGORITHM:
 485
 486         The usual way to do this is to check if the determinant is
 487         zero, but we need the characteristic polynomial for the
 488         determinant. The minimal polynomial is a lot easier to get,
 489         so we use Corollary 2 in Chapter V of Koecher to check
 490         whether or not the paren't algebra's zero element is a root
 491         of this element's minimal polynomial.
 492
 493         Beware that we can't use the superclass method, because it
 494         relies on the algebra being associative.
 495
 496         SETUP::
 497
 498             sage: from mjo.eja.eja_algebra import random_eja
 499
 500         TESTS:
 501
 502         The identity element is always invertible::
 503
 504             sage: set_random_seed()
 505             sage: J = random_eja()
 506             sage: J.one().is_invertible()
 507             True
 508
 509         The zero element is never invertible in a non-trivial algebra::
 510
 511             sage: set_random_seed()
 512             sage: J = random_eja()
 513             sage: (not J.is_trivial()) and J.zero().is_invertible()
 514             False
 515
 516         """
 517         if self.is_zero():
 518             if self.parent().is_trivial():
 519                 return True
 520             else:
 521                 return False
 522
 523         # In fact, we only need to know if the constant term is non-zero,
 524         # so we can pass in the field's zero element instead.
 525         zero = self.base_ring().zero()
 526         p = self.minimal_polynomial()
 527         return not (p(zero) == zero)
 528
 529
 530     def is_nilpotent(self):
 531         """
 532         Return whether or not some power of this element is zero.
 533
 534         ALGORITHM:
 535
 536         We use Theorem 5 in Chapter III of Koecher, which says that
 537         an element ``x`` is nilpotent if and only if ``x.operator()``
 538         is nilpotent. And it is a basic fact of linear algebra that
 539         an operator on an `n`-dimensional space is nilpotent if and
 540         only if, when raised to the `n`th power, it equals the zero
 541         operator (for example, see Axler Corollary 8.8).
 542
 543         SETUP::
 544
 545             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 546             ....:                                  random_eja)
 547
 548         EXAMPLES::
 549
 550             sage: J = JordanSpinEJA(3)
 551             sage: x = sum(J.gens())
 552             sage: x.is_nilpotent()
 553             False
 554
 555         TESTS:
 556
 557         The identity element is never nilpotent::
 558
 559             sage: set_random_seed()
 560             sage: random_eja().one().is_nilpotent()
 561             False
 562
 563         The additive identity is always nilpotent::
 564
 565             sage: set_random_seed()
 566             sage: random_eja().zero().is_nilpotent()
 567             True
 568
 569         """
 570         P = self.parent()
 571         zero_operator = P.zero().operator()
 572         return self.operator()**P.dimension() == zero_operator
 573
 574
 575     def is_regular(self):
 576         """
 577         Return whether or not this is a regular element.
 578
 579         SETUP::
 580
 581             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 582             ....:                                  random_eja)
 583
 584         EXAMPLES:
 585
 586         The identity element always has degree one, but any element
 587         linearly-independent from it is regular::
 588
 589             sage: J = JordanSpinEJA(5)
 590             sage: J.one().is_regular()
 591             False
 592             sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 593             sage: for x in J.gens():
 594             ....:     (J.one() + x).is_regular()
 595             False
 596             True
 597             True
 598             True
 599             True
 600
 601         TESTS:
 602
 603         The zero element should never be regular, unless the parent
 604         algebra has dimension one::
 605
 606             sage: set_random_seed()
 607             sage: J = random_eja()
 608             sage: J.dimension() == 1 or not J.zero().is_regular()
 609             True
 610
 611         The unit element isn't regular unless the algebra happens to
 612         consist of only its scalar multiples::
 613
 614             sage: set_random_seed()
 615             sage: J = random_eja()
 616             sage: J.dimension() == 1 or not J.one().is_regular()
 617             True
 618
 619         """
 620         return self.degree() == self.parent().rank()
 621
 622
 623     def degree(self):
 624         """
 625         Return the degree of this element, which is defined to be
 626         the degree of its minimal polynomial.
 627
 628         ALGORITHM:
 629
 630         For now, we skip the messy minimal polynomial computation
 631         and instead return the dimension of the vector space spanned
 632         by the powers of this element. The latter is a bit more
 633         straightforward to compute.
 634
 635         SETUP::
 636
 637             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 638             ....:                                  random_eja)
 639
 640         EXAMPLES::
 641
 642             sage: J = JordanSpinEJA(4)
 643             sage: J.one().degree()
 644             1
 645             sage: e0,e1,e2,e3 = J.gens()
 646             sage: (e0 - e1).degree()
 647             2
 648
 649         In the spin factor algebra (of rank two), all elements that
 650         aren't multiples of the identity are regular::
 651
 652             sage: set_random_seed()
 653             sage: J = JordanSpinEJA.random_instance()
 654             sage: x = J.random_element()
 655             sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 656             True
 657
 658         TESTS:
 659
 660         The zero and unit elements are both of degree one::
 661
 662             sage: set_random_seed()
 663             sage: J = random_eja()
 664             sage: J.zero().degree()
 665             1
 666             sage: J.one().degree()
 667             1
 668
 669         Our implementation agrees with the definition::
 670
 671             sage: set_random_seed()
 672             sage: x = random_eja().random_element()
 673             sage: x.degree() == x.minimal_polynomial().degree()
 674             True
 675
 676         """
 677         if self.is_zero() and not self.parent().is_trivial():
 678             # The minimal polynomial of zero in a nontrivial algebra
 679             # is "t"; in a trivial algebra it's "1" by convention
 680             # (it's an empty product).
 681             return 1
 682         return self.subalgebra_generated_by().dimension()
 683
 684
 685     def left_matrix(self):
 686         """
 687         Our parent class defines ``left_matrix`` and ``matrix``
 688         methods whose names are misleading. We don't want them.
 689         """
 690         raise NotImplementedError("use operator().matrix() instead")
 691
 692     matrix = left_matrix
 693
 694
 695     def minimal_polynomial(self):
 696         """
 697         Return the minimal polynomial of this element,
 698         as a function of the variable `t`.
 699
 700         ALGORITHM:
 701
 702         We restrict ourselves to the associative subalgebra
 703         generated by this element, and then return the minimal
 704         polynomial of this element's operator matrix (in that
 705         subalgebra). This works by Baes Proposition 2.3.16.
 706
 707         SETUP::
 708
 709             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 710             ....:                                  RealSymmetricEJA,
 711             ....:                                  random_eja)
 712
 713         TESTS:
 714
 715         The minimal polynomial of the identity and zero elements are
 716         always the same::
 717
 718             sage: set_random_seed()
 719             sage: J = random_eja()
 720             sage: J.one().minimal_polynomial()
 721             t - 1
 722             sage: J.zero().minimal_polynomial()
 723             t
 724
 725         The degree of an element is (by one definition) the degree
 726         of its minimal polynomial::
 727
 728             sage: set_random_seed()
 729             sage: x = random_eja().random_element()
 730             sage: x.degree() == x.minimal_polynomial().degree()
 731             True
 732
 733         The minimal polynomial and the characteristic polynomial coincide
 734         and are known (see Alizadeh, Example 11.11) for all elements of
 735         the spin factor algebra that aren't scalar multiples of the
 736         identity. We require the dimension of the algebra to be at least
 737         two here so that said elements actually exist::
 738
 739             sage: set_random_seed()
 740             sage: n_max = max(2, JordanSpinEJA._max_test_case_size())
 741             sage: n = ZZ.random_element(2, n_max)
 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_max = RealSymmetricEJA._max_test_case_size()
 767             sage: n = ZZ.random_element(1, n_max)
 768             sage: J1 = RealSymmetricEJA(n,QQ)
 769             sage: J2 = RealSymmetricEJA(n,QQ,False)
 770             sage: X = random_matrix(QQ,n)
 771             sage: X = X*X.transpose()
 772             sage: x1 = J1(X)
 773             sage: x2 = J2(X)
 774             sage: x1.minimal_polynomial() == x2.minimal_polynomial()
 775             True
 776
 777         """
 778         if self.is_zero():
 779             # We would generate a zero-dimensional subalgebra
 780             # where the minimal polynomial would be constant.
 781             # That might be correct, but only if *this* algebra
 782             # is trivial too.
 783             if not self.parent().is_trivial():
 784                 # Pretty sure we know what the minimal polynomial of
 785                 # the zero operator is going to be. This ensures
 786                 # consistency of e.g. the polynomial variable returned
 787                 # in the "normal" case without us having to think about it.
 788                 return self.operator().minimal_polynomial()
 789
 790         A = self.subalgebra_generated_by()
 791         return A(self).operator().minimal_polynomial()
 792
 793
 794
 795     def natural_representation(self):
 796         """
 797         Return a more-natural representation of this element.
 798
 799         Every finite-dimensional Euclidean Jordan Algebra is a
 800         direct sum of five simple algebras, four of which comprise
 801         Hermitian matrices. This method returns the original
 802         "natural" representation of this element as a Hermitian
 803         matrix, if it has one. If not, you get the usual representation.
 804
 805         SETUP::
 806
 807             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 808             ....:                                  QuaternionHermitianEJA)
 809
 810         EXAMPLES::
 811
 812             sage: J = ComplexHermitianEJA(3)
 813             sage: J.one()
 814             e0 + e3 + e8
 815             sage: J.one().natural_representation()
 816             [1 0 0 0 0 0]
 817             [0 1 0 0 0 0]
 818             [0 0 1 0 0 0]
 819             [0 0 0 1 0 0]
 820             [0 0 0 0 1 0]
 821             [0 0 0 0 0 1]
 822
 823         ::
 824
 825             sage: J = QuaternionHermitianEJA(3)
 826             sage: J.one()
 827             e0 + e5 + e14
 828             sage: J.one().natural_representation()
 829             [1 0 0 0 0 0 0 0 0 0 0 0]
 830             [0 1 0 0 0 0 0 0 0 0 0 0]
 831             [0 0 1 0 0 0 0 0 0 0 0 0]
 832             [0 0 0 1 0 0 0 0 0 0 0 0]
 833             [0 0 0 0 1 0 0 0 0 0 0 0]
 834             [0 0 0 0 0 1 0 0 0 0 0 0]
 835             [0 0 0 0 0 0 1 0 0 0 0 0]
 836             [0 0 0 0 0 0 0 1 0 0 0 0]
 837             [0 0 0 0 0 0 0 0 1 0 0 0]
 838             [0 0 0 0 0 0 0 0 0 1 0 0]
 839             [0 0 0 0 0 0 0 0 0 0 1 0]
 840             [0 0 0 0 0 0 0 0 0 0 0 1]
 841
 842         """
 843         B = self.parent().natural_basis()
 844         W = self.parent().natural_basis_space()
 845         return W.linear_combination(zip(B,self.to_vector()))
 846
 847
 848     def norm(self):
 849         """
 850         The norm of this element with respect to :meth:`inner_product`.
 851
 852         SETUP::
 853
 854             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 855             ....:                                  RealCartesianProductEJA)
 856
 857         EXAMPLES::
 858
 859             sage: J = RealCartesianProductEJA(2)
 860             sage: x = sum(J.gens())
 861             sage: x.norm()
 862             sqrt(2)
 863
 864         ::
 865
 866             sage: J = JordanSpinEJA(4)
 867             sage: x = sum(J.gens())
 868             sage: x.norm()
 869             2
 870
 871         """
 872         return self.inner_product(self).sqrt()
 873
 874
 875     def operator(self):
 876         """
 877         Return the left-multiplication-by-this-element
 878         operator on the ambient algebra.
 879
 880         SETUP::
 881
 882             sage: from mjo.eja.eja_algebra import random_eja
 883
 884         TESTS::
 885
 886             sage: set_random_seed()
 887             sage: J = random_eja()
 888             sage: x = J.random_element()
 889             sage: y = J.random_element()
 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 = 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()