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