mjo/eja/eja_element.py

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