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 = eja_ln(2)
 137                 sage: e0,e1 = J.gens()
 138                 sage: x = e0 + e1
 139                 sage: x.det()
 140                 0
 141                 sage: J = eja_ln(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 = eja_ln(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 = eja_ln(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 = eja_ln(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 = eja_ln(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):
 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: n = ZZ.random_element(1,10).abs()
 336                 sage: J = eja_ln(n)
 337                 sage: x = J.random_element()
 338                 sage: x_vec = x.vector()
 339                 sage: x0 = x_vec[0]
 340                 sage: x_bar = x_vec[1:]
 341                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 342                 sage: B = 2*x0*x_bar.row()
 343                 sage: C = 2*x0*x_bar.column()
 344                 sage: D = identity_matrix(QQ, n-1)
 345                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 346                 sage: D = D + 2*x_bar.tensor_product(x_bar)
 347                 sage: Q = block_matrix(2,2,[A,B,C,D])
 348                 sage: Q == x.quadratic_representation()
 349                 True
 350
 351             """
 352             return 2*(self.matrix()**2) - (self**2).matrix()
 353
 354
 355         def span_of_powers(self):
 356             """
 357             Return the vector space spanned by successive powers of
 358             this element.
 359             """
 360             # The dimension of the subalgebra can't be greater than
 361             # the big algebra, so just put everything into a list
 362             # and let span() get rid of the excess.
 363             V = self.vector().parent()
 364             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 365
 366
 367         def subalgebra_generated_by(self):
 368             """
 369             Return the associative subalgebra of the parent EJA generated
 370             by this element.
 371
 372             TESTS::
 373
 374                 sage: set_random_seed()
 375                 sage: x = random_eja().random_element()
 376                 sage: x.subalgebra_generated_by().is_associative()
 377                 True
 378
 379             Squaring in the subalgebra should be the same thing as
 380             squaring in the superalgebra::
 381
 382                 sage: set_random_seed()
 383                 sage: x = random_eja().random_element()
 384                 sage: u = x.subalgebra_generated_by().random_element()
 385                 sage: u.matrix()*u.vector() == (u**2).vector()
 386                 True
 387
 388             """
 389             # First get the subspace spanned by the powers of myself...
 390             V = self.span_of_powers()
 391             F = self.base_ring()
 392
 393             # Now figure out the entries of the right-multiplication
 394             # matrix for the successive basis elements b0, b1,... of
 395             # that subspace.
 396             mats = []
 397             for b_right in V.basis():
 398                 eja_b_right = self.parent()(b_right)
 399                 b_right_rows = []
 400                 # The first row of the right-multiplication matrix by
 401                 # b1 is what we get if we apply that matrix to b1. The
 402                 # second row of the right multiplication matrix by b1
 403                 # is what we get when we apply that matrix to b2...
 404                 #
 405                 # IMPORTANT: this assumes that all vectors are COLUMN
 406                 # vectors, unlike our superclass (which uses row vectors).
 407                 for b_left in V.basis():
 408                     eja_b_left = self.parent()(b_left)
 409                     # Multiply in the original EJA, but then get the
 410                     # coordinates from the subalgebra in terms of its
 411                     # basis.
 412                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 413                     b_right_rows.append(this_row)
 414                 b_right_matrix = matrix(F, b_right_rows)
 415                 mats.append(b_right_matrix)
 416
 417             # It's an algebra of polynomials in one element, and EJAs
 418             # are power-associative.
 419             #
 420             # TODO: choose generator names intelligently.
 421             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 422
 423
 424         def subalgebra_idempotent(self):
 425             """
 426             Find an idempotent in the associative subalgebra I generate
 427             using Proposition 2.3.5 in Baes.
 428
 429             TESTS::
 430
 431                 sage: set_random_seed()
 432                 sage: J = eja_rn(5)
 433                 sage: c = J.random_element().subalgebra_idempotent()
 434                 sage: c^2 == c
 435                 True
 436                 sage: J = eja_ln(5)
 437                 sage: c = J.random_element().subalgebra_idempotent()
 438                 sage: c^2 == c
 439                 True
 440
 441             """
 442             if self.is_nilpotent():
 443                 raise ValueError("this only works with non-nilpotent elements!")
 444
 445             V = self.span_of_powers()
 446             J = self.subalgebra_generated_by()
 447             # Mis-design warning: the basis used for span_of_powers()
 448             # and subalgebra_generated_by() must be the same, and in
 449             # the same order!
 450             u = J(V.coordinates(self.vector()))
 451
 452             # The image of the matrix of left-u^m-multiplication
 453             # will be minimal for some natural number s...
 454             s = 0
 455             minimal_dim = V.dimension()
 456             for i in xrange(1, V.dimension()):
 457                 this_dim = (u**i).matrix().image().dimension()
 458                 if this_dim < minimal_dim:
 459                     minimal_dim = this_dim
 460                     s = i
 461
 462             # Now minimal_matrix should correspond to the smallest
 463             # non-zero subspace in Baes's (or really, Koecher's)
 464             # proposition.
 465             #
 466             # However, we need to restrict the matrix to work on the
 467             # subspace... or do we? Can't we just solve, knowing that
 468             # A(c) = u^(s+1) should have a solution in the big space,
 469             # too?
 470             #
 471             # Beware, solve_right() means that we're using COLUMN vectors.
 472             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 473             u_next = u**(s+1)
 474             A = u_next.matrix()
 475             c_coordinates = A.solve_right(u_next.vector())
 476
 477             # Now c_coordinates is the idempotent we want, but it's in
 478             # the coordinate system of the subalgebra.
 479             #
 480             # We need the basis for J, but as elements of the parent algebra.
 481             #
 482             basis = [self.parent(v) for v in V.basis()]
 483             return self.parent().linear_combination(zip(c_coordinates, basis))
 484
 485
 486         def trace(self):
 487             """
 488             Return my trace, the sum of my eigenvalues.
 489
 490             EXAMPLES::
 491
 492                 sage: J = eja_ln(3)
 493                 sage: e0,e1,e2 = J.gens()
 494                 sage: x = e0 + e1 + e2
 495                 sage: x.trace()
 496                 2
 497
 498             """
 499             cs = self.characteristic_polynomial().coefficients(sparse=False)
 500             if len(cs) >= 2:
 501                 return -1*cs[-2]
 502             else:
 503                 raise ValueError('charpoly had fewer than 2 coefficients')
 504
 505
 506         def trace_inner_product(self, other):
 507             """
 508             Return the trace inner product of myself and ``other``.
 509             """
 510             if not other in self.parent():
 511                 raise ArgumentError("'other' must live in the same algebra")
 512
 513             return (self*other).trace()
 514
 515
 516 def eja_rn(dimension, field=QQ):
 517     """
 518     Return the Euclidean Jordan Algebra corresponding to the set
 519     `R^n` under the Hadamard product.
 520
 521     EXAMPLES:
 522
 523     This multiplication table can be verified by hand::
 524
 525         sage: J = eja_rn(3)
 526         sage: e0,e1,e2 = J.gens()
 527         sage: e0*e0
 528         e0
 529         sage: e0*e1
 530         0
 531         sage: e0*e2
 532         0
 533         sage: e1*e1
 534         e1
 535         sage: e1*e2
 536         0
 537         sage: e2*e2
 538         e2
 539
 540     """
 541     # The FiniteDimensionalAlgebra constructor takes a list of
 542     # matrices, the ith representing right multiplication by the ith
 543     # basis element in the vector space. So if e_1 = (1,0,0), then
 544     # right (Hadamard) multiplication of x by e_1 picks out the first
 545     # component of x; and likewise for the ith basis element e_i.
 546     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
 547            for i in xrange(dimension) ]
 548
 549     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 550
 551
 552 def eja_ln(dimension, field=QQ):
 553     """
 554     Return the Jordan algebra corresponding to the Lorentz "ice cream"
 555     cone of the given ``dimension``.
 556
 557     EXAMPLES:
 558
 559     This multiplication table can be verified by hand::
 560
 561         sage: J = eja_ln(4)
 562         sage: e0,e1,e2,e3 = J.gens()
 563         sage: e0*e0
 564         e0
 565         sage: e0*e1
 566         e1
 567         sage: e0*e2
 568         e2
 569         sage: e0*e3
 570         e3
 571         sage: e1*e2
 572         0
 573         sage: e1*e3
 574         0
 575         sage: e2*e3
 576         0
 577
 578     In one dimension, this is the reals under multiplication::
 579
 580       sage: J1 = eja_ln(1)
 581       sage: J2 = eja_rn(1)
 582       sage: J1 == J2
 583       True
 584
 585     """
 586     Qs = []
 587     id_matrix = identity_matrix(field,dimension)
 588     for i in xrange(dimension):
 589         ei = id_matrix.column(i)
 590         Qi = zero_matrix(field,dimension)
 591         Qi.set_row(0, ei)
 592         Qi.set_column(0, ei)
 593         Qi += diagonal_matrix(dimension, [ei[0]]*dimension)
 594         # The addition of the diagonal matrix adds an extra ei[0] in the
 595         # upper-left corner of the matrix.
 596         Qi[0,0] = Qi[0,0] * ~field(2)
 597         Qs.append(Qi)
 598
 599     # The rank of the spin factor algebra is two, UNLESS we're in a
 600     # one-dimensional ambient space (the rank is bounded by the
 601     # ambient dimension).
 602     rank = min(dimension,2)
 603     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=rank)
 604
 605
 606 def eja_sn(dimension, field=QQ):
 607     """
 608     Return the simple Jordan algebra of ``dimension``-by-``dimension``
 609     symmetric matrices over ``field``.
 610
 611     EXAMPLES::
 612
 613         sage: J = eja_sn(2)
 614         sage: e0, e1, e2 = J.gens()
 615         sage: e0*e0
 616         e0
 617         sage: e1*e1
 618         e0 + e2
 619         sage: e2*e2
 620         e2
 621
 622     """
 623     S = _real_symmetric_basis(dimension, field=field)
 624     Qs = _multiplication_table_from_matrix_basis(S)
 625
 626     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 627
 628
 629 def random_eja():
 630     """
 631     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 632
 633     ALGORITHM:
 634
 635     For now, we choose a random natural number ``n`` (greater than zero)
 636     and then give you back one of the following:
 637
 638       * The cartesian product of the rational numbers ``n`` times; this is
 639         ``QQ^n`` with the Hadamard product.
 640
 641       * The Jordan spin algebra on ``QQ^n``.
 642
 643       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
 644         product.
 645
 646     Later this might be extended to return Cartesian products of the
 647     EJAs above.
 648
 649     TESTS::
 650
 651         sage: random_eja()
 652         Euclidean Jordan algebra of degree...
 653
 654     """
 655     n = ZZ.random_element(1,10).abs()
 656     constructor = choice([eja_rn, eja_ln, eja_sn])
 657     return constructor(dimension=n, field=QQ)
 658
 659
 660
 661 def _real_symmetric_basis(n, field=QQ):
 662     """
 663     Return a basis for the space of real symmetric n-by-n matrices.
 664     """
 665     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
 666     # coordinates.
 667     S = []
 668     for i in xrange(n):
 669         for j in xrange(i+1):
 670             Eij = matrix(field, n, lambda k,l: k==i and l==j)
 671             if i == j:
 672                 Sij = Eij
 673             else:
 674                 # Beware, orthogonal but not normalized!
 675                 Sij = Eij + Eij.transpose()
 676             S.append(Sij)
 677     return S
 678
 679
 680 def _multiplication_table_from_matrix_basis(basis):
 681     """
 682     At least three of the five simple Euclidean Jordan algebras have the
 683     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
 684     multiplication on the right is matrix multiplication. Given a basis
 685     for the underlying matrix space, this function returns a
 686     multiplication table (obtained by looping through the basis
 687     elements) for an algebra of those matrices.
 688     """
 689     # In S^2, for example, we nominally have four coordinates even
 690     # though the space is of dimension three only. The vector space V
 691     # is supposed to hold the entire long vector, and the subspace W
 692     # of V will be spanned by the vectors that arise from symmetric
 693     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
 694     field = basis[0].base_ring()
 695     dimension = basis[0].nrows()
 696
 697     def mat2vec(m):
 698         return vector(field, m.list())
 699
 700     def vec2mat(v):
 701         return matrix(field, dimension, v.list())
 702
 703     V = VectorSpace(field, dimension**2)
 704     W = V.span( mat2vec(s) for s in basis )
 705
 706     # Taking the span above reorders our basis (thanks, jerk!) so we
 707     # need to put our "matrix basis" in the same order as the
 708     # (reordered) vector basis.
 709     S = [ vec2mat(b) for b in W.basis() ]
 710
 711     Qs = []
 712     for s in S:
 713         # Brute force the multiplication-by-s matrix by looping
 714         # through all elements of the basis and doing the computation
 715         # to find out what the corresponding row should be. BEWARE:
 716         # these multiplication tables won't be symmetric! It therefore
 717         # becomes REALLY IMPORTANT that the underlying algebra
 718         # constructor uses ROW vectors and not COLUMN vectors. That's
 719         # why we're computing rows here and not columns.
 720         Q_rows = []
 721         for t in S:
 722             this_row = mat2vec((s*t + t*s)/2)
 723             Q_rows.append(W.coordinates(this_row))
 724         Q = matrix(field, W.dimension(), Q_rows)
 725         Qs.append(Q)
 726
 727     return Qs
 728
 729
 730 def _embed_complex_matrix(M):
 731     """
 732     Embed the n-by-n complex matrix ``M`` into the space of real
 733     matrices of size 2n-by-2n via the map the sends each entry `z = a +
 734     bi` to the block matrix ``[[a,b],[-b,a]]``.
 735
 736     EXAMPLES::
 737
 738         sage: F = QuadraticField(-1,'i')
 739         sage: x1 = F(4 - 2*i)
 740         sage: x2 = F(1 + 2*i)
 741         sage: x3 = F(-i)
 742         sage: x4 = F(6)
 743         sage: M = matrix(F,2,[x1,x2,x3,x4])
 744         sage: _embed_complex_matrix(M)
 745         [ 4  2| 1 -2]
 746         [-2  4| 2  1]
 747         [-----+-----]
 748         [ 0  1| 6  0]
 749         [-1  0| 0  6]
 750
 751     """
 752     n = M.nrows()
 753     if M.ncols() != n:
 754         raise ArgumentError("the matrix 'M' must be square")
 755     field = M.base_ring()
 756     blocks = []
 757     for z in M.list():
 758         a = z.real()
 759         b = z.imag()
 760         blocks.append(matrix(field, 2, [[a,-b],[b,a]]))
 761
 762     # We can drop the imaginaries here.
 763     return block_matrix(field.base_ring(), n, blocks)
 764
 765
 766 def _unembed_complex_matrix(M):
 767     """
 768     The inverse of _embed_complex_matrix().
 769
 770     EXAMPLES::
 771
 772         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
 773         ....:                 [-2,  1,  -4,  3],
 774         ....:                 [ 9,  10, 11, 12],
 775         ....:                 [-10, 9, -12, 11] ])
 776         sage: _unembed_complex_matrix(A)
 777         [  -2*i + 1   -4*i + 3]
 778         [ -10*i + 9 -12*i + 11]
 779     """
 780     n = ZZ(M.nrows())
 781     if M.ncols() != n:
 782         raise ArgumentError("the matrix 'M' must be square")
 783     if not n.mod(2).is_zero():
 784         raise ArgumentError("the matrix 'M' must be a complex embedding")
 785
 786     F = QuadraticField(-1, 'i')
 787     i = F.gen()
 788
 789     # Go top-left to bottom-right (reading order), converting every
 790     # 2-by-2 block we see to a single complex element.
 791     elements = []
 792     for k in xrange(n/2):
 793         for j in xrange(n/2):
 794             submat = M[2*k:2*k+2,2*j:2*j+2]
 795             if submat[0,0] != submat[1,1]:
 796                 raise ArgumentError('bad real submatrix')
 797             if submat[0,1] != -submat[1,0]:
 798                 raise ArgumentError('bad imag submatrix')
 799             z = submat[0,0] + submat[1,0]*i
 800             elements.append(z)
 801
 802     return matrix(F, n/2, elements)
 803
 804
 805 def RealSymmetricSimpleEJA(n):
 806     """
 807     The rank-n simple EJA consisting of real symmetric n-by-n
 808     matrices, the usual symmetric Jordan product, and the trace inner
 809     product. It has dimension `(n^2 + n)/2` over the reals.
 810     """
 811     pass
 812
 813 def ComplexHermitianSimpleEJA(n, field=QQ):
 814     """
 815     The rank-n simple EJA consisting of complex Hermitian n-by-n
 816     matrices over the real numbers, the usual symmetric Jordan product,
 817     and the real-part-of-trace inner product. It has dimension `n^2 over
 818     the reals.
 819     """
 820     F = QuadraticField(-1, 'i')
 821     i = F.gen()
 822     S = _real_symmetric_basis(n, field=F)
 823     T = []
 824     for s in S:
 825         T.append(s)
 826         T.append(i*s)
 827     embed_T = [ _embed_complex_matrix(t) for t in T ]
 828     Qs = _multiplication_table_from_matrix_basis(embed_T)
 829     return FiniteDimensionalEuclideanJordanAlgebra(field, Qs, rank=n)
 830
 831 def QuaternionHermitianSimpleEJA(n):
 832     """
 833     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
 834     matrices, the usual symmetric Jordan product, and the
 835     real-part-of-trace inner product. It has dimension `2n^2 - n` over
 836     the reals.
 837     """
 838     pass
 839
 840 def OctonionHermitianSimpleEJA(n):
 841     """
 842     This shit be crazy. It has dimension 27 over the reals.
 843     """
 844     n = 3
 845     pass
 846
 847 def JordanSpinSimpleEJA(n):
 848     """
 849     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
 850     with the usual inner product and jordan product ``x*y =
 851     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
 852     the reals.
 853     """
 854     pass