mjo/eja/euclidean_jordan_algebra.py

   1 """
   2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
   3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
   4 are used in optimization, and have some additional nice methods beyond
   5 what can be supported in a general Jordan Algebra.
   6 """
   7
   8 from sage.categories.magmatic_algebras import MagmaticAlgebras
   9 from sage.structure.element import is_Matrix
  10 from sage.structure.category_object import normalize_names
  11
  12 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra
  13 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element import FiniteDimensionalAlgebraElement
  14
  15 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
  16     @staticmethod
  17     def __classcall_private__(cls,
  18                               field,
  19                               mult_table,
  20                               names='e',
  21                               assume_associative=False,
  22                               category=None,
  23                               rank=None,
  24                               natural_basis=None,
  25                               inner_product=None):
  26         n = len(mult_table)
  27         mult_table = [b.base_extend(field) for b in mult_table]
  28         for b in mult_table:
  29             b.set_immutable()
  30             if not (is_Matrix(b) and b.dimensions() == (n, n)):
  31                 raise ValueError("input is not a multiplication table")
  32         mult_table = tuple(mult_table)
  33
  34         cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
  35         cat.or_subcategory(category)
  36         if assume_associative:
  37             cat = cat.Associative()
  38
  39         names = normalize_names(n, names)
  40
  41         fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls)
  42         return fda.__classcall__(cls,
  43                                  field,
  44                                  mult_table,
  45                                  assume_associative=assume_associative,
  46                                  names=names,
  47                                  category=cat,
  48                                  rank=rank,
  49                                  natural_basis=natural_basis,
  50                                  inner_product=inner_product)
  51
  52
  53     def __init__(self, field,
  54                  mult_table,
  55                  names='e',
  56                  assume_associative=False,
  57                  category=None,
  58                  rank=None,
  59                  natural_basis=None,
  60                  inner_product=None):
  61         """
  62         EXAMPLES:
  63
  64         By definition, Jordan multiplication commutes::
  65
  66             sage: set_random_seed()
  67             sage: J = random_eja()
  68             sage: x = J.random_element()
  69             sage: y = J.random_element()
  70             sage: x*y == y*x
  71             True
  72
  73         """
  74         self._rank = rank
  75         self._natural_basis = natural_basis
  76         self._inner_product = inner_product
  77         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
  78         fda.__init__(field,
  79                      mult_table,
  80                      names=names,
  81                      category=category)
  82
  83
  84     def _repr_(self):
  85         """
  86         Return a string representation of ``self``.
  87         """
  88         fmt = "Euclidean Jordan algebra of degree {} over {}"
  89         return fmt.format(self.degree(), self.base_ring())
  90
  91
  92     def inner_product(self, x, y):
  93         """
  94         The inner product associated with this Euclidean Jordan algebra.
  95
  96         Will default to the trace inner product if nothing else.
  97         """
  98         if (not x in self) or (not y in self):
  99             raise ArgumentError("arguments must live in this algebra")
 100         if self._inner_product is None:
 101             return x.trace_inner_product(y)
 102         else:
 103             return self._inner_product(x,y)
 104
 105
 106     def natural_basis(self):
 107         """
 108         Return a more-natural representation of this algebra's basis.
 109
 110         Every finite-dimensional Euclidean Jordan Algebra is a direct
 111         sum of five simple algebras, four of which comprise Hermitian
 112         matrices. This method returns the original "natural" basis
 113         for our underlying vector space. (Typically, the natural basis
 114         is used to construct the multiplication table in the first place.)
 115
 116         Note that this will always return a matrix. The standard basis
 117         in `R^n` will be returned as `n`-by-`1` column matrices.
 118
 119         EXAMPLES::
 120
 121             sage: J = RealSymmetricSimpleEJA(2)
 122             sage: J.basis()
 123             Family (e0, e1, e2)
 124             sage: J.natural_basis()
 125             (
 126             [1 0]  [0 1]  [0 0]
 127             [0 0], [1 0], [0 1]
 128             )
 129
 130         ::
 131
 132             sage: J = JordanSpinSimpleEJA(2)
 133             sage: J.basis()
 134             Family (e0, e1)
 135             sage: J.natural_basis()
 136             (
 137             [1]  [0]
 138             [0], [1]
 139             )
 140
 141         """
 142         if self._natural_basis is None:
 143             return tuple( b.vector().column() for b in self.basis() )
 144         else:
 145             return self._natural_basis
 146
 147
 148     def rank(self):
 149         """
 150         Return the rank of this EJA.
 151         """
 152         if self._rank is None:
 153             raise ValueError("no rank specified at genesis")
 154         else:
 155             return self._rank
 156
 157
 158     class Element(FiniteDimensionalAlgebraElement):
 159         """
 160         An element of a Euclidean Jordan algebra.
 161         """
 162
 163         def __pow__(self, n):
 164             """
 165             Return ``self`` raised to the power ``n``.
 166
 167             Jordan algebras are always power-associative; see for
 168             example Faraut and Koranyi, Proposition II.1.2 (ii).
 169
 170             .. WARNING:
 171
 172                 We have to override this because our superclass uses row vectors
 173                 instead of column vectors! We, on the other hand, assume column
 174                 vectors everywhere.
 175
 176             EXAMPLES::
 177
 178                 sage: set_random_seed()
 179                 sage: x = random_eja().random_element()
 180                 sage: x.operator_matrix()*x.vector() == (x^2).vector()
 181                 True
 182
 183             A few examples of power-associativity::
 184
 185                 sage: set_random_seed()
 186                 sage: x = random_eja().random_element()
 187                 sage: x*(x*x)*(x*x) == x^5
 188                 True
 189                 sage: (x*x)*(x*x*x) == x^5
 190                 True
 191
 192             We also know that powers operator-commute (Koecher, Chapter
 193             III, Corollary 1)::
 194
 195                 sage: set_random_seed()
 196                 sage: x = random_eja().random_element()
 197                 sage: m = ZZ.random_element(0,10)
 198                 sage: n = ZZ.random_element(0,10)
 199                 sage: Lxm = (x^m).operator_matrix()
 200                 sage: Lxn = (x^n).operator_matrix()
 201                 sage: Lxm*Lxn == Lxn*Lxm
 202                 True
 203
 204             """
 205             A = self.parent()
 206             if n == 0:
 207                 return A.one()
 208             elif n == 1:
 209                 return self
 210             else:
 211                 return A( (self.operator_matrix()**(n-1))*self.vector() )
 212
 213
 214         def characteristic_polynomial(self):
 215             """
 216             Return my characteristic polynomial (if I'm a regular
 217             element).
 218
 219             Eventually this should be implemented in terms of the parent
 220             algebra's characteristic polynomial that works for ALL
 221             elements.
 222             """
 223             if self.is_regular():
 224                 return self.minimal_polynomial()
 225             else:
 226                 raise NotImplementedError('irregular element')
 227
 228
 229         def inner_product(self, other):
 230             """
 231             Return the parent algebra's inner product of myself and ``other``.
 232
 233             EXAMPLES:
 234
 235             The inner product in the Jordan spin algebra is the usual
 236             inner product on `R^n` (this example only works because the
 237             basis for the Jordan algebra is the standard basis in `R^n`)::
 238
 239                 sage: J = JordanSpinSimpleEJA(3)
 240                 sage: x = vector(QQ,[1,2,3])
 241                 sage: y = vector(QQ,[4,5,6])
 242                 sage: x.inner_product(y)
 243                 32
 244                 sage: J(x).inner_product(J(y))
 245                 32
 246
 247             TESTS:
 248
 249             Ensure that we can always compute an inner product, and that
 250             it gives us back a real number::
 251
 252                 sage: set_random_seed()
 253                 sage: J = random_eja()
 254                 sage: x = J.random_element()
 255                 sage: y = J.random_element()
 256                 sage: x.inner_product(y) in RR
 257                 True
 258
 259             """
 260             P = self.parent()
 261             if not other in P:
 262                 raise ArgumentError("'other' must live in the same algebra")
 263
 264             return P.inner_product(self, other)
 265
 266
 267         def operator_commutes_with(self, other):
 268             """
 269             Return whether or not this element operator-commutes
 270             with ``other``.
 271
 272             EXAMPLES:
 273
 274             The definition of a Jordan algebra says that any element
 275             operator-commutes with its square::
 276
 277                 sage: set_random_seed()
 278                 sage: x = random_eja().random_element()
 279                 sage: x.operator_commutes_with(x^2)
 280                 True
 281
 282             TESTS:
 283
 284             Test Lemma 1 from Chapter III of Koecher::
 285
 286                 sage: set_random_seed()
 287                 sage: J = random_eja()
 288                 sage: u = J.random_element()
 289                 sage: v = J.random_element()
 290                 sage: lhs = u.operator_commutes_with(u*v)
 291                 sage: rhs = v.operator_commutes_with(u^2)
 292                 sage: lhs == rhs
 293                 True
 294
 295             """
 296             if not other in self.parent():
 297                 raise ArgumentError("'other' must live in the same algebra")
 298
 299             A = self.operator_matrix()
 300             B = other.operator_matrix()
 301             return (A*B == B*A)
 302
 303
 304         def det(self):
 305             """
 306             Return my determinant, the product of my eigenvalues.
 307
 308             EXAMPLES::
 309
 310                 sage: J = JordanSpinSimpleEJA(2)
 311                 sage: e0,e1 = J.gens()
 312                 sage: x = e0 + e1
 313                 sage: x.det()
 314                 0
 315                 sage: J = JordanSpinSimpleEJA(3)
 316                 sage: e0,e1,e2 = J.gens()
 317                 sage: x = e0 + e1 + e2
 318                 sage: x.det()
 319                 -1
 320
 321             """
 322             cs = self.characteristic_polynomial().coefficients(sparse=False)
 323             r = len(cs) - 1
 324             if r >= 0:
 325                 return cs[0] * (-1)**r
 326             else:
 327                 raise ValueError('charpoly had no coefficients')
 328
 329
 330         def inverse(self):
 331             """
 332             Return the Jordan-multiplicative inverse of this element.
 333
 334             We can't use the superclass method because it relies on the
 335             algebra being associative.
 336
 337             EXAMPLES:
 338
 339             The inverse in the spin factor algebra is given in Alizadeh's
 340             Example 11.11::
 341
 342                 sage: set_random_seed()
 343                 sage: n = ZZ.random_element(1,10)
 344                 sage: J = JordanSpinSimpleEJA(n)
 345                 sage: x = J.random_element()
 346                 sage: while x.is_zero():
 347                 ....:     x = J.random_element()
 348                 sage: x_vec = x.vector()
 349                 sage: x0 = x_vec[0]
 350                 sage: x_bar = x_vec[1:]
 351                 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
 352                 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
 353                 sage: x_inverse = coeff*inv_vec
 354                 sage: x.inverse() == J(x_inverse)
 355                 True
 356
 357             TESTS:
 358
 359             The identity element is its own inverse::
 360
 361                 sage: set_random_seed()
 362                 sage: J = random_eja()
 363                 sage: J.one().inverse() == J.one()
 364                 True
 365
 366             If an element has an inverse, it acts like one. TODO: this
 367             can be a lot less ugly once ``is_invertible`` doesn't crash
 368             on irregular elements::
 369
 370                 sage: set_random_seed()
 371                 sage: J = random_eja()
 372                 sage: x = J.random_element()
 373                 sage: try:
 374                 ....:     x.inverse()*x == J.one()
 375                 ....: except:
 376                 ....:     True
 377                 True
 378
 379             """
 380             if self.parent().is_associative():
 381                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 382                 return elt.inverse()
 383
 384             # TODO: we can do better once the call to is_invertible()
 385             # doesn't crash on irregular elements.
 386             #if not self.is_invertible():
 387             #    raise ArgumentError('element is not invertible')
 388
 389             # We do this a little different than the usual recursive
 390             # call to a finite-dimensional algebra element, because we
 391             # wind up with an inverse that lives in the subalgebra and
 392             # we need information about the parent to convert it back.
 393             V = self.span_of_powers()
 394             assoc_subalg = self.subalgebra_generated_by()
 395             # Mis-design warning: the basis used for span_of_powers()
 396             # and subalgebra_generated_by() must be the same, and in
 397             # the same order!
 398             elt = assoc_subalg(V.coordinates(self.vector()))
 399
 400             # This will be in the subalgebra's coordinates...
 401             fda_elt = FiniteDimensionalAlgebraElement(assoc_subalg, elt)
 402             subalg_inverse = fda_elt.inverse()
 403
 404             # So we have to convert back...
 405             basis = [ self.parent(v) for v in V.basis() ]
 406             pairs = zip(subalg_inverse.vector(), basis)
 407             return self.parent().linear_combination(pairs)
 408
 409
 410         def is_invertible(self):
 411             """
 412             Return whether or not this element is invertible.
 413
 414             We can't use the superclass method because it relies on
 415             the algebra being associative.
 416             """
 417             return not self.det().is_zero()
 418
 419
 420         def is_nilpotent(self):
 421             """
 422             Return whether or not some power of this element is zero.
 423
 424             The superclass method won't work unless we're in an
 425             associative algebra, and we aren't. However, we generate
 426             an assocoative subalgebra and we're nilpotent there if and
 427             only if we're nilpotent here (probably).
 428
 429             TESTS:
 430
 431             The identity element is never nilpotent::
 432
 433                 sage: set_random_seed()
 434                 sage: random_eja().one().is_nilpotent()
 435                 False
 436
 437             The additive identity is always nilpotent::
 438
 439                 sage: set_random_seed()
 440                 sage: random_eja().zero().is_nilpotent()
 441                 True
 442
 443             """
 444             # The element we're going to call "is_nilpotent()" on.
 445             # Either myself, interpreted as an element of a finite-
 446             # dimensional algebra, or an element of an associative
 447             # subalgebra.
 448             elt = None
 449
 450             if self.parent().is_associative():
 451                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 452             else:
 453                 V = self.span_of_powers()
 454                 assoc_subalg = self.subalgebra_generated_by()
 455                 # Mis-design warning: the basis used for span_of_powers()
 456                 # and subalgebra_generated_by() must be the same, and in
 457                 # the same order!
 458                 elt = assoc_subalg(V.coordinates(self.vector()))
 459
 460             # Recursive call, but should work since elt lives in an
 461             # associative algebra.
 462             return elt.is_nilpotent()
 463
 464
 465         def is_regular(self):
 466             """
 467             Return whether or not this is a regular element.
 468
 469             EXAMPLES:
 470
 471             The identity element always has degree one, but any element
 472             linearly-independent from it is regular::
 473
 474                 sage: J = JordanSpinSimpleEJA(5)
 475                 sage: J.one().is_regular()
 476                 False
 477                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 478                 sage: for x in J.gens():
 479                 ....:     (J.one() + x).is_regular()
 480                 False
 481                 True
 482                 True
 483                 True
 484                 True
 485
 486             """
 487             return self.degree() == self.parent().rank()
 488
 489
 490         def degree(self):
 491             """
 492             Compute the degree of this element the straightforward way
 493             according to the definition; by appending powers to a list
 494             and figuring out its dimension (that is, whether or not
 495             they're linearly dependent).
 496
 497             EXAMPLES::
 498
 499                 sage: J = JordanSpinSimpleEJA(4)
 500                 sage: J.one().degree()
 501                 1
 502                 sage: e0,e1,e2,e3 = J.gens()
 503                 sage: (e0 - e1).degree()
 504                 2
 505
 506             In the spin factor algebra (of rank two), all elements that
 507             aren't multiples of the identity are regular::
 508
 509                 sage: set_random_seed()
 510                 sage: n = ZZ.random_element(1,10)
 511                 sage: J = JordanSpinSimpleEJA(n)
 512                 sage: x = J.random_element()
 513                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 514                 True
 515
 516             """
 517             return self.span_of_powers().dimension()
 518
 519
 520         def minimal_polynomial(self):
 521             """
 522             EXAMPLES::
 523
 524                 sage: set_random_seed()
 525                 sage: x = random_eja().random_element()
 526                 sage: x.degree() == x.minimal_polynomial().degree()
 527                 True
 528
 529             ::
 530
 531                 sage: set_random_seed()
 532                 sage: x = random_eja().random_element()
 533                 sage: x.degree() == x.minimal_polynomial().degree()
 534                 True
 535
 536             The minimal polynomial and the characteristic polynomial coincide
 537             and are known (see Alizadeh, Example 11.11) for all elements of
 538             the spin factor algebra that aren't scalar multiples of the
 539             identity::
 540
 541                 sage: set_random_seed()
 542                 sage: n = ZZ.random_element(2,10)
 543                 sage: J = JordanSpinSimpleEJA(n)
 544                 sage: y = J.random_element()
 545                 sage: while y == y.coefficient(0)*J.one():
 546                 ....:     y = J.random_element()
 547                 sage: y0 = y.vector()[0]
 548                 sage: y_bar = y.vector()[1:]
 549                 sage: actual = y.minimal_polynomial()
 550                 sage: x = SR.symbol('x', domain='real')
 551                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 552                 sage: bool(actual == expected)
 553                 True
 554
 555             """
 556             # The element we're going to call "minimal_polynomial()" on.
 557             # Either myself, interpreted as an element of a finite-
 558             # dimensional algebra, or an element of an associative
 559             # subalgebra.
 560             elt = None
 561
 562             if self.parent().is_associative():
 563                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 564             else:
 565                 V = self.span_of_powers()
 566                 assoc_subalg = self.subalgebra_generated_by()
 567                 # Mis-design warning: the basis used for span_of_powers()
 568                 # and subalgebra_generated_by() must be the same, and in
 569                 # the same order!
 570                 elt = assoc_subalg(V.coordinates(self.vector()))
 571
 572             # Recursive call, but should work since elt lives in an
 573             # associative algebra.
 574             return elt.minimal_polynomial()
 575
 576
 577         def natural_representation(self):
 578             """
 579             Return a more-natural representation of this element.
 580
 581             Every finite-dimensional Euclidean Jordan Algebra is a
 582             direct sum of five simple algebras, four of which comprise
 583             Hermitian matrices. This method returns the original
 584             "natural" representation of this element as a Hermitian
 585             matrix, if it has one. If not, you get the usual representation.
 586
 587             EXAMPLES::
 588
 589                 sage: J = ComplexHermitianSimpleEJA(3)
 590                 sage: J.one()
 591                 e0 + e5 + e8
 592                 sage: J.one().natural_representation()
 593                 [1 0 0 0 0 0]
 594                 [0 1 0 0 0 0]
 595                 [0 0 1 0 0 0]
 596                 [0 0 0 1 0 0]
 597                 [0 0 0 0 1 0]
 598                 [0 0 0 0 0 1]
 599
 600             """
 601             B = self.parent().natural_basis()
 602             W = B[0].matrix_space()
 603             return W.linear_combination(zip(self.vector(), B))
 604
 605
 606         def operator_matrix(self):
 607             """
 608             Return the matrix that represents left- (or right-)
 609             multiplication by this element in the parent algebra.
 610
 611             We have to override this because the superclass method
 612             returns a matrix that acts on row vectors (that is, on
 613             the right).
 614
 615             EXAMPLES:
 616
 617             Test the first polarization identity from my notes, Koecher Chapter
 618             III, or from Baes (2.3)::
 619
 620                 sage: set_random_seed()
 621                 sage: J = random_eja()
 622                 sage: x = J.random_element()
 623                 sage: y = J.random_element()
 624                 sage: Lx = x.operator_matrix()
 625                 sage: Ly = y.operator_matrix()
 626                 sage: Lxx = (x*x).operator_matrix()
 627                 sage: Lxy = (x*y).operator_matrix()
 628                 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
 629                 True
 630
 631             Test the second polarization identity from my notes or from
 632             Baes (2.4)::
 633
 634                 sage: set_random_seed()
 635                 sage: J = random_eja()
 636                 sage: x = J.random_element()
 637                 sage: y = J.random_element()
 638                 sage: z = J.random_element()
 639                 sage: Lx = x.operator_matrix()
 640                 sage: Ly = y.operator_matrix()
 641                 sage: Lz = z.operator_matrix()
 642                 sage: Lzy = (z*y).operator_matrix()
 643                 sage: Lxy = (x*y).operator_matrix()
 644                 sage: Lxz = (x*z).operator_matrix()
 645                 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
 646                 True
 647
 648             Test the third polarization identity from my notes or from
 649             Baes (2.5)::
 650
 651                 sage: set_random_seed()
 652                 sage: J = random_eja()
 653                 sage: u = J.random_element()
 654                 sage: y = J.random_element()
 655                 sage: z = J.random_element()
 656                 sage: Lu = u.operator_matrix()
 657                 sage: Ly = y.operator_matrix()
 658                 sage: Lz = z.operator_matrix()
 659                 sage: Lzy = (z*y).operator_matrix()
 660                 sage: Luy = (u*y).operator_matrix()
 661                 sage: Luz = (u*z).operator_matrix()
 662                 sage: Luyz = (u*(y*z)).operator_matrix()
 663                 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
 664                 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
 665                 sage: bool(lhs == rhs)
 666                 True
 667
 668             """
 669             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 670             return fda_elt.matrix().transpose()
 671
 672
 673         def quadratic_representation(self, other=None):
 674             """
 675             Return the quadratic representation of this element.
 676
 677             EXAMPLES:
 678
 679             The explicit form in the spin factor algebra is given by
 680             Alizadeh's Example 11.12::
 681
 682                 sage: set_random_seed()
 683                 sage: n = ZZ.random_element(1,10)
 684                 sage: J = JordanSpinSimpleEJA(n)
 685                 sage: x = J.random_element()
 686                 sage: x_vec = x.vector()
 687                 sage: x0 = x_vec[0]
 688                 sage: x_bar = x_vec[1:]
 689                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 690                 sage: B = 2*x0*x_bar.row()
 691                 sage: C = 2*x0*x_bar.column()
 692                 sage: D = identity_matrix(QQ, n-1)
 693                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 694                 sage: D = D + 2*x_bar.tensor_product(x_bar)
 695                 sage: Q = block_matrix(2,2,[A,B,C,D])
 696                 sage: Q == x.quadratic_representation()
 697                 True
 698
 699             Test all of the properties from Theorem 11.2 in Alizadeh::
 700
 701                 sage: set_random_seed()
 702                 sage: J = random_eja()
 703                 sage: x = J.random_element()
 704                 sage: y = J.random_element()
 705
 706             Property 1:
 707
 708                 sage: actual = x.quadratic_representation(y)
 709                 sage: expected = ( (x+y).quadratic_representation()
 710                 ....:              -x.quadratic_representation()
 711                 ....:              -y.quadratic_representation() ) / 2
 712                 sage: actual == expected
 713                 True
 714
 715             Property 2:
 716
 717                 sage: alpha = QQ.random_element()
 718                 sage: actual = (alpha*x).quadratic_representation()
 719                 sage: expected = (alpha^2)*x.quadratic_representation()
 720                 sage: actual == expected
 721                 True
 722
 723             Property 5:
 724
 725                 sage: Qy = y.quadratic_representation()
 726                 sage: actual = J(Qy*x.vector()).quadratic_representation()
 727                 sage: expected = Qy*x.quadratic_representation()*Qy
 728                 sage: actual == expected
 729                 True
 730
 731             Property 6:
 732
 733                 sage: k = ZZ.random_element(1,10)
 734                 sage: actual = (x^k).quadratic_representation()
 735                 sage: expected = (x.quadratic_representation())^k
 736                 sage: actual == expected
 737                 True
 738
 739             """
 740             if other is None:
 741                 other=self
 742             elif not other in self.parent():
 743                 raise ArgumentError("'other' must live in the same algebra")
 744
 745             L = self.operator_matrix()
 746             M = other.operator_matrix()
 747             return ( L*M + M*L - (self*other).operator_matrix() )
 748
 749
 750         def span_of_powers(self):
 751             """
 752             Return the vector space spanned by successive powers of
 753             this element.
 754             """
 755             # The dimension of the subalgebra can't be greater than
 756             # the big algebra, so just put everything into a list
 757             # and let span() get rid of the excess.
 758             V = self.vector().parent()
 759             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 760
 761
 762         def subalgebra_generated_by(self):
 763             """
 764             Return the associative subalgebra of the parent EJA generated
 765             by this element.
 766
 767             TESTS::
 768
 769                 sage: set_random_seed()
 770                 sage: x = random_eja().random_element()
 771                 sage: x.subalgebra_generated_by().is_associative()
 772                 True
 773
 774             Squaring in the subalgebra should be the same thing as
 775             squaring in the superalgebra::
 776
 777                 sage: set_random_seed()
 778                 sage: x = random_eja().random_element()
 779                 sage: u = x.subalgebra_generated_by().random_element()
 780                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
 781                 True
 782
 783             """
 784             # First get the subspace spanned by the powers of myself...
 785             V = self.span_of_powers()
 786             F = self.base_ring()
 787
 788             # Now figure out the entries of the right-multiplication
 789             # matrix for the successive basis elements b0, b1,... of
 790             # that subspace.
 791             mats = []
 792             for b_right in V.basis():
 793                 eja_b_right = self.parent()(b_right)
 794                 b_right_rows = []
 795                 # The first row of the right-multiplication matrix by
 796                 # b1 is what we get if we apply that matrix to b1. The
 797                 # second row of the right multiplication matrix by b1
 798                 # is what we get when we apply that matrix to b2...
 799                 #
 800                 # IMPORTANT: this assumes that all vectors are COLUMN
 801                 # vectors, unlike our superclass (which uses row vectors).
 802                 for b_left in V.basis():
 803                     eja_b_left = self.parent()(b_left)
 804                     # Multiply in the original EJA, but then get the
 805                     # coordinates from the subalgebra in terms of its
 806                     # basis.
 807                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 808                     b_right_rows.append(this_row)
 809                 b_right_matrix = matrix(F, b_right_rows)
 810                 mats.append(b_right_matrix)
 811
 812             # It's an algebra of polynomials in one element, and EJAs
 813             # are power-associative.
 814             #
 815             # TODO: choose generator names intelligently.
 816             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 817
 818
 819         def subalgebra_idempotent(self):
 820             """
 821             Find an idempotent in the associative subalgebra I generate
 822             using Proposition 2.3.5 in Baes.
 823
 824             TESTS::
 825
 826                 sage: set_random_seed()
 827                 sage: J = eja_rn(5)
 828                 sage: c = J.random_element().subalgebra_idempotent()
 829                 sage: c^2 == c
 830                 True
 831                 sage: J = JordanSpinSimpleEJA(5)
 832                 sage: c = J.random_element().subalgebra_idempotent()
 833                 sage: c^2 == c
 834                 True
 835
 836             """
 837             if self.is_nilpotent():
 838                 raise ValueError("this only works with non-nilpotent elements!")
 839
 840             V = self.span_of_powers()
 841             J = self.subalgebra_generated_by()
 842             # Mis-design warning: the basis used for span_of_powers()
 843             # and subalgebra_generated_by() must be the same, and in
 844             # the same order!
 845             u = J(V.coordinates(self.vector()))
 846
 847             # The image of the matrix of left-u^m-multiplication
 848             # will be minimal for some natural number s...
 849             s = 0
 850             minimal_dim = V.dimension()
 851             for i in xrange(1, V.dimension()):
 852                 this_dim = (u**i).operator_matrix().image().dimension()
 853                 if this_dim < minimal_dim:
 854                     minimal_dim = this_dim
 855                     s = i
 856
 857             # Now minimal_matrix should correspond to the smallest
 858             # non-zero subspace in Baes's (or really, Koecher's)
 859             # proposition.
 860             #
 861             # However, we need to restrict the matrix to work on the
 862             # subspace... or do we? Can't we just solve, knowing that
 863             # A(c) = u^(s+1) should have a solution in the big space,
 864             # too?
 865             #
 866             # Beware, solve_right() means that we're using COLUMN vectors.
 867             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 868             u_next = u**(s+1)
 869             A = u_next.operator_matrix()
 870             c_coordinates = A.solve_right(u_next.vector())
 871
 872             # Now c_coordinates is the idempotent we want, but it's in
 873             # the coordinate system of the subalgebra.
 874             #
 875             # We need the basis for J, but as elements of the parent algebra.
 876             #
 877             basis = [self.parent(v) for v in V.basis()]
 878             return self.parent().linear_combination(zip(c_coordinates, basis))
 879
 880
 881         def trace(self):
 882             """
 883             Return my trace, the sum of my eigenvalues.
 884
 885             EXAMPLES::
 886
 887                 sage: J = JordanSpinSimpleEJA(3)
 888                 sage: e0,e1,e2 = J.gens()
 889                 sage: x = e0 + e1 + e2
 890                 sage: x.trace()
 891                 2
 892
 893             """
 894             cs = self.characteristic_polynomial().coefficients(sparse=False)
 895             if len(cs) >= 2:
 896                 return -1*cs[-2]
 897             else:
 898                 raise ValueError('charpoly had fewer than 2 coefficients')
 899
 900
 901         def trace_inner_product(self, other):
 902             """
 903             Return the trace inner product of myself and ``other``.
 904             """
 905             if not other in self.parent():
 906                 raise ArgumentError("'other' must live in the same algebra")
 907
 908             return (self*other).trace()
 909
 910
 911 def eja_rn(dimension, field=QQ):
 912     """
 913     Return the Euclidean Jordan Algebra corresponding to the set
 914     `R^n` under the Hadamard product.
 915
 916     EXAMPLES:
 917
 918     This multiplication table can be verified by hand::
 919
 920         sage: J = eja_rn(3)
 921         sage: e0,e1,e2 = J.gens()
 922         sage: e0*e0
 923         e0
 924         sage: e0*e1
 925         0
 926         sage: e0*e2
 927         0
 928         sage: e1*e1
 929         e1
 930         sage: e1*e2
 931         0
 932         sage: e2*e2
 933         e2
 934
 935     """
 936     # The FiniteDimensionalAlgebra constructor takes a list of
 937     # matrices, the ith representing right multiplication by the ith
 938     # basis element in the vector space. So if e_1 = (1,0,0), then
 939     # right (Hadamard) multiplication of x by e_1 picks out the first
 940     # component of x; and likewise for the ith basis element e_i.
 941     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
 942            for i in xrange(dimension) ]
 943
 944     # The usual inner product on R^n.
 945     ip = lambda x, y: x.vector().inner_product(y.vector())
 946
 947     return FiniteDimensionalEuclideanJordanAlgebra(field,
 948                                                    Qs,
 949                                                    rank=dimension,
 950                                                    inner_product=ip)
 951
 952
 953
 954 def random_eja():
 955     """
 956     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 957
 958     ALGORITHM:
 959
 960     For now, we choose a random natural number ``n`` (greater than zero)
 961     and then give you back one of the following:
 962
 963       * The cartesian product of the rational numbers ``n`` times; this is
 964         ``QQ^n`` with the Hadamard product.
 965
 966       * The Jordan spin algebra on ``QQ^n``.
 967
 968       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
 969         product.
 970
 971     Later this might be extended to return Cartesian products of the
 972     EJAs above.
 973
 974     TESTS::
 975
 976         sage: random_eja()
 977         Euclidean Jordan algebra of degree...
 978
 979     """
 980     n = ZZ.random_element(1,5)
 981     constructor = choice([eja_rn,
 982                           JordanSpinSimpleEJA,
 983                           RealSymmetricSimpleEJA,
 984                           ComplexHermitianSimpleEJA])
 985     return constructor(n, field=QQ)
 986
 987
 988
 989 def _real_symmetric_basis(n, field=QQ):
 990     """
 991     Return a basis for the space of real symmetric n-by-n matrices.
 992     """
 993     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
 994     # coordinates.
 995     S = []
 996     for i in xrange(n):
 997         for j in xrange(i+1):
 998             Eij = matrix(field, n, lambda k,l: k==i and l==j)
 999             if i == j:
1000                 Sij = Eij
1001             else:
1002                 # Beware, orthogonal but not normalized!
1003                 Sij = Eij + Eij.transpose()
1004             S.append(Sij)
1005     return tuple(S)
1006
1007
1008 def _complex_hermitian_basis(n, field=QQ):
1009     """
1010     Returns a basis for the space of complex Hermitian n-by-n matrices.
1011
1012     TESTS::
1013
1014         sage: set_random_seed()
1015         sage: n = ZZ.random_element(1,5)
1016         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1017         True
1018
1019     """
1020     F = QuadraticField(-1, 'I')
1021     I = F.gen()
1022
1023     # This is like the symmetric case, but we need to be careful:
1024     #
1025     #   * We want conjugate-symmetry, not just symmetry.
1026     #   * The diagonal will (as a result) be real.
1027     #
1028     S = []
1029     for i in xrange(n):
1030         for j in xrange(i+1):
1031             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1032             if i == j:
1033                 Sij = _embed_complex_matrix(Eij)
1034                 S.append(Sij)
1035             else:
1036                 # Beware, orthogonal but not normalized! The second one
1037                 # has a minus because it's conjugated.
1038                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1039                 S.append(Sij_real)
1040                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1041                 S.append(Sij_imag)
1042     return tuple(S)
1043
1044
1045 def _multiplication_table_from_matrix_basis(basis):
1046     """
1047     At least three of the five simple Euclidean Jordan algebras have the
1048     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1049     multiplication on the right is matrix multiplication. Given a basis
1050     for the underlying matrix space, this function returns a
1051     multiplication table (obtained by looping through the basis
1052     elements) for an algebra of those matrices. A reordered copy
1053     of the basis is also returned to work around the fact that
1054     the ``span()`` in this function will change the order of the basis
1055     from what we think it is, to... something else.
1056     """
1057     # In S^2, for example, we nominally have four coordinates even
1058     # though the space is of dimension three only. The vector space V
1059     # is supposed to hold the entire long vector, and the subspace W
1060     # of V will be spanned by the vectors that arise from symmetric
1061     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1062     field = basis[0].base_ring()
1063     dimension = basis[0].nrows()
1064
1065     def mat2vec(m):
1066         return vector(field, m.list())
1067
1068     def vec2mat(v):
1069         return matrix(field, dimension, v.list())
1070
1071     V = VectorSpace(field, dimension**2)
1072     W = V.span( mat2vec(s) for s in basis )
1073
1074     # Taking the span above reorders our basis (thanks, jerk!) so we
1075     # need to put our "matrix basis" in the same order as the
1076     # (reordered) vector basis.
1077     S = tuple( vec2mat(b) for b in W.basis() )
1078
1079     Qs = []
1080     for s in S:
1081         # Brute force the multiplication-by-s matrix by looping
1082         # through all elements of the basis and doing the computation
1083         # to find out what the corresponding row should be. BEWARE:
1084         # these multiplication tables won't be symmetric! It therefore
1085         # becomes REALLY IMPORTANT that the underlying algebra
1086         # constructor uses ROW vectors and not COLUMN vectors. That's
1087         # why we're computing rows here and not columns.
1088         Q_rows = []
1089         for t in S:
1090             this_row = mat2vec((s*t + t*s)/2)
1091             Q_rows.append(W.coordinates(this_row))
1092         Q = matrix(field, W.dimension(), Q_rows)
1093         Qs.append(Q)
1094
1095     return (Qs, S)
1096
1097
1098 def _embed_complex_matrix(M):
1099     """
1100     Embed the n-by-n complex matrix ``M`` into the space of real
1101     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1102     bi` to the block matrix ``[[a,b],[-b,a]]``.
1103
1104     EXAMPLES::
1105
1106         sage: F = QuadraticField(-1,'i')
1107         sage: x1 = F(4 - 2*i)
1108         sage: x2 = F(1 + 2*i)
1109         sage: x3 = F(-i)
1110         sage: x4 = F(6)
1111         sage: M = matrix(F,2,[x1,x2,x3,x4])
1112         sage: _embed_complex_matrix(M)
1113         [ 4  2| 1 -2]
1114         [-2  4| 2  1]
1115         [-----+-----]
1116         [ 0  1| 6  0]
1117         [-1  0| 0  6]
1118
1119     """
1120     n = M.nrows()
1121     if M.ncols() != n:
1122         raise ArgumentError("the matrix 'M' must be square")
1123     field = M.base_ring()
1124     blocks = []
1125     for z in M.list():
1126         a = z.real()
1127         b = z.imag()
1128         blocks.append(matrix(field, 2, [[a,-b],[b,a]]))
1129
1130     # We can drop the imaginaries here.
1131     return block_matrix(field.base_ring(), n, blocks)
1132
1133
1134 def _unembed_complex_matrix(M):
1135     """
1136     The inverse of _embed_complex_matrix().
1137
1138     EXAMPLES::
1139
1140         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1141         ....:                 [-2,  1,  -4,  3],
1142         ....:                 [ 9,  10, 11, 12],
1143         ....:                 [-10, 9, -12, 11] ])
1144         sage: _unembed_complex_matrix(A)
1145         [  -2*i + 1   -4*i + 3]
1146         [ -10*i + 9 -12*i + 11]
1147     """
1148     n = ZZ(M.nrows())
1149     if M.ncols() != n:
1150         raise ArgumentError("the matrix 'M' must be square")
1151     if not n.mod(2).is_zero():
1152         raise ArgumentError("the matrix 'M' must be a complex embedding")
1153
1154     F = QuadraticField(-1, 'i')
1155     i = F.gen()
1156
1157     # Go top-left to bottom-right (reading order), converting every
1158     # 2-by-2 block we see to a single complex element.
1159     elements = []
1160     for k in xrange(n/2):
1161         for j in xrange(n/2):
1162             submat = M[2*k:2*k+2,2*j:2*j+2]
1163             if submat[0,0] != submat[1,1]:
1164                 raise ArgumentError('bad real submatrix')
1165             if submat[0,1] != -submat[1,0]:
1166                 raise ArgumentError('bad imag submatrix')
1167             z = submat[0,0] + submat[1,0]*i
1168             elements.append(z)
1169
1170     return matrix(F, n/2, elements)
1171
1172
1173 def RealSymmetricSimpleEJA(n, field=QQ):
1174     """
1175     The rank-n simple EJA consisting of real symmetric n-by-n
1176     matrices, the usual symmetric Jordan product, and the trace inner
1177     product. It has dimension `(n^2 + n)/2` over the reals.
1178
1179     EXAMPLES::
1180
1181         sage: J = RealSymmetricSimpleEJA(2)
1182         sage: e0, e1, e2 = J.gens()
1183         sage: e0*e0
1184         e0
1185         sage: e1*e1
1186         e0 + e2
1187         sage: e2*e2
1188         e2
1189
1190     TESTS:
1191
1192     The degree of this algebra is `(n^2 + n) / 2`::
1193
1194         sage: set_random_seed()
1195         sage: n = ZZ.random_element(1,5)
1196         sage: J = RealSymmetricSimpleEJA(n)
1197         sage: J.degree() == (n^2 + n)/2
1198         True
1199
1200     """
1201     S = _real_symmetric_basis(n, field=field)
1202     (Qs, T) = _multiplication_table_from_matrix_basis(S)
1203
1204     return FiniteDimensionalEuclideanJordanAlgebra(field,
1205                                                    Qs,
1206                                                    rank=n,
1207                                                    natural_basis=T)
1208
1209
1210 def ComplexHermitianSimpleEJA(n, field=QQ):
1211     """
1212     The rank-n simple EJA consisting of complex Hermitian n-by-n
1213     matrices over the real numbers, the usual symmetric Jordan product,
1214     and the real-part-of-trace inner product. It has dimension `n^2` over
1215     the reals.
1216
1217     TESTS:
1218
1219     The degree of this algebra is `n^2`::
1220
1221         sage: set_random_seed()
1222         sage: n = ZZ.random_element(1,5)
1223         sage: J = ComplexHermitianSimpleEJA(n)
1224         sage: J.degree() == n^2
1225         True
1226
1227     """
1228     S = _complex_hermitian_basis(n)
1229     (Qs, T) = _multiplication_table_from_matrix_basis(S)
1230     return FiniteDimensionalEuclideanJordanAlgebra(field,
1231                                                    Qs,
1232                                                    rank=n,
1233                                                    natural_basis=T)
1234
1235
1236 def QuaternionHermitianSimpleEJA(n):
1237     """
1238     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1239     matrices, the usual symmetric Jordan product, and the
1240     real-part-of-trace inner product. It has dimension `2n^2 - n` over
1241     the reals.
1242     """
1243     pass
1244
1245 def OctonionHermitianSimpleEJA(n):
1246     """
1247     This shit be crazy. It has dimension 27 over the reals.
1248     """
1249     n = 3
1250     pass
1251
1252 def JordanSpinSimpleEJA(n, field=QQ):
1253     """
1254     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1255     with the usual inner product and jordan product ``x*y =
1256     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1257     the reals.
1258
1259     EXAMPLES:
1260
1261     This multiplication table can be verified by hand::
1262
1263         sage: J = JordanSpinSimpleEJA(4)
1264         sage: e0,e1,e2,e3 = J.gens()
1265         sage: e0*e0
1266         e0
1267         sage: e0*e1
1268         e1
1269         sage: e0*e2
1270         e2
1271         sage: e0*e3
1272         e3
1273         sage: e1*e2
1274         0
1275         sage: e1*e3
1276         0
1277         sage: e2*e3
1278         0
1279
1280     """
1281     Qs = []
1282     id_matrix = identity_matrix(field, n)
1283     for i in xrange(n):
1284         ei = id_matrix.column(i)
1285         Qi = zero_matrix(field, n)
1286         Qi.set_row(0, ei)
1287         Qi.set_column(0, ei)
1288         Qi += diagonal_matrix(n, [ei[0]]*n)
1289         # The addition of the diagonal matrix adds an extra ei[0] in the
1290         # upper-left corner of the matrix.
1291         Qi[0,0] = Qi[0,0] * ~field(2)
1292         Qs.append(Qi)
1293
1294     # The usual inner product on R^n.
1295     ip = lambda x, y: x.vector().inner_product(y.vector())
1296
1297     # The rank of the spin factor algebra is two, UNLESS we're in a
1298     # one-dimensional ambient space (the rank is bounded by the
1299     # ambient dimension).
1300     return FiniteDimensionalEuclideanJordanAlgebra(field,
1301                                                    Qs,
1302                                                    rank=min(n,2),
1303                                                    inner_product=ip)