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