mjo/eja/eja_element.py

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