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