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