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