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         def minimal_polynomial(self):
 465             """
 466             EXAMPLES::
 467
 468                 sage: set_random_seed()
 469                 sage: x = random_eja().random_element()
 470                 sage: x.degree() == x.minimal_polynomial().degree()
 471                 True
 472
 473             ::
 474
 475                 sage: set_random_seed()
 476                 sage: x = random_eja().random_element()
 477                 sage: x.degree() == x.minimal_polynomial().degree()
 478                 True
 479
 480             The minimal polynomial and the characteristic polynomial coincide
 481             and are known (see Alizadeh, Example 11.11) for all elements of
 482             the spin factor algebra that aren't scalar multiples of the
 483             identity::
 484
 485                 sage: set_random_seed()
 486                 sage: n = ZZ.random_element(2,10)
 487                 sage: J = JordanSpinSimpleEJA(n)
 488                 sage: y = J.random_element()
 489                 sage: while y == y.coefficient(0)*J.one():
 490                 ....:     y = J.random_element()
 491                 sage: y0 = y.vector()[0]
 492                 sage: y_bar = y.vector()[1:]
 493                 sage: actual = y.minimal_polynomial()
 494                 sage: x = SR.symbol('x', domain='real')
 495                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 496                 sage: bool(actual == expected)
 497                 True
 498
 499             """
 500             # The element we're going to call "minimal_polynomial()" on.
 501             # Either myself, interpreted as an element of a finite-
 502             # dimensional algebra, or an element of an associative
 503             # subalgebra.
 504             elt = None
 505
 506             if self.parent().is_associative():
 507                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 508             else:
 509                 V = self.span_of_powers()
 510                 assoc_subalg = self.subalgebra_generated_by()
 511                 # Mis-design warning: the basis used for span_of_powers()
 512                 # and subalgebra_generated_by() must be the same, and in
 513                 # the same order!
 514                 elt = assoc_subalg(V.coordinates(self.vector()))
 515
 516             # Recursive call, but should work since elt lives in an
 517             # associative algebra.
 518             return elt.minimal_polynomial()
 519
 520
 521         def natural_representation(self):
 522             """
 523             Return a more-natural representation of this element.
 524
 525             Every finite-dimensional Euclidean Jordan Algebra is a
 526             direct sum of five simple algebras, four of which comprise
 527             Hermitian matrices. This method returns the original
 528             "natural" representation of this element as a Hermitian
 529             matrix, if it has one. If not, you get the usual representation.
 530
 531             EXAMPLES::
 532
 533                 sage: J = ComplexHermitianSimpleEJA(3)
 534                 sage: J.one()
 535                 e0 + e5 + e8
 536                 sage: J.one().natural_representation()
 537                 [1 0 0 0 0 0]
 538                 [0 1 0 0 0 0]
 539                 [0 0 1 0 0 0]
 540                 [0 0 0 1 0 0]
 541                 [0 0 0 0 1 0]
 542                 [0 0 0 0 0 1]
 543
 544             """
 545             B = self.parent().natural_basis()
 546             W = B[0].matrix_space()
 547             return W.linear_combination(zip(self.vector(), B))
 548
 549
 550         def operator_matrix(self):
 551             """
 552             Return the matrix that represents left- (or right-)
 553             multiplication by this element in the parent algebra.
 554
 555             We have to override this because the superclass method
 556             returns a matrix that acts on row vectors (that is, on
 557             the right).
 558
 559             EXAMPLES:
 560
 561             Test the first polarization identity from my notes, Koecher Chapter
 562             III, or from Baes (2.3)::
 563
 564                 sage: set_random_seed()
 565                 sage: J = random_eja()
 566                 sage: x = J.random_element()
 567                 sage: y = J.random_element()
 568                 sage: Lx = x.operator_matrix()
 569                 sage: Ly = y.operator_matrix()
 570                 sage: Lxx = (x*x).operator_matrix()
 571                 sage: Lxy = (x*y).operator_matrix()
 572                 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
 573                 True
 574
 575             Test the second polarization identity from my notes or from
 576             Baes (2.4)::
 577
 578                 sage: set_random_seed()
 579                 sage: J = random_eja()
 580                 sage: x = J.random_element()
 581                 sage: y = J.random_element()
 582                 sage: z = J.random_element()
 583                 sage: Lx = x.operator_matrix()
 584                 sage: Ly = y.operator_matrix()
 585                 sage: Lz = z.operator_matrix()
 586                 sage: Lzy = (z*y).operator_matrix()
 587                 sage: Lxy = (x*y).operator_matrix()
 588                 sage: Lxz = (x*z).operator_matrix()
 589                 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
 590                 True
 591
 592             Test the third polarization identity from my notes or from
 593             Baes (2.5)::
 594
 595                 sage: set_random_seed()
 596                 sage: J = random_eja()
 597                 sage: u = J.random_element()
 598                 sage: y = J.random_element()
 599                 sage: z = J.random_element()
 600                 sage: Lu = u.operator_matrix()
 601                 sage: Ly = y.operator_matrix()
 602                 sage: Lz = z.operator_matrix()
 603                 sage: Lzy = (z*y).operator_matrix()
 604                 sage: Luy = (u*y).operator_matrix()
 605                 sage: Luz = (u*z).operator_matrix()
 606                 sage: Luyz = (u*(y*z)).operator_matrix()
 607                 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
 608                 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
 609                 sage: bool(lhs == rhs)
 610                 True
 611
 612             """
 613             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 614             return fda_elt.matrix().transpose()
 615
 616
 617         def quadratic_representation(self, other=None):
 618             """
 619             Return the quadratic representation of this element.
 620
 621             EXAMPLES:
 622
 623             The explicit form in the spin factor algebra is given by
 624             Alizadeh's Example 11.12::
 625
 626                 sage: set_random_seed()
 627                 sage: n = ZZ.random_element(1,10)
 628                 sage: J = JordanSpinSimpleEJA(n)
 629                 sage: x = J.random_element()
 630                 sage: x_vec = x.vector()
 631                 sage: x0 = x_vec[0]
 632                 sage: x_bar = x_vec[1:]
 633                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 634                 sage: B = 2*x0*x_bar.row()
 635                 sage: C = 2*x0*x_bar.column()
 636                 sage: D = identity_matrix(QQ, n-1)
 637                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 638                 sage: D = D + 2*x_bar.tensor_product(x_bar)
 639                 sage: Q = block_matrix(2,2,[A,B,C,D])
 640                 sage: Q == x.quadratic_representation()
 641                 True
 642
 643             Test all of the properties from Theorem 11.2 in Alizadeh::
 644
 645                 sage: set_random_seed()
 646                 sage: J = random_eja()
 647                 sage: x = J.random_element()
 648                 sage: y = J.random_element()
 649
 650             Property 1:
 651
 652                 sage: actual = x.quadratic_representation(y)
 653                 sage: expected = ( (x+y).quadratic_representation()
 654                 ....:              -x.quadratic_representation()
 655                 ....:              -y.quadratic_representation() ) / 2
 656                 sage: actual == expected
 657                 True
 658
 659             Property 2:
 660
 661                 sage: alpha = QQ.random_element()
 662                 sage: actual = (alpha*x).quadratic_representation()
 663                 sage: expected = (alpha^2)*x.quadratic_representation()
 664                 sage: actual == expected
 665                 True
 666
 667             Property 5:
 668
 669                 sage: Qy = y.quadratic_representation()
 670                 sage: actual = J(Qy*x.vector()).quadratic_representation()
 671                 sage: expected = Qy*x.quadratic_representation()*Qy
 672                 sage: actual == expected
 673                 True
 674
 675             Property 6:
 676
 677                 sage: k = ZZ.random_element(1,10)
 678                 sage: actual = (x^k).quadratic_representation()
 679                 sage: expected = (x.quadratic_representation())^k
 680                 sage: actual == expected
 681                 True
 682
 683             """
 684             if other is None:
 685                 other=self
 686             elif not other in self.parent():
 687                 raise ArgumentError("'other' must live in the same algebra")
 688
 689             L = self.operator_matrix()
 690             M = other.operator_matrix()
 691             return ( L*M + M*L - (self*other).operator_matrix() )
 692
 693
 694         def span_of_powers(self):
 695             """
 696             Return the vector space spanned by successive powers of
 697             this element.
 698             """
 699             # The dimension of the subalgebra can't be greater than
 700             # the big algebra, so just put everything into a list
 701             # and let span() get rid of the excess.
 702             V = self.vector().parent()
 703             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 704
 705
 706         def subalgebra_generated_by(self):
 707             """
 708             Return the associative subalgebra of the parent EJA generated
 709             by this element.
 710
 711             TESTS::
 712
 713                 sage: set_random_seed()
 714                 sage: x = random_eja().random_element()
 715                 sage: x.subalgebra_generated_by().is_associative()
 716                 True
 717
 718             Squaring in the subalgebra should be the same thing as
 719             squaring in the superalgebra::
 720
 721                 sage: set_random_seed()
 722                 sage: x = random_eja().random_element()
 723                 sage: u = x.subalgebra_generated_by().random_element()
 724                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
 725                 True
 726
 727             """
 728             # First get the subspace spanned by the powers of myself...
 729             V = self.span_of_powers()
 730             F = self.base_ring()
 731
 732             # Now figure out the entries of the right-multiplication
 733             # matrix for the successive basis elements b0, b1,... of
 734             # that subspace.
 735             mats = []
 736             for b_right in V.basis():
 737                 eja_b_right = self.parent()(b_right)
 738                 b_right_rows = []
 739                 # The first row of the right-multiplication matrix by
 740                 # b1 is what we get if we apply that matrix to b1. The
 741                 # second row of the right multiplication matrix by b1
 742                 # is what we get when we apply that matrix to b2...
 743                 #
 744                 # IMPORTANT: this assumes that all vectors are COLUMN
 745                 # vectors, unlike our superclass (which uses row vectors).
 746                 for b_left in V.basis():
 747                     eja_b_left = self.parent()(b_left)
 748                     # Multiply in the original EJA, but then get the
 749                     # coordinates from the subalgebra in terms of its
 750                     # basis.
 751                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 752                     b_right_rows.append(this_row)
 753                 b_right_matrix = matrix(F, b_right_rows)
 754                 mats.append(b_right_matrix)
 755
 756             # It's an algebra of polynomials in one element, and EJAs
 757             # are power-associative.
 758             #
 759             # TODO: choose generator names intelligently.
 760             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 761
 762
 763         def subalgebra_idempotent(self):
 764             """
 765             Find an idempotent in the associative subalgebra I generate
 766             using Proposition 2.3.5 in Baes.
 767
 768             TESTS::
 769
 770                 sage: set_random_seed()
 771                 sage: J = eja_rn(5)
 772                 sage: c = J.random_element().subalgebra_idempotent()
 773                 sage: c^2 == c
 774                 True
 775                 sage: J = JordanSpinSimpleEJA(5)
 776                 sage: c = J.random_element().subalgebra_idempotent()
 777                 sage: c^2 == c
 778                 True
 779
 780             """
 781             if self.is_nilpotent():
 782                 raise ValueError("this only works with non-nilpotent elements!")
 783
 784             V = self.span_of_powers()
 785             J = self.subalgebra_generated_by()
 786             # Mis-design warning: the basis used for span_of_powers()
 787             # and subalgebra_generated_by() must be the same, and in
 788             # the same order!
 789             u = J(V.coordinates(self.vector()))
 790
 791             # The image of the matrix of left-u^m-multiplication
 792             # will be minimal for some natural number s...
 793             s = 0
 794             minimal_dim = V.dimension()
 795             for i in xrange(1, V.dimension()):
 796                 this_dim = (u**i).operator_matrix().image().dimension()
 797                 if this_dim < minimal_dim:
 798                     minimal_dim = this_dim
 799                     s = i
 800
 801             # Now minimal_matrix should correspond to the smallest
 802             # non-zero subspace in Baes's (or really, Koecher's)
 803             # proposition.
 804             #
 805             # However, we need to restrict the matrix to work on the
 806             # subspace... or do we? Can't we just solve, knowing that
 807             # A(c) = u^(s+1) should have a solution in the big space,
 808             # too?
 809             #
 810             # Beware, solve_right() means that we're using COLUMN vectors.
 811             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 812             u_next = u**(s+1)
 813             A = u_next.operator_matrix()
 814             c_coordinates = A.solve_right(u_next.vector())
 815
 816             # Now c_coordinates is the idempotent we want, but it's in
 817             # the coordinate system of the subalgebra.
 818             #
 819             # We need the basis for J, but as elements of the parent algebra.
 820             #
 821             basis = [self.parent(v) for v in V.basis()]
 822             return self.parent().linear_combination(zip(c_coordinates, basis))
 823
 824
 825         def trace(self):
 826             """
 827             Return my trace, the sum of my eigenvalues.
 828
 829             EXAMPLES::
 830
 831                 sage: J = JordanSpinSimpleEJA(3)
 832                 sage: e0,e1,e2 = J.gens()
 833                 sage: x = e0 + e1 + e2
 834                 sage: x.trace()
 835                 2
 836
 837             """
 838             cs = self.characteristic_polynomial().coefficients(sparse=False)
 839             if len(cs) >= 2:
 840                 return -1*cs[-2]
 841             else:
 842                 raise ValueError('charpoly had fewer than 2 coefficients')
 843
 844
 845         def trace_inner_product(self, other):
 846             """
 847             Return the trace inner product of myself and ``other``.
 848             """
 849             if not other in self.parent():
 850                 raise ArgumentError("'other' must live in the same algebra")
 851
 852             return (self*other).trace()
 853
 854
 855 def eja_rn(dimension, field=QQ):
 856     """
 857     Return the Euclidean Jordan Algebra corresponding to the set
 858     `R^n` under the Hadamard product.
 859
 860     EXAMPLES:
 861
 862     This multiplication table can be verified by hand::
 863
 864         sage: J = eja_rn(3)
 865         sage: e0,e1,e2 = J.gens()
 866         sage: e0*e0
 867         e0
 868         sage: e0*e1
 869         0
 870         sage: e0*e2
 871         0
 872         sage: e1*e1
 873         e1
 874         sage: e1*e2
 875         0
 876         sage: e2*e2
 877         e2
 878
 879     """
 880     # The FiniteDimensionalAlgebra constructor takes a list of
 881     # matrices, the ith representing right multiplication by the ith
 882     # basis element in the vector space. So if e_1 = (1,0,0), then
 883     # right (Hadamard) multiplication of x by e_1 picks out the first
 884     # component of x; and likewise for the ith basis element e_i.
 885     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
 886            for i in xrange(dimension) ]
 887
 888     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 889
 890
 891
 892 def random_eja():
 893     """
 894     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 895
 896     ALGORITHM:
 897
 898     For now, we choose a random natural number ``n`` (greater than zero)
 899     and then give you back one of the following:
 900
 901       * The cartesian product of the rational numbers ``n`` times; this is
 902         ``QQ^n`` with the Hadamard product.
 903
 904       * The Jordan spin algebra on ``QQ^n``.
 905
 906       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
 907         product.
 908
 909     Later this might be extended to return Cartesian products of the
 910     EJAs above.
 911
 912     TESTS::
 913
 914         sage: random_eja()
 915         Euclidean Jordan algebra of degree...
 916
 917     """
 918     n = ZZ.random_element(1,5)
 919     constructor = choice([eja_rn,
 920                           JordanSpinSimpleEJA,
 921                           RealSymmetricSimpleEJA,
 922                           ComplexHermitianSimpleEJA])
 923     return constructor(n, field=QQ)
 924
 925
 926
 927 def _real_symmetric_basis(n, field=QQ):
 928     """
 929     Return a basis for the space of real symmetric n-by-n matrices.
 930     """
 931     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
 932     # coordinates.
 933     S = []
 934     for i in xrange(n):
 935         for j in xrange(i+1):
 936             Eij = matrix(field, n, lambda k,l: k==i and l==j)
 937             if i == j:
 938                 Sij = Eij
 939             else:
 940                 # Beware, orthogonal but not normalized!
 941                 Sij = Eij + Eij.transpose()
 942             S.append(Sij)
 943     return tuple(S)
 944
 945
 946 def _complex_hermitian_basis(n, field=QQ):
 947     """
 948     Returns a basis for the space of complex Hermitian n-by-n matrices.
 949
 950     TESTS::
 951
 952         sage: set_random_seed()
 953         sage: n = ZZ.random_element(1,5)
 954         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
 955         True
 956
 957     """
 958     F = QuadraticField(-1, 'I')
 959     I = F.gen()
 960
 961     # This is like the symmetric case, but we need to be careful:
 962     #
 963     #   * We want conjugate-symmetry, not just symmetry.
 964     #   * The diagonal will (as a result) be real.
 965     #
 966     S = []
 967     for i in xrange(n):
 968         for j in xrange(i+1):
 969             Eij = matrix(field, n, lambda k,l: k==i and l==j)
 970             if i == j:
 971                 Sij = _embed_complex_matrix(Eij)
 972                 S.append(Sij)
 973             else:
 974                 # Beware, orthogonal but not normalized! The second one
 975                 # has a minus because it's conjugated.
 976                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
 977                 S.append(Sij_real)
 978                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
 979                 S.append(Sij_imag)
 980     return tuple(S)
 981
 982
 983 def _multiplication_table_from_matrix_basis(basis):
 984     """
 985     At least three of the five simple Euclidean Jordan algebras have the
 986     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
 987     multiplication on the right is matrix multiplication. Given a basis
 988     for the underlying matrix space, this function returns a
 989     multiplication table (obtained by looping through the basis
 990     elements) for an algebra of those matrices. A reordered copy
 991     of the basis is also returned to work around the fact that
 992     the ``span()`` in this function will change the order of the basis
 993     from what we think it is, to... something else.
 994     """
 995     # In S^2, for example, we nominally have four coordinates even
 996     # though the space is of dimension three only. The vector space V
 997     # is supposed to hold the entire long vector, and the subspace W
 998     # of V will be spanned by the vectors that arise from symmetric
 999     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1000     field = basis[0].base_ring()
1001     dimension = basis[0].nrows()
1002
1003     def mat2vec(m):
1004         return vector(field, m.list())
1005
1006     def vec2mat(v):
1007         return matrix(field, dimension, v.list())
1008
1009     V = VectorSpace(field, dimension**2)
1010     W = V.span( mat2vec(s) for s in basis )
1011
1012     # Taking the span above reorders our basis (thanks, jerk!) so we
1013     # need to put our "matrix basis" in the same order as the
1014     # (reordered) vector basis.
1015     S = tuple( vec2mat(b) for b in W.basis() )
1016
1017     Qs = []
1018     for s in S:
1019         # Brute force the multiplication-by-s matrix by looping
1020         # through all elements of the basis and doing the computation
1021         # to find out what the corresponding row should be. BEWARE:
1022         # these multiplication tables won't be symmetric! It therefore
1023         # becomes REALLY IMPORTANT that the underlying algebra
1024         # constructor uses ROW vectors and not COLUMN vectors. That's
1025         # why we're computing rows here and not columns.
1026         Q_rows = []
1027         for t in S:
1028             this_row = mat2vec((s*t + t*s)/2)
1029             Q_rows.append(W.coordinates(this_row))
1030         Q = matrix(field, W.dimension(), Q_rows)
1031         Qs.append(Q)
1032
1033     return (Qs, S)
1034
1035
1036 def _embed_complex_matrix(M):
1037     """
1038     Embed the n-by-n complex matrix ``M`` into the space of real
1039     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1040     bi` to the block matrix ``[[a,b],[-b,a]]``.
1041
1042     EXAMPLES::
1043
1044         sage: F = QuadraticField(-1,'i')
1045         sage: x1 = F(4 - 2*i)
1046         sage: x2 = F(1 + 2*i)
1047         sage: x3 = F(-i)
1048         sage: x4 = F(6)
1049         sage: M = matrix(F,2,[x1,x2,x3,x4])
1050         sage: _embed_complex_matrix(M)
1051         [ 4  2| 1 -2]
1052         [-2  4| 2  1]
1053         [-----+-----]
1054         [ 0  1| 6  0]
1055         [-1  0| 0  6]
1056
1057     """
1058     n = M.nrows()
1059     if M.ncols() != n:
1060         raise ArgumentError("the matrix 'M' must be square")
1061     field = M.base_ring()
1062     blocks = []
1063     for z in M.list():
1064         a = z.real()
1065         b = z.imag()
1066         blocks.append(matrix(field, 2, [[a,-b],[b,a]]))
1067
1068     # We can drop the imaginaries here.
1069     return block_matrix(field.base_ring(), n, blocks)
1070
1071
1072 def _unembed_complex_matrix(M):
1073     """
1074     The inverse of _embed_complex_matrix().
1075
1076     EXAMPLES::
1077
1078         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1079         ....:                 [-2,  1,  -4,  3],
1080         ....:                 [ 9,  10, 11, 12],
1081         ....:                 [-10, 9, -12, 11] ])
1082         sage: _unembed_complex_matrix(A)
1083         [  -2*i + 1   -4*i + 3]
1084         [ -10*i + 9 -12*i + 11]
1085     """
1086     n = ZZ(M.nrows())
1087     if M.ncols() != n:
1088         raise ArgumentError("the matrix 'M' must be square")
1089     if not n.mod(2).is_zero():
1090         raise ArgumentError("the matrix 'M' must be a complex embedding")
1091
1092     F = QuadraticField(-1, 'i')
1093     i = F.gen()
1094
1095     # Go top-left to bottom-right (reading order), converting every
1096     # 2-by-2 block we see to a single complex element.
1097     elements = []
1098     for k in xrange(n/2):
1099         for j in xrange(n/2):
1100             submat = M[2*k:2*k+2,2*j:2*j+2]
1101             if submat[0,0] != submat[1,1]:
1102                 raise ArgumentError('bad real submatrix')
1103             if submat[0,1] != -submat[1,0]:
1104                 raise ArgumentError('bad imag submatrix')
1105             z = submat[0,0] + submat[1,0]*i
1106             elements.append(z)
1107
1108     return matrix(F, n/2, elements)
1109
1110
1111 def RealSymmetricSimpleEJA(n, field=QQ):
1112     """
1113     The rank-n simple EJA consisting of real symmetric n-by-n
1114     matrices, the usual symmetric Jordan product, and the trace inner
1115     product. It has dimension `(n^2 + n)/2` over the reals.
1116
1117     EXAMPLES::
1118
1119         sage: J = RealSymmetricSimpleEJA(2)
1120         sage: e0, e1, e2 = J.gens()
1121         sage: e0*e0
1122         e0
1123         sage: e1*e1
1124         e0 + e2
1125         sage: e2*e2
1126         e2
1127
1128     TESTS:
1129
1130     The degree of this algebra is `(n^2 + n) / 2`::
1131
1132         sage: set_random_seed()
1133         sage: n = ZZ.random_element(1,5)
1134         sage: J = RealSymmetricSimpleEJA(n)
1135         sage: J.degree() == (n^2 + n)/2
1136         True
1137
1138     """
1139     S = _real_symmetric_basis(n, field=field)
1140     (Qs, T) = _multiplication_table_from_matrix_basis(S)
1141
1142     return FiniteDimensionalEuclideanJordanAlgebra(field,
1143                                                    Qs,
1144                                                    rank=n,
1145                                                    natural_basis=T)
1146
1147
1148 def ComplexHermitianSimpleEJA(n, field=QQ):
1149     """
1150     The rank-n simple EJA consisting of complex Hermitian n-by-n
1151     matrices over the real numbers, the usual symmetric Jordan product,
1152     and the real-part-of-trace inner product. It has dimension `n^2` over
1153     the reals.
1154
1155     TESTS:
1156
1157     The degree of this algebra is `n^2`::
1158
1159         sage: set_random_seed()
1160         sage: n = ZZ.random_element(1,5)
1161         sage: J = ComplexHermitianSimpleEJA(n)
1162         sage: J.degree() == n^2
1163         True
1164
1165     """
1166     S = _complex_hermitian_basis(n)
1167     (Qs, T) = _multiplication_table_from_matrix_basis(S)
1168     return FiniteDimensionalEuclideanJordanAlgebra(field,
1169                                                    Qs,
1170                                                    rank=n,
1171                                                    natural_basis=T)
1172
1173
1174 def QuaternionHermitianSimpleEJA(n):
1175     """
1176     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1177     matrices, the usual symmetric Jordan product, and the
1178     real-part-of-trace inner product. It has dimension `2n^2 - n` over
1179     the reals.
1180     """
1181     pass
1182
1183 def OctonionHermitianSimpleEJA(n):
1184     """
1185     This shit be crazy. It has dimension 27 over the reals.
1186     """
1187     n = 3
1188     pass
1189
1190 def JordanSpinSimpleEJA(n, field=QQ):
1191     """
1192     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1193     with the usual inner product and jordan product ``x*y =
1194     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1195     the reals.
1196
1197     EXAMPLES:
1198
1199     This multiplication table can be verified by hand::
1200
1201         sage: J = JordanSpinSimpleEJA(4)
1202         sage: e0,e1,e2,e3 = J.gens()
1203         sage: e0*e0
1204         e0
1205         sage: e0*e1
1206         e1
1207         sage: e0*e2
1208         e2
1209         sage: e0*e3
1210         e3
1211         sage: e1*e2
1212         0
1213         sage: e1*e3
1214         0
1215         sage: e2*e3
1216         0
1217
1218     In one dimension, this is the reals under multiplication::
1219
1220       sage: J1 = JordanSpinSimpleEJA(1)
1221       sage: J2 = eja_rn(1)
1222       sage: J1 == J2
1223       True
1224
1225     """
1226     Qs = []
1227     id_matrix = identity_matrix(field, n)
1228     for i in xrange(n):
1229         ei = id_matrix.column(i)
1230         Qi = zero_matrix(field, n)
1231         Qi.set_row(0, ei)
1232         Qi.set_column(0, ei)
1233         Qi += diagonal_matrix(n, [ei[0]]*n)
1234         # The addition of the diagonal matrix adds an extra ei[0] in the
1235         # upper-left corner of the matrix.
1236         Qi[0,0] = Qi[0,0] * ~field(2)
1237         Qs.append(Qi)
1238
1239     # The rank of the spin factor algebra is two, UNLESS we're in a
1240     # one-dimensional ambient space (the rank is bounded by the
1241     # ambient dimension).
1242     return FiniteDimensionalEuclideanJordanAlgebra(field, Qs, rank=min(n,2))