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_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_minimal_idempotent(self):
 550         """
 551         Return whether or not this element is a minimal idempotent.
 552
 553
 554         An element of a Euclidean Jordan algebra is a minimal idempotent
 555         if it :meth:`is_idempotent` and if its Peirce subalgebra
 556         corresponding to the eigenvalue ``1`` has dimension ``1`` (Baes,
 557         Proposition 2.7.17).
 558
 559         SETUP::
 560
 561             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 562             ....:                                  RealSymmetricEJA,
 563             ....:                                  random_eja)
 564
 565         WARNING::
 566
 567         This method is sloooooow.
 568
 569         EXAMPLES:
 570
 571         The spectral decomposition of a non-regular element should always
 572         contain at least one non-minimal idempotent::
 573
 574             sage: J = RealSymmetricEJA(3, AA)
 575             sage: x = sum(J.gens())
 576             sage: x.is_regular()
 577             False
 578             sage: [ c.is_minimal_idempotent()
 579             ....:   for (l,c) in x.spectral_decomposition() ]
 580             [False, True]
 581
 582         On the other hand, the spectral decomposition of a regular
 583         element should always be in terms of minimal idempotents::
 584
 585             sage: J = JordanSpinEJA(4, AA)
 586             sage: x = sum( i*J.gens()[i] for i in range(len(J.gens())) )
 587             sage: x.is_regular()
 588             True
 589             sage: [ c.is_minimal_idempotent()
 590             ....:   for (l,c) in x.spectral_decomposition() ]
 591             [True, True]
 592
 593         TESTS:
 594
 595         The identity element is minimal only in an EJA of rank one::
 596
 597             sage: set_random_seed()
 598             sage: J = random_eja()
 599             sage: J.rank() == 1 or not J.one().is_minimal_idempotent()
 600             True
 601
 602         A non-idempotent cannot be a minimal idempotent::
 603
 604             sage: set_random_seed()
 605             sage: J = JordanSpinEJA(4)
 606             sage: x = J.random_element()
 607             sage: (not x.is_idempotent()) and x.is_minimal_idempotent()
 608             False
 609
 610         Proposition 2.7.19 in Baes says that an element is a minimal
 611         idempotent if and only if it's idempotent with trace equal to
 612         unity::
 613
 614             sage: set_random_seed()
 615             sage: J = JordanSpinEJA(4)
 616             sage: x = J.random_element()
 617             sage: expected = (x.is_idempotent() and x.trace() == 1)
 618             sage: actual = x.is_minimal_idempotent()
 619             sage: actual == expected
 620             True
 621
 622         """
 623         # TODO: when the Peirce decomposition is implemented for real,
 624         # we can use that instead of finding this eigenspace manually.
 625         #
 626         # Trivial eigenspaces don't appear in the list, so we default to the
 627         # trivial one and override it if there's a nontrivial space in the
 628         # list.
 629         if not self.is_idempotent():
 630             return False
 631
 632         J1 = VectorSpace(self.parent().base_ring(), 0)
 633         for (eigval, eigspace) in self.operator().matrix().left_eigenspaces():
 634             if eigval == 1:
 635                 J1 = eigspace
 636         return (J1.dimension() == 1)
 637
 638
 639     def is_nilpotent(self):
 640         """
 641         Return whether or not some power of this element is zero.
 642
 643         ALGORITHM:
 644
 645         We use Theorem 5 in Chapter III of Koecher, which says that
 646         an element ``x`` is nilpotent if and only if ``x.operator()``
 647         is nilpotent. And it is a basic fact of linear algebra that
 648         an operator on an `n`-dimensional space is nilpotent if and
 649         only if, when raised to the `n`th power, it equals the zero
 650         operator (for example, see Axler Corollary 8.8).
 651
 652         SETUP::
 653
 654             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 655             ....:                                  random_eja)
 656
 657         EXAMPLES::
 658
 659             sage: J = JordanSpinEJA(3)
 660             sage: x = sum(J.gens())
 661             sage: x.is_nilpotent()
 662             False
 663
 664         TESTS:
 665
 666         The identity element is never nilpotent, except in a trivial EJA::
 667
 668             sage: set_random_seed()
 669             sage: J = random_eja()
 670             sage: J.one().is_nilpotent() and not J.is_trivial()
 671             False
 672
 673         The additive identity is always nilpotent::
 674
 675             sage: set_random_seed()
 676             sage: random_eja().zero().is_nilpotent()
 677             True
 678
 679         """
 680         P = self.parent()
 681         zero_operator = P.zero().operator()
 682         return self.operator()**P.dimension() == zero_operator
 683
 684
 685     def is_regular(self):
 686         """
 687         Return whether or not this is a regular element.
 688
 689         SETUP::
 690
 691             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 692             ....:                                  random_eja)
 693
 694         EXAMPLES:
 695
 696         The identity element always has degree one, but any element
 697         linearly-independent from it is regular::
 698
 699             sage: J = JordanSpinEJA(5)
 700             sage: J.one().is_regular()
 701             False
 702             sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 703             sage: for x in J.gens():
 704             ....:     (J.one() + x).is_regular()
 705             False
 706             True
 707             True
 708             True
 709             True
 710
 711         TESTS:
 712
 713         The zero element should never be regular, unless the parent
 714         algebra has dimension less than or equal to one::
 715
 716             sage: set_random_seed()
 717             sage: J = random_eja()
 718             sage: J.dimension() <= 1 or not J.zero().is_regular()
 719             True
 720
 721         The unit element isn't regular unless the algebra happens to
 722         consist of only its scalar multiples::
 723
 724             sage: set_random_seed()
 725             sage: J = random_eja()
 726             sage: J.dimension() <= 1 or not J.one().is_regular()
 727             True
 728
 729         """
 730         return self.degree() == self.parent().rank()
 731
 732
 733     def degree(self):
 734         """
 735         Return the degree of this element, which is defined to be
 736         the degree of its minimal polynomial.
 737
 738         ALGORITHM:
 739
 740         For now, we skip the messy minimal polynomial computation
 741         and instead return the dimension of the vector space spanned
 742         by the powers of this element. The latter is a bit more
 743         straightforward to compute.
 744
 745         SETUP::
 746
 747             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 748             ....:                                  random_eja)
 749
 750         EXAMPLES::
 751
 752             sage: J = JordanSpinEJA(4)
 753             sage: J.one().degree()
 754             1
 755             sage: e0,e1,e2,e3 = J.gens()
 756             sage: (e0 - e1).degree()
 757             2
 758
 759         In the spin factor algebra (of rank two), all elements that
 760         aren't multiples of the identity are regular::
 761
 762             sage: set_random_seed()
 763             sage: J = JordanSpinEJA.random_instance()
 764             sage: x = J.random_element()
 765             sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 766             True
 767
 768         TESTS:
 769
 770         The zero and unit elements are both of degree one in nontrivial
 771         algebras::
 772
 773             sage: set_random_seed()
 774             sage: J = random_eja()
 775             sage: d = J.zero().degree()
 776             sage: (J.is_trivial() and d == 0) or d == 1
 777             True
 778             sage: d = J.one().degree()
 779             sage: (J.is_trivial() and d == 0) or d == 1
 780             True
 781
 782         Our implementation agrees with the definition::
 783
 784             sage: set_random_seed()
 785             sage: x = random_eja().random_element()
 786             sage: x.degree() == x.minimal_polynomial().degree()
 787             True
 788
 789         """
 790         if self.is_zero() and not self.parent().is_trivial():
 791             # The minimal polynomial of zero in a nontrivial algebra
 792             # is "t"; in a trivial algebra it's "1" by convention
 793             # (it's an empty product).
 794             return 1
 795         return self.subalgebra_generated_by().dimension()
 796
 797
 798     def left_matrix(self):
 799         """
 800         Our parent class defines ``left_matrix`` and ``matrix``
 801         methods whose names are misleading. We don't want them.
 802         """
 803         raise NotImplementedError("use operator().matrix() instead")
 804
 805     matrix = left_matrix
 806
 807
 808     def minimal_polynomial(self):
 809         """
 810         Return the minimal polynomial of this element,
 811         as a function of the variable `t`.
 812
 813         ALGORITHM:
 814
 815         We restrict ourselves to the associative subalgebra
 816         generated by this element, and then return the minimal
 817         polynomial of this element's operator matrix (in that
 818         subalgebra). This works by Baes Proposition 2.3.16.
 819
 820         SETUP::
 821
 822             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 823             ....:                                  RealSymmetricEJA,
 824             ....:                                  TrivialEJA,
 825             ....:                                  random_eja)
 826
 827         EXAMPLES:
 828
 829         Keeping in mind that the polynomial ``1`` evaluates the identity
 830         element (also the zero element) of the trivial algebra, it is clear
 831         that the polynomial ``1`` is the minimal polynomial of the only
 832         element in a trivial algebra::
 833
 834             sage: J = TrivialEJA()
 835             sage: J.one().minimal_polynomial()
 836             1
 837             sage: J.zero().minimal_polynomial()
 838             1
 839
 840         TESTS:
 841
 842         The minimal polynomial of the identity and zero elements are
 843         always the same::
 844
 845             sage: set_random_seed()
 846             sage: J = random_eja(nontrivial=True)
 847             sage: J.one().minimal_polynomial()
 848             t - 1
 849             sage: J.zero().minimal_polynomial()
 850             t
 851
 852         The degree of an element is (by one definition) the degree
 853         of its minimal polynomial::
 854
 855             sage: set_random_seed()
 856             sage: x = random_eja().random_element()
 857             sage: x.degree() == x.minimal_polynomial().degree()
 858             True
 859
 860         The minimal polynomial and the characteristic polynomial coincide
 861         and are known (see Alizadeh, Example 11.11) for all elements of
 862         the spin factor algebra that aren't scalar multiples of the
 863         identity. We require the dimension of the algebra to be at least
 864         two here so that said elements actually exist::
 865
 866             sage: set_random_seed()
 867             sage: n_max = max(2, JordanSpinEJA._max_test_case_size())
 868             sage: n = ZZ.random_element(2, n_max)
 869             sage: J = JordanSpinEJA(n)
 870             sage: y = J.random_element()
 871             sage: while y == y.coefficient(0)*J.one():
 872             ....:     y = J.random_element()
 873             sage: y0 = y.to_vector()[0]
 874             sage: y_bar = y.to_vector()[1:]
 875             sage: actual = y.minimal_polynomial()
 876             sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
 877             sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
 878             sage: bool(actual == expected)
 879             True
 880
 881         The minimal polynomial should always kill its element::
 882
 883             sage: set_random_seed()
 884             sage: x = random_eja().random_element()
 885             sage: p = x.minimal_polynomial()
 886             sage: x.apply_univariate_polynomial(p)
 887             0
 888
 889         The minimal polynomial is invariant under a change of basis,
 890         and in particular, a re-scaling of the basis::
 891
 892             sage: set_random_seed()
 893             sage: n_max = RealSymmetricEJA._max_test_case_size()
 894             sage: n = ZZ.random_element(1, n_max)
 895             sage: J1 = RealSymmetricEJA(n,QQ)
 896             sage: J2 = RealSymmetricEJA(n,QQ,normalize_basis=False)
 897             sage: X = random_matrix(QQ,n)
 898             sage: X = X*X.transpose()
 899             sage: x1 = J1(X)
 900             sage: x2 = J2(X)
 901             sage: x1.minimal_polynomial() == x2.minimal_polynomial()
 902             True
 903
 904         """
 905         if self.is_zero():
 906             # We would generate a zero-dimensional subalgebra
 907             # where the minimal polynomial would be constant.
 908             # That might be correct, but only if *this* algebra
 909             # is trivial too.
 910             if not self.parent().is_trivial():
 911                 # Pretty sure we know what the minimal polynomial of
 912                 # the zero operator is going to be. This ensures
 913                 # consistency of e.g. the polynomial variable returned
 914                 # in the "normal" case without us having to think about it.
 915                 return self.operator().minimal_polynomial()
 916
 917         A = self.subalgebra_generated_by()
 918         return A(self).operator().minimal_polynomial()
 919
 920
 921
 922     def natural_representation(self):
 923         """
 924         Return a more-natural representation of this element.
 925
 926         Every finite-dimensional Euclidean Jordan Algebra is a
 927         direct sum of five simple algebras, four of which comprise
 928         Hermitian matrices. This method returns the original
 929         "natural" representation of this element as a Hermitian
 930         matrix, if it has one. If not, you get the usual representation.
 931
 932         SETUP::
 933
 934             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 935             ....:                                  QuaternionHermitianEJA)
 936
 937         EXAMPLES::
 938
 939             sage: J = ComplexHermitianEJA(3)
 940             sage: J.one()
 941             e0 + e3 + e8
 942             sage: J.one().natural_representation()
 943             [1 0 0 0 0 0]
 944             [0 1 0 0 0 0]
 945             [0 0 1 0 0 0]
 946             [0 0 0 1 0 0]
 947             [0 0 0 0 1 0]
 948             [0 0 0 0 0 1]
 949
 950         ::
 951
 952             sage: J = QuaternionHermitianEJA(3)
 953             sage: J.one()
 954             e0 + e5 + e14
 955             sage: J.one().natural_representation()
 956             [1 0 0 0 0 0 0 0 0 0 0 0]
 957             [0 1 0 0 0 0 0 0 0 0 0 0]
 958             [0 0 1 0 0 0 0 0 0 0 0 0]
 959             [0 0 0 1 0 0 0 0 0 0 0 0]
 960             [0 0 0 0 1 0 0 0 0 0 0 0]
 961             [0 0 0 0 0 1 0 0 0 0 0 0]
 962             [0 0 0 0 0 0 1 0 0 0 0 0]
 963             [0 0 0 0 0 0 0 1 0 0 0 0]
 964             [0 0 0 0 0 0 0 0 1 0 0 0]
 965             [0 0 0 0 0 0 0 0 0 1 0 0]
 966             [0 0 0 0 0 0 0 0 0 0 1 0]
 967             [0 0 0 0 0 0 0 0 0 0 0 1]
 968
 969         """
 970         B = self.parent().natural_basis()
 971         W = self.parent().natural_basis_space()
 972         return W.linear_combination(izip(B,self.to_vector()))
 973
 974
 975     def norm(self):
 976         """
 977         The norm of this element with respect to :meth:`inner_product`.
 978
 979         SETUP::
 980
 981             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 982             ....:                                  RealCartesianProductEJA)
 983
 984         EXAMPLES::
 985
 986             sage: J = RealCartesianProductEJA(2)
 987             sage: x = sum(J.gens())
 988             sage: x.norm()
 989             sqrt(2)
 990
 991         ::
 992
 993             sage: J = JordanSpinEJA(4)
 994             sage: x = sum(J.gens())
 995             sage: x.norm()
 996             2
 997
 998         """
 999         return self.inner_product(self).sqrt()
1000
1001
1002     def operator(self):
1003         """
1004         Return the left-multiplication-by-this-element
1005         operator on the ambient algebra.
1006
1007         SETUP::
1008
1009             sage: from mjo.eja.eja_algebra import random_eja
1010
1011         TESTS::
1012
1013             sage: set_random_seed()
1014             sage: J = random_eja()
1015             sage: x,y = J.random_elements(2)
1016             sage: x.operator()(y) == x*y
1017             True
1018             sage: y.operator()(x) == x*y
1019             True
1020
1021         """
1022         P = self.parent()
1023         left_mult_by_self = lambda y: self*y
1024         L = P.module_morphism(function=left_mult_by_self, codomain=P)
1025         return FiniteDimensionalEuclideanJordanAlgebraOperator(
1026                  P,
1027                  P,
1028                  L.matrix() )
1029
1030
1031     def quadratic_representation(self, other=None):
1032         """
1033         Return the quadratic representation of this element.
1034
1035         SETUP::
1036
1037             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1038             ....:                                  random_eja)
1039
1040         EXAMPLES:
1041
1042         The explicit form in the spin factor algebra is given by
1043         Alizadeh's Example 11.12::
1044
1045             sage: set_random_seed()
1046             sage: x = JordanSpinEJA.random_instance().random_element()
1047             sage: x_vec = x.to_vector()
1048             sage: n = x_vec.degree()
1049             sage: x0 = x_vec[0]
1050             sage: x_bar = x_vec[1:]
1051             sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1052             sage: B = 2*x0*x_bar.row()
1053             sage: C = 2*x0*x_bar.column()
1054             sage: D = matrix.identity(QQ, n-1)
1055             sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1056             sage: D = D + 2*x_bar.tensor_product(x_bar)
1057             sage: Q = matrix.block(2,2,[A,B,C,D])
1058             sage: Q == x.quadratic_representation().matrix()
1059             True
1060
1061         Test all of the properties from Theorem 11.2 in Alizadeh::
1062
1063             sage: set_random_seed()
1064             sage: J = random_eja()
1065             sage: x,y = J.random_elements(2)
1066             sage: Lx = x.operator()
1067             sage: Lxx = (x*x).operator()
1068             sage: Qx = x.quadratic_representation()
1069             sage: Qy = y.quadratic_representation()
1070             sage: Qxy = x.quadratic_representation(y)
1071             sage: Qex = J.one().quadratic_representation(x)
1072             sage: n = ZZ.random_element(10)
1073             sage: Qxn = (x^n).quadratic_representation()
1074
1075         Property 1:
1076
1077             sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1078             True
1079
1080         Property 2 (multiply on the right for :trac:`28272`):
1081
1082             sage: alpha = J.base_ring().random_element()
1083             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
1084             True
1085
1086         Property 3:
1087
1088             sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
1089             True
1090
1091             sage: not x.is_invertible() or (
1092             ....:   ~Qx
1093             ....:   ==
1094             ....:   x.inverse().quadratic_representation() )
1095             True
1096
1097             sage: Qxy(J.one()) == x*y
1098             True
1099
1100         Property 4:
1101
1102             sage: not x.is_invertible() or (
1103             ....:   x.quadratic_representation(x.inverse())*Qx
1104             ....:   == Qx*x.quadratic_representation(x.inverse()) )
1105             True
1106
1107             sage: not x.is_invertible() or (
1108             ....:   x.quadratic_representation(x.inverse())*Qx
1109             ....:   ==
1110             ....:   2*Lx*Qex - Qx )
1111             True
1112
1113             sage: 2*Lx*Qex - Qx == Lxx
1114             True
1115
1116         Property 5:
1117
1118             sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1119             True
1120
1121         Property 6:
1122
1123             sage: Qxn == (Qx)^n
1124             True
1125
1126         Property 7:
1127
1128             sage: not x.is_invertible() or (
1129             ....:   Qx*x.inverse().operator() == Lx )
1130             True
1131
1132         Property 8:
1133
1134             sage: not x.operator_commutes_with(y) or (
1135             ....:   Qx(y)^n == Qxn(y^n) )
1136             True
1137
1138         """
1139         if other is None:
1140             other=self
1141         elif not other in self.parent():
1142             raise TypeError("'other' must live in the same algebra")
1143
1144         L = self.operator()
1145         M = other.operator()
1146         return ( L*M + M*L - (self*other).operator() )
1147
1148
1149
1150     def spectral_decomposition(self):
1151         """
1152         Return the unique spectral decomposition of this element.
1153
1154         ALGORITHM:
1155
1156         Following Faraut and Korányi's Theorem III.1.1, we restrict this
1157         element's left-multiplication-by operator to the subalgebra it
1158         generates. We then compute the spectral decomposition of that
1159         operator, and the spectral projectors we get back must be the
1160         left-multiplication-by operators for the idempotents we
1161         seek. Thus applying them to the identity element gives us those
1162         idempotents.
1163
1164         Since the eigenvalues are required to be distinct, we take
1165         the spectral decomposition of the zero element to be zero
1166         times the identity element of the algebra (which is idempotent,
1167         obviously).
1168
1169         SETUP::
1170
1171             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1172
1173         EXAMPLES:
1174
1175         The spectral decomposition of the identity is ``1`` times itself,
1176         and the spectral decomposition of zero is ``0`` times the identity::
1177
1178             sage: J = RealSymmetricEJA(3,AA)
1179             sage: J.one()
1180             e0 + e2 + e5
1181             sage: J.one().spectral_decomposition()
1182             [(1, e0 + e2 + e5)]
1183             sage: J.zero().spectral_decomposition()
1184             [(0, e0 + e2 + e5)]
1185
1186         TESTS::
1187
1188             sage: J = RealSymmetricEJA(4,AA)
1189             sage: x = sum(J.gens())
1190             sage: sd = x.spectral_decomposition()
1191             sage: l0 = sd[0][0]
1192             sage: l1 = sd[1][0]
1193             sage: c0 = sd[0][1]
1194             sage: c1 = sd[1][1]
1195             sage: c0.inner_product(c1) == 0
1196             True
1197             sage: c0.is_idempotent()
1198             True
1199             sage: c1.is_idempotent()
1200             True
1201             sage: c0 + c1 == J.one()
1202             True
1203             sage: l0*c0 + l1*c1 == x
1204             True
1205
1206         """
1207         P = self.parent()
1208         A = self.subalgebra_generated_by(orthonormalize_basis=True)
1209         result = []
1210         for (evalue, proj) in A(self).operator().spectral_decomposition():
1211             result.append( (evalue, proj(A.one()).superalgebra_element()) )
1212         return result
1213
1214     def subalgebra_generated_by(self, orthonormalize_basis=False):
1215         """
1216         Return the associative subalgebra of the parent EJA generated
1217         by this element.
1218
1219         Since our parent algebra is unital, we want "subalgebra" to mean
1220         "unital subalgebra" as well; thus the subalgebra that an element
1221         generates will itself be a Euclidean Jordan algebra after
1222         restricting the algebra operations appropriately. This is the
1223         subalgebra that Faraut and Korányi work with in section II.2, for
1224         example.
1225
1226         SETUP::
1227
1228             sage: from mjo.eja.eja_algebra import random_eja
1229
1230         TESTS:
1231
1232         This subalgebra, being composed of only powers, is associative::
1233
1234             sage: set_random_seed()
1235             sage: x0 = random_eja().random_element()
1236             sage: A = x0.subalgebra_generated_by()
1237             sage: x,y,z = A.random_elements(3)
1238             sage: (x*y)*z == x*(y*z)
1239             True
1240
1241         Squaring in the subalgebra should work the same as in
1242         the superalgebra::
1243
1244             sage: set_random_seed()
1245             sage: x = random_eja().random_element()
1246             sage: A = x.subalgebra_generated_by()
1247             sage: A(x^2) == A(x)*A(x)
1248             True
1249
1250         By definition, the subalgebra generated by the zero element is
1251         the one-dimensional algebra generated by the identity
1252         element... unless the original algebra was trivial, in which
1253         case the subalgebra is trivial too::
1254
1255             sage: set_random_seed()
1256             sage: A = random_eja().zero().subalgebra_generated_by()
1257             sage: (A.is_trivial() and A.dimension() == 0) or A.dimension() == 1
1258             True
1259
1260         """
1261         return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
1262
1263
1264     def subalgebra_idempotent(self):
1265         """
1266         Find an idempotent in the associative subalgebra I generate
1267         using Proposition 2.3.5 in Baes.
1268
1269         SETUP::
1270
1271             sage: from mjo.eja.eja_algebra import random_eja
1272
1273         TESTS::
1274
1275             sage: set_random_seed()
1276             sage: J = random_eja()
1277             sage: x = J.random_element()
1278             sage: while x.is_nilpotent():
1279             ....:     x = J.random_element()
1280             sage: c = x.subalgebra_idempotent()
1281             sage: c^2 == c
1282             True
1283
1284         """
1285         if self.is_nilpotent():
1286             raise ValueError("this only works with non-nilpotent elements!")
1287
1288         J = self.subalgebra_generated_by()
1289         u = J(self)
1290
1291         # The image of the matrix of left-u^m-multiplication
1292         # will be minimal for some natural number s...
1293         s = 0
1294         minimal_dim = J.dimension()
1295         for i in xrange(1, minimal_dim):
1296             this_dim = (u**i).operator().matrix().image().dimension()
1297             if this_dim < minimal_dim:
1298                 minimal_dim = this_dim
1299                 s = i
1300
1301         # Now minimal_matrix should correspond to the smallest
1302         # non-zero subspace in Baes's (or really, Koecher's)
1303         # proposition.
1304         #
1305         # However, we need to restrict the matrix to work on the
1306         # subspace... or do we? Can't we just solve, knowing that
1307         # A(c) = u^(s+1) should have a solution in the big space,
1308         # too?
1309         #
1310         # Beware, solve_right() means that we're using COLUMN vectors.
1311         # Our FiniteDimensionalAlgebraElement superclass uses rows.
1312         u_next = u**(s+1)
1313         A = u_next.operator().matrix()
1314         c = J.from_vector(A.solve_right(u_next.to_vector()))
1315
1316         # Now c is the idempotent we want, but it still lives in the subalgebra.
1317         return c.superalgebra_element()
1318
1319
1320     def trace(self):
1321         """
1322         Return my trace, the sum of my eigenvalues.
1323
1324         In a trivial algebra, however you want to look at it, the trace is
1325         an empty sum for which we declare the result to be zero.
1326
1327         SETUP::
1328
1329             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1330             ....:                                  RealCartesianProductEJA,
1331             ....:                                  TrivialEJA,
1332             ....:                                  random_eja)
1333
1334         EXAMPLES::
1335
1336             sage: J = TrivialEJA()
1337             sage: J.zero().trace()
1338             0
1339
1340         ::
1341             sage: J = JordanSpinEJA(3)
1342             sage: x = sum(J.gens())
1343             sage: x.trace()
1344             2
1345
1346         ::
1347
1348             sage: J = RealCartesianProductEJA(5)
1349             sage: J.one().trace()
1350             5
1351
1352         TESTS:
1353
1354         The trace of an element is a real number::
1355
1356             sage: set_random_seed()
1357             sage: J = random_eja()
1358             sage: J.random_element().trace() in RLF
1359             True
1360
1361         """
1362         P = self.parent()
1363         r = P.rank()
1364
1365         if r == 0:
1366             # Special case for the trivial algebra where
1367             # the trace is an empty sum.
1368             return P.base_ring().zero()
1369
1370         p = P._charpoly_coeff(r-1)
1371         # The _charpoly_coeff function already adds the factor of
1372         # -1 to ensure that _charpoly_coeff(r-1) is really what
1373         # appears in front of t^{r-1} in the charpoly. However,
1374         # we want the negative of THAT for the trace.
1375         return -p(*self.to_vector())
1376
1377
1378     def trace_inner_product(self, other):
1379         """
1380         Return the trace inner product of myself and ``other``.
1381
1382         SETUP::
1383
1384             sage: from mjo.eja.eja_algebra import random_eja
1385
1386         TESTS:
1387
1388         The trace inner product is commutative, bilinear, and associative::
1389
1390             sage: set_random_seed()
1391             sage: J = random_eja()
1392             sage: x,y,z = J.random_elements(3)
1393             sage: # commutative
1394             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1395             True
1396             sage: # bilinear
1397             sage: a = J.base_ring().random_element();
1398             sage: actual = (a*(x+z)).trace_inner_product(y)
1399             sage: expected = ( a*x.trace_inner_product(y) +
1400             ....:              a*z.trace_inner_product(y) )
1401             sage: actual == expected
1402             True
1403             sage: actual = x.trace_inner_product(a*(y+z))
1404             sage: expected = ( a*x.trace_inner_product(y) +
1405             ....:              a*x.trace_inner_product(z) )
1406             sage: actual == expected
1407             True
1408             sage: # associative
1409             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1410             True
1411
1412         """
1413         if not other in self.parent():
1414             raise TypeError("'other' must live in the same algebra")
1415
1416         return (self*other).trace()
1417
1418
1419     def trace_norm(self):
1420         """
1421         The norm of this element with respect to :meth:`trace_inner_product`.
1422
1423         SETUP::
1424
1425             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1426             ....:                                  RealCartesianProductEJA)
1427
1428         EXAMPLES::
1429
1430             sage: J = RealCartesianProductEJA(2)
1431             sage: x = sum(J.gens())
1432             sage: x.trace_norm()
1433             sqrt(2)
1434
1435         ::
1436
1437             sage: J = JordanSpinEJA(4)
1438             sage: x = sum(J.gens())
1439             sage: x.trace_norm()
1440             2*sqrt(2)
1441
1442         """
1443         return self.trace_inner_product(self).sqrt()