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         TESTS:
 434
 435         The identity element is its own inverse::
 436
 437             sage: set_random_seed()
 438             sage: J = random_eja()
 439             sage: J.one().inverse() == J.one()
 440             True
 441
 442         If an element has an inverse, it acts like one::
 443
 444             sage: set_random_seed()
 445             sage: J = random_eja()
 446             sage: x = J.random_element()
 447             sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
 448             True
 449
 450         The inverse of the inverse is what we started with::
 451
 452             sage: set_random_seed()
 453             sage: J = random_eja()
 454             sage: x = J.random_element()
 455             sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
 456             True
 457
 458         The zero element is never invertible::
 459
 460             sage: set_random_seed()
 461             sage: J = random_eja().zero().inverse()
 462             Traceback (most recent call last):
 463             ...
 464             ValueError: element is not invertible
 465
 466         Proposition II.2.3 in Faraut and Korányi says that the inverse
 467         of an element is the inverse of its left-multiplication operator
 468         applied to the algebra's identity, when that inverse exists::
 469
 470             sage: set_random_seed()
 471             sage: J = random_eja()
 472             sage: x = J.random_element()
 473             sage: (not x.operator().is_invertible()) or (
 474             ....:    x.operator().inverse()(J.one()) == x.inverse() )
 475             True
 476
 477         Proposition II.2.4 in Faraut and Korányi gives a formula for
 478         the inverse based on the characteristic polynomial and the
 479         Cayley-Hamilton theorem for Euclidean Jordan algebras::
 480
 481             sage: set_random_seed()
 482             sage: J = ComplexHermitianEJA(3)
 483             sage: x = J.random_element()
 484             sage: while not x.is_invertible():
 485             ....:     x = J.random_element()
 486             sage: r = J.rank()
 487             sage: a = x.characteristic_polynomial().coefficients(sparse=False)
 488             sage: expected  = (-1)^(r+1)/x.det()
 489             sage: expected *= sum( a[i+1]*x^i for i in range(r) )
 490             sage: x.inverse() == expected
 491             True
 492
 493         """
 494         if not self.is_invertible():
 495             raise ValueError("element is not invertible")
 496
 497         return (~self.quadratic_representation())(self)
 498
 499
 500     def is_invertible(self):
 501         """
 502         Return whether or not this element is invertible.
 503
 504         ALGORITHM:
 505
 506         The usual way to do this is to check if the determinant is
 507         zero, but we need the characteristic polynomial for the
 508         determinant. The minimal polynomial is a lot easier to get,
 509         so we use Corollary 2 in Chapter V of Koecher to check
 510         whether or not the paren't algebra's zero element is a root
 511         of this element's minimal polynomial.
 512
 513         Beware that we can't use the superclass method, because it
 514         relies on the algebra being associative.
 515
 516         SETUP::
 517
 518             sage: from mjo.eja.eja_algebra import random_eja
 519
 520         TESTS:
 521
 522         The identity element is always invertible::
 523
 524             sage: set_random_seed()
 525             sage: J = random_eja()
 526             sage: J.one().is_invertible()
 527             True
 528
 529         The zero element is never invertible in a non-trivial algebra::
 530
 531             sage: set_random_seed()
 532             sage: J = random_eja()
 533             sage: (not J.is_trivial()) and J.zero().is_invertible()
 534             False
 535
 536         """
 537         if self.is_zero():
 538             if self.parent().is_trivial():
 539                 return True
 540             else:
 541                 return False
 542
 543         # In fact, we only need to know if the constant term is non-zero,
 544         # so we can pass in the field's zero element instead.
 545         zero = self.base_ring().zero()
 546         p = self.minimal_polynomial()
 547         return not (p(zero) == zero)
 548
 549
 550     def is_nilpotent(self):
 551         """
 552         Return whether or not some power of this element is zero.
 553
 554         ALGORITHM:
 555
 556         We use Theorem 5 in Chapter III of Koecher, which says that
 557         an element ``x`` is nilpotent if and only if ``x.operator()``
 558         is nilpotent. And it is a basic fact of linear algebra that
 559         an operator on an `n`-dimensional space is nilpotent if and
 560         only if, when raised to the `n`th power, it equals the zero
 561         operator (for example, see Axler Corollary 8.8).
 562
 563         SETUP::
 564
 565             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 566             ....:                                  random_eja)
 567
 568         EXAMPLES::
 569
 570             sage: J = JordanSpinEJA(3)
 571             sage: x = sum(J.gens())
 572             sage: x.is_nilpotent()
 573             False
 574
 575         TESTS:
 576
 577         The identity element is never nilpotent::
 578
 579             sage: set_random_seed()
 580             sage: random_eja().one().is_nilpotent()
 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 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::
 681
 682             sage: set_random_seed()
 683             sage: J = random_eja()
 684             sage: J.zero().degree()
 685             1
 686             sage: J.one().degree()
 687             1
 688
 689         Our implementation agrees with the definition::
 690
 691             sage: set_random_seed()
 692             sage: x = random_eja().random_element()
 693             sage: x.degree() == x.minimal_polynomial().degree()
 694             True
 695
 696         """
 697         if self.is_zero() and not self.parent().is_trivial():
 698             # The minimal polynomial of zero in a nontrivial algebra
 699             # is "t"; in a trivial algebra it's "1" by convention
 700             # (it's an empty product).
 701             return 1
 702         return self.subalgebra_generated_by().dimension()
 703
 704
 705     def left_matrix(self):
 706         """
 707         Our parent class defines ``left_matrix`` and ``matrix``
 708         methods whose names are misleading. We don't want them.
 709         """
 710         raise NotImplementedError("use operator().matrix() instead")
 711
 712     matrix = left_matrix
 713
 714
 715     def minimal_polynomial(self):
 716         """
 717         Return the minimal polynomial of this element,
 718         as a function of the variable `t`.
 719
 720         ALGORITHM:
 721
 722         We restrict ourselves to the associative subalgebra
 723         generated by this element, and then return the minimal
 724         polynomial of this element's operator matrix (in that
 725         subalgebra). This works by Baes Proposition 2.3.16.
 726
 727         SETUP::
 728
 729             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 730             ....:                                  RealSymmetricEJA,
 731             ....:                                  random_eja)
 732
 733         TESTS:
 734
 735         The minimal polynomial of the identity and zero elements are
 736         always the same::
 737
 738             sage: set_random_seed()
 739             sage: J = random_eja()
 740             sage: J.one().minimal_polynomial()
 741             t - 1
 742             sage: J.zero().minimal_polynomial()
 743             t
 744
 745         The degree of an element is (by one definition) the degree
 746         of its minimal polynomial::
 747
 748             sage: set_random_seed()
 749             sage: x = random_eja().random_element()
 750             sage: x.degree() == x.minimal_polynomial().degree()
 751             True
 752
 753         The minimal polynomial and the characteristic polynomial coincide
 754         and are known (see Alizadeh, Example 11.11) for all elements of
 755         the spin factor algebra that aren't scalar multiples of the
 756         identity. We require the dimension of the algebra to be at least
 757         two here so that said elements actually exist::
 758
 759             sage: set_random_seed()
 760             sage: n_max = max(2, JordanSpinEJA._max_test_case_size())
 761             sage: n = ZZ.random_element(2, n_max)
 762             sage: J = JordanSpinEJA(n)
 763             sage: y = J.random_element()
 764             sage: while y == y.coefficient(0)*J.one():
 765             ....:     y = J.random_element()
 766             sage: y0 = y.to_vector()[0]
 767             sage: y_bar = y.to_vector()[1:]
 768             sage: actual = y.minimal_polynomial()
 769             sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
 770             sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
 771             sage: bool(actual == expected)
 772             True
 773
 774         The minimal polynomial should always kill its element::
 775
 776             sage: set_random_seed()
 777             sage: x = random_eja().random_element()
 778             sage: p = x.minimal_polynomial()
 779             sage: x.apply_univariate_polynomial(p)
 780             0
 781
 782         The minimal polynomial is invariant under a change of basis,
 783         and in particular, a re-scaling of the basis::
 784
 785             sage: set_random_seed()
 786             sage: n_max = RealSymmetricEJA._max_test_case_size()
 787             sage: n = ZZ.random_element(1, n_max)
 788             sage: J1 = RealSymmetricEJA(n,QQ)
 789             sage: J2 = RealSymmetricEJA(n,QQ,normalize_basis=False)
 790             sage: X = random_matrix(QQ,n)
 791             sage: X = X*X.transpose()
 792             sage: x1 = J1(X)
 793             sage: x2 = J2(X)
 794             sage: x1.minimal_polynomial() == x2.minimal_polynomial()
 795             True
 796
 797         """
 798         if self.is_zero():
 799             # We would generate a zero-dimensional subalgebra
 800             # where the minimal polynomial would be constant.
 801             # That might be correct, but only if *this* algebra
 802             # is trivial too.
 803             if not self.parent().is_trivial():
 804                 # Pretty sure we know what the minimal polynomial of
 805                 # the zero operator is going to be. This ensures
 806                 # consistency of e.g. the polynomial variable returned
 807                 # in the "normal" case without us having to think about it.
 808                 return self.operator().minimal_polynomial()
 809
 810         A = self.subalgebra_generated_by()
 811         return A(self).operator().minimal_polynomial()
 812
 813
 814
 815     def natural_representation(self):
 816         """
 817         Return a more-natural representation of this element.
 818
 819         Every finite-dimensional Euclidean Jordan Algebra is a
 820         direct sum of five simple algebras, four of which comprise
 821         Hermitian matrices. This method returns the original
 822         "natural" representation of this element as a Hermitian
 823         matrix, if it has one. If not, you get the usual representation.
 824
 825         SETUP::
 826
 827             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 828             ....:                                  QuaternionHermitianEJA)
 829
 830         EXAMPLES::
 831
 832             sage: J = ComplexHermitianEJA(3)
 833             sage: J.one()
 834             e0 + e3 + e8
 835             sage: J.one().natural_representation()
 836             [1 0 0 0 0 0]
 837             [0 1 0 0 0 0]
 838             [0 0 1 0 0 0]
 839             [0 0 0 1 0 0]
 840             [0 0 0 0 1 0]
 841             [0 0 0 0 0 1]
 842
 843         ::
 844
 845             sage: J = QuaternionHermitianEJA(3)
 846             sage: J.one()
 847             e0 + e5 + e14
 848             sage: J.one().natural_representation()
 849             [1 0 0 0 0 0 0 0 0 0 0 0]
 850             [0 1 0 0 0 0 0 0 0 0 0 0]
 851             [0 0 1 0 0 0 0 0 0 0 0 0]
 852             [0 0 0 1 0 0 0 0 0 0 0 0]
 853             [0 0 0 0 1 0 0 0 0 0 0 0]
 854             [0 0 0 0 0 1 0 0 0 0 0 0]
 855             [0 0 0 0 0 0 1 0 0 0 0 0]
 856             [0 0 0 0 0 0 0 1 0 0 0 0]
 857             [0 0 0 0 0 0 0 0 1 0 0 0]
 858             [0 0 0 0 0 0 0 0 0 1 0 0]
 859             [0 0 0 0 0 0 0 0 0 0 1 0]
 860             [0 0 0 0 0 0 0 0 0 0 0 1]
 861
 862         """
 863         B = self.parent().natural_basis()
 864         W = self.parent().natural_basis_space()
 865         return W.linear_combination(izip(B,self.to_vector()))
 866
 867
 868     def norm(self):
 869         """
 870         The norm of this element with respect to :meth:`inner_product`.
 871
 872         SETUP::
 873
 874             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 875             ....:                                  RealCartesianProductEJA)
 876
 877         EXAMPLES::
 878
 879             sage: J = RealCartesianProductEJA(2)
 880             sage: x = sum(J.gens())
 881             sage: x.norm()
 882             sqrt(2)
 883
 884         ::
 885
 886             sage: J = JordanSpinEJA(4)
 887             sage: x = sum(J.gens())
 888             sage: x.norm()
 889             2
 890
 891         """
 892         return self.inner_product(self).sqrt()
 893
 894
 895     def operator(self):
 896         """
 897         Return the left-multiplication-by-this-element
 898         operator on the ambient algebra.
 899
 900         SETUP::
 901
 902             sage: from mjo.eja.eja_algebra import random_eja
 903
 904         TESTS::
 905
 906             sage: set_random_seed()
 907             sage: J = random_eja()
 908             sage: x,y = J.random_elements(2)
 909             sage: x.operator()(y) == x*y
 910             True
 911             sage: y.operator()(x) == x*y
 912             True
 913
 914         """
 915         P = self.parent()
 916         left_mult_by_self = lambda y: self*y
 917         L = P.module_morphism(function=left_mult_by_self, codomain=P)
 918         return FiniteDimensionalEuclideanJordanAlgebraOperator(
 919                  P,
 920                  P,
 921                  L.matrix() )
 922
 923
 924     def quadratic_representation(self, other=None):
 925         """
 926         Return the quadratic representation of this element.
 927
 928         SETUP::
 929
 930             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 931             ....:                                  random_eja)
 932
 933         EXAMPLES:
 934
 935         The explicit form in the spin factor algebra is given by
 936         Alizadeh's Example 11.12::
 937
 938             sage: set_random_seed()
 939             sage: x = JordanSpinEJA.random_instance().random_element()
 940             sage: x_vec = x.to_vector()
 941             sage: n = x_vec.degree()
 942             sage: x0 = x_vec[0]
 943             sage: x_bar = x_vec[1:]
 944             sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 945             sage: B = 2*x0*x_bar.row()
 946             sage: C = 2*x0*x_bar.column()
 947             sage: D = matrix.identity(QQ, n-1)
 948             sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 949             sage: D = D + 2*x_bar.tensor_product(x_bar)
 950             sage: Q = matrix.block(2,2,[A,B,C,D])
 951             sage: Q == x.quadratic_representation().matrix()
 952             True
 953
 954         Test all of the properties from Theorem 11.2 in Alizadeh::
 955
 956             sage: set_random_seed()
 957             sage: J = random_eja()
 958             sage: x,y = J.random_elements(2)
 959             sage: Lx = x.operator()
 960             sage: Lxx = (x*x).operator()
 961             sage: Qx = x.quadratic_representation()
 962             sage: Qy = y.quadratic_representation()
 963             sage: Qxy = x.quadratic_representation(y)
 964             sage: Qex = J.one().quadratic_representation(x)
 965             sage: n = ZZ.random_element(10)
 966             sage: Qxn = (x^n).quadratic_representation()
 967
 968         Property 1:
 969
 970             sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
 971             True
 972
 973         Property 2 (multiply on the right for :trac:`28272`):
 974
 975             sage: alpha = J.base_ring().random_element()
 976             sage: (alpha*x).quadratic_representation() == Qx*(alpha^2)
 977             True
 978
 979         Property 3:
 980
 981             sage: not x.is_invertible() or ( Qx(x.inverse()) == x )
 982             True
 983
 984             sage: not x.is_invertible() or (
 985             ....:   ~Qx
 986             ....:   ==
 987             ....:   x.inverse().quadratic_representation() )
 988             True
 989
 990             sage: Qxy(J.one()) == x*y
 991             True
 992
 993         Property 4:
 994
 995             sage: not x.is_invertible() or (
 996             ....:   x.quadratic_representation(x.inverse())*Qx
 997             ....:   == Qx*x.quadratic_representation(x.inverse()) )
 998             True
 999
1000             sage: not x.is_invertible() or (
1001             ....:   x.quadratic_representation(x.inverse())*Qx
1002             ....:   ==
1003             ....:   2*Lx*Qex - Qx )
1004             True
1005
1006             sage: 2*Lx*Qex - Qx == Lxx
1007             True
1008
1009         Property 5:
1010
1011             sage: Qy(x).quadratic_representation() == Qy*Qx*Qy
1012             True
1013
1014         Property 6:
1015
1016             sage: Qxn == (Qx)^n
1017             True
1018
1019         Property 7:
1020
1021             sage: not x.is_invertible() or (
1022             ....:   Qx*x.inverse().operator() == Lx )
1023             True
1024
1025         Property 8:
1026
1027             sage: not x.operator_commutes_with(y) or (
1028             ....:   Qx(y)^n == Qxn(y^n) )
1029             True
1030
1031         """
1032         if other is None:
1033             other=self
1034         elif not other in self.parent():
1035             raise TypeError("'other' must live in the same algebra")
1036
1037         L = self.operator()
1038         M = other.operator()
1039         return ( L*M + M*L - (self*other).operator() )
1040
1041
1042
1043     def spectral_decomposition(self):
1044         """
1045         Return the unique spectral decomposition of this element.
1046
1047         ALGORITHM:
1048
1049         Following Faraut and Korányi's Theorem III.1.1, we restrict this
1050         element's left-multiplication-by operator to the subalgebra it
1051         generates. We then compute the spectral decomposition of that
1052         operator, and the spectral projectors we get back must be the
1053         left-multiplication-by operators for the idempotents we
1054         seek. Thus applying them to the identity element gives us those
1055         idempotents.
1056
1057         Since the eigenvalues are required to be distinct, we take
1058         the spectral decomposition of the zero element to be zero
1059         times the identity element of the algebra (which is idempotent,
1060         obviously).
1061
1062         SETUP::
1063
1064             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1065
1066         EXAMPLES:
1067
1068         The spectral decomposition of the identity is ``1`` times itself,
1069         and the spectral decomposition of zero is ``0`` times the identity::
1070
1071             sage: J = RealSymmetricEJA(3,AA)
1072             sage: J.one()
1073             e0 + e2 + e5
1074             sage: J.one().spectral_decomposition()
1075             [(1, e0 + e2 + e5)]
1076             sage: J.zero().spectral_decomposition()
1077             [(0, e0 + e2 + e5)]
1078
1079         TESTS::
1080
1081             sage: J = RealSymmetricEJA(4,AA)
1082             sage: x = sum(J.gens())
1083             sage: sd = x.spectral_decomposition()
1084             sage: l0 = sd[0][0]
1085             sage: l1 = sd[1][0]
1086             sage: c0 = sd[0][1]
1087             sage: c1 = sd[1][1]
1088             sage: c0.inner_product(c1) == 0
1089             True
1090             sage: c0.is_idempotent()
1091             True
1092             sage: c1.is_idempotent()
1093             True
1094             sage: c0 + c1 == J.one()
1095             True
1096             sage: l0*c0 + l1*c1 == x
1097             True
1098
1099         """
1100         P = self.parent()
1101         A = self.subalgebra_generated_by(orthonormalize_basis=True)
1102         result = []
1103         for (evalue, proj) in A(self).operator().spectral_decomposition():
1104             result.append( (evalue, proj(A.one()).superalgebra_element()) )
1105         return result
1106
1107     def subalgebra_generated_by(self, orthonormalize_basis=False):
1108         """
1109         Return the associative subalgebra of the parent EJA generated
1110         by this element.
1111
1112         Since our parent algebra is unital, we want "subalgebra" to mean
1113         "unital subalgebra" as well; thus the subalgebra that an element
1114         generates will itself be a Euclidean Jordan algebra after
1115         restricting the algebra operations appropriately. This is the
1116         subalgebra that Faraut and Korányi work with in section II.2, for
1117         example.
1118
1119         SETUP::
1120
1121             sage: from mjo.eja.eja_algebra import random_eja
1122
1123         TESTS:
1124
1125         This subalgebra, being composed of only powers, is associative::
1126
1127             sage: set_random_seed()
1128             sage: x0 = random_eja().random_element()
1129             sage: A = x0.subalgebra_generated_by()
1130             sage: x,y,z = A.random_elements(3)
1131             sage: (x*y)*z == x*(y*z)
1132             True
1133
1134         Squaring in the subalgebra should work the same as in
1135         the superalgebra::
1136
1137             sage: set_random_seed()
1138             sage: x = random_eja().random_element()
1139             sage: A = x.subalgebra_generated_by()
1140             sage: A(x^2) == A(x)*A(x)
1141             True
1142
1143         By definition, the subalgebra generated by the zero element is the
1144         one-dimensional algebra generated by the identity element::
1145
1146             sage: set_random_seed()
1147             sage: A = random_eja().zero().subalgebra_generated_by()
1148             sage: A.dimension()
1149             1
1150
1151         """
1152         return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
1153
1154
1155     def subalgebra_idempotent(self):
1156         """
1157         Find an idempotent in the associative subalgebra I generate
1158         using Proposition 2.3.5 in Baes.
1159
1160         SETUP::
1161
1162             sage: from mjo.eja.eja_algebra import random_eja
1163
1164         TESTS::
1165
1166             sage: set_random_seed()
1167             sage: J = random_eja()
1168             sage: x = J.random_element()
1169             sage: while x.is_nilpotent():
1170             ....:     x = J.random_element()
1171             sage: c = x.subalgebra_idempotent()
1172             sage: c^2 == c
1173             True
1174
1175         """
1176         if self.is_nilpotent():
1177             raise ValueError("this only works with non-nilpotent elements!")
1178
1179         J = self.subalgebra_generated_by()
1180         u = J(self)
1181
1182         # The image of the matrix of left-u^m-multiplication
1183         # will be minimal for some natural number s...
1184         s = 0
1185         minimal_dim = J.dimension()
1186         for i in xrange(1, minimal_dim):
1187             this_dim = (u**i).operator().matrix().image().dimension()
1188             if this_dim < minimal_dim:
1189                 minimal_dim = this_dim
1190                 s = i
1191
1192         # Now minimal_matrix should correspond to the smallest
1193         # non-zero subspace in Baes's (or really, Koecher's)
1194         # proposition.
1195         #
1196         # However, we need to restrict the matrix to work on the
1197         # subspace... or do we? Can't we just solve, knowing that
1198         # A(c) = u^(s+1) should have a solution in the big space,
1199         # too?
1200         #
1201         # Beware, solve_right() means that we're using COLUMN vectors.
1202         # Our FiniteDimensionalAlgebraElement superclass uses rows.
1203         u_next = u**(s+1)
1204         A = u_next.operator().matrix()
1205         c = J.from_vector(A.solve_right(u_next.to_vector()))
1206
1207         # Now c is the idempotent we want, but it still lives in the subalgebra.
1208         return c.superalgebra_element()
1209
1210
1211     def trace(self):
1212         """
1213         Return my trace, the sum of my eigenvalues.
1214
1215         SETUP::
1216
1217             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1218             ....:                                  RealCartesianProductEJA,
1219             ....:                                  random_eja)
1220
1221         EXAMPLES::
1222
1223             sage: J = JordanSpinEJA(3)
1224             sage: x = sum(J.gens())
1225             sage: x.trace()
1226             2
1227
1228         ::
1229
1230             sage: J = RealCartesianProductEJA(5)
1231             sage: J.one().trace()
1232             5
1233
1234         TESTS:
1235
1236         The trace of an element is a real number::
1237
1238             sage: set_random_seed()
1239             sage: J = random_eja()
1240             sage: J.random_element().trace() in RLF
1241             True
1242
1243         """
1244         P = self.parent()
1245         r = P.rank()
1246         p = P._charpoly_coeff(r-1)
1247         # The _charpoly_coeff function already adds the factor of
1248         # -1 to ensure that _charpoly_coeff(r-1) is really what
1249         # appears in front of t^{r-1} in the charpoly. However,
1250         # we want the negative of THAT for the trace.
1251         return -p(*self.to_vector())
1252
1253
1254     def trace_inner_product(self, other):
1255         """
1256         Return the trace inner product of myself and ``other``.
1257
1258         SETUP::
1259
1260             sage: from mjo.eja.eja_algebra import random_eja
1261
1262         TESTS:
1263
1264         The trace inner product is commutative, bilinear, and associative::
1265
1266             sage: set_random_seed()
1267             sage: J = random_eja()
1268             sage: x,y,z = J.random_elements(3)
1269             sage: # commutative
1270             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1271             True
1272             sage: # bilinear
1273             sage: a = J.base_ring().random_element();
1274             sage: actual = (a*(x+z)).trace_inner_product(y)
1275             sage: expected = ( a*x.trace_inner_product(y) +
1276             ....:              a*z.trace_inner_product(y) )
1277             sage: actual == expected
1278             True
1279             sage: actual = x.trace_inner_product(a*(y+z))
1280             sage: expected = ( a*x.trace_inner_product(y) +
1281             ....:              a*x.trace_inner_product(z) )
1282             sage: actual == expected
1283             True
1284             sage: # associative
1285             sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1286             True
1287
1288         """
1289         if not other in self.parent():
1290             raise TypeError("'other' must live in the same algebra")
1291
1292         return (self*other).trace()
1293
1294
1295     def trace_norm(self):
1296         """
1297         The norm of this element with respect to :meth:`trace_inner_product`.
1298
1299         SETUP::
1300
1301             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
1302             ....:                                  RealCartesianProductEJA)
1303
1304         EXAMPLES::
1305
1306             sage: J = RealCartesianProductEJA(2)
1307             sage: x = sum(J.gens())
1308             sage: x.trace_norm()
1309             sqrt(2)
1310
1311         ::
1312
1313             sage: J = JordanSpinEJA(4)
1314             sage: x = sum(J.gens())
1315             sage: x.trace_norm()
1316             2*sqrt(2)
1317
1318         """
1319         return self.trace_inner_product(self).sqrt()