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