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