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):
 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 = JordanSpinSimpleEJA(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 = JordanSpinSimpleEJA(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 = JordanSpinSimpleEJA(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
 553 def random_eja():
 554     """
 555     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 556
 557     ALGORITHM:
 558
 559     For now, we choose a random natural number ``n`` (greater than zero)
 560     and then give you back one of the following:
 561
 562       * The cartesian product of the rational numbers ``n`` times; this is
 563         ``QQ^n`` with the Hadamard product.
 564
 565       * The Jordan spin algebra on ``QQ^n``.
 566
 567       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
 568         product.
 569
 570     Later this might be extended to return Cartesian products of the
 571     EJAs above.
 572
 573     TESTS::
 574
 575         sage: random_eja()
 576         Euclidean Jordan algebra of degree...
 577
 578     """
 579     n = ZZ.random_element(1,5).abs()
 580     constructor = choice([eja_rn,
 581                           JordanSpinSimpleEJA,
 582                           RealSymmetricSimpleEJA,
 583                           ComplexHermitianSimpleEJA])
 584     return constructor(n, field=QQ)
 585
 586
 587
 588 def _real_symmetric_basis(n, field=QQ):
 589     """
 590     Return a basis for the space of real symmetric n-by-n matrices.
 591     """
 592     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
 593     # coordinates.
 594     S = []
 595     for i in xrange(n):
 596         for j in xrange(i+1):
 597             Eij = matrix(field, n, lambda k,l: k==i and l==j)
 598             if i == j:
 599                 Sij = Eij
 600             else:
 601                 # Beware, orthogonal but not normalized!
 602                 Sij = Eij + Eij.transpose()
 603             S.append(Sij)
 604     return S
 605
 606
 607 def _multiplication_table_from_matrix_basis(basis):
 608     """
 609     At least three of the five simple Euclidean Jordan algebras have the
 610     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
 611     multiplication on the right is matrix multiplication. Given a basis
 612     for the underlying matrix space, this function returns a
 613     multiplication table (obtained by looping through the basis
 614     elements) for an algebra of those matrices.
 615     """
 616     # In S^2, for example, we nominally have four coordinates even
 617     # though the space is of dimension three only. The vector space V
 618     # is supposed to hold the entire long vector, and the subspace W
 619     # of V will be spanned by the vectors that arise from symmetric
 620     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
 621     field = basis[0].base_ring()
 622     dimension = basis[0].nrows()
 623
 624     def mat2vec(m):
 625         return vector(field, m.list())
 626
 627     def vec2mat(v):
 628         return matrix(field, dimension, v.list())
 629
 630     V = VectorSpace(field, dimension**2)
 631     W = V.span( mat2vec(s) for s in basis )
 632
 633     # Taking the span above reorders our basis (thanks, jerk!) so we
 634     # need to put our "matrix basis" in the same order as the
 635     # (reordered) vector basis.
 636     S = [ vec2mat(b) for b in W.basis() ]
 637
 638     Qs = []
 639     for s in S:
 640         # Brute force the multiplication-by-s matrix by looping
 641         # through all elements of the basis and doing the computation
 642         # to find out what the corresponding row should be. BEWARE:
 643         # these multiplication tables won't be symmetric! It therefore
 644         # becomes REALLY IMPORTANT that the underlying algebra
 645         # constructor uses ROW vectors and not COLUMN vectors. That's
 646         # why we're computing rows here and not columns.
 647         Q_rows = []
 648         for t in S:
 649             this_row = mat2vec((s*t + t*s)/2)
 650             Q_rows.append(W.coordinates(this_row))
 651         Q = matrix(field, W.dimension(), Q_rows)
 652         Qs.append(Q)
 653
 654     return Qs
 655
 656
 657 def _embed_complex_matrix(M):
 658     """
 659     Embed the n-by-n complex matrix ``M`` into the space of real
 660     matrices of size 2n-by-2n via the map the sends each entry `z = a +
 661     bi` to the block matrix ``[[a,b],[-b,a]]``.
 662
 663     EXAMPLES::
 664
 665         sage: F = QuadraticField(-1,'i')
 666         sage: x1 = F(4 - 2*i)
 667         sage: x2 = F(1 + 2*i)
 668         sage: x3 = F(-i)
 669         sage: x4 = F(6)
 670         sage: M = matrix(F,2,[x1,x2,x3,x4])
 671         sage: _embed_complex_matrix(M)
 672         [ 4  2| 1 -2]
 673         [-2  4| 2  1]
 674         [-----+-----]
 675         [ 0  1| 6  0]
 676         [-1  0| 0  6]
 677
 678     """
 679     n = M.nrows()
 680     if M.ncols() != n:
 681         raise ArgumentError("the matrix 'M' must be square")
 682     field = M.base_ring()
 683     blocks = []
 684     for z in M.list():
 685         a = z.real()
 686         b = z.imag()
 687         blocks.append(matrix(field, 2, [[a,-b],[b,a]]))
 688
 689     # We can drop the imaginaries here.
 690     return block_matrix(field.base_ring(), n, blocks)
 691
 692
 693 def _unembed_complex_matrix(M):
 694     """
 695     The inverse of _embed_complex_matrix().
 696
 697     EXAMPLES::
 698
 699         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
 700         ....:                 [-2,  1,  -4,  3],
 701         ....:                 [ 9,  10, 11, 12],
 702         ....:                 [-10, 9, -12, 11] ])
 703         sage: _unembed_complex_matrix(A)
 704         [  -2*i + 1   -4*i + 3]
 705         [ -10*i + 9 -12*i + 11]
 706     """
 707     n = ZZ(M.nrows())
 708     if M.ncols() != n:
 709         raise ArgumentError("the matrix 'M' must be square")
 710     if not n.mod(2).is_zero():
 711         raise ArgumentError("the matrix 'M' must be a complex embedding")
 712
 713     F = QuadraticField(-1, 'i')
 714     i = F.gen()
 715
 716     # Go top-left to bottom-right (reading order), converting every
 717     # 2-by-2 block we see to a single complex element.
 718     elements = []
 719     for k in xrange(n/2):
 720         for j in xrange(n/2):
 721             submat = M[2*k:2*k+2,2*j:2*j+2]
 722             if submat[0,0] != submat[1,1]:
 723                 raise ArgumentError('bad real submatrix')
 724             if submat[0,1] != -submat[1,0]:
 725                 raise ArgumentError('bad imag submatrix')
 726             z = submat[0,0] + submat[1,0]*i
 727             elements.append(z)
 728
 729     return matrix(F, n/2, elements)
 730
 731
 732 def RealSymmetricSimpleEJA(n, field=QQ):
 733     """
 734     The rank-n simple EJA consisting of real symmetric n-by-n
 735     matrices, the usual symmetric Jordan product, and the trace inner
 736     product. It has dimension `(n^2 + n)/2` over the reals.
 737
 738     EXAMPLES::
 739
 740         sage: J = RealSymmetricSimpleEJA(2)
 741         sage: e0, e1, e2 = J.gens()
 742         sage: e0*e0
 743         e0
 744         sage: e1*e1
 745         e0 + e2
 746         sage: e2*e2
 747         e2
 748     """
 749     S = _real_symmetric_basis(n, field=field)
 750     Qs = _multiplication_table_from_matrix_basis(S)
 751
 752     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=n)
 753
 754
 755 def ComplexHermitianSimpleEJA(n, field=QQ):
 756     """
 757     The rank-n simple EJA consisting of complex Hermitian n-by-n
 758     matrices over the real numbers, the usual symmetric Jordan product,
 759     and the real-part-of-trace inner product. It has dimension `n^2 over
 760     the reals.
 761     """
 762     F = QuadraticField(-1, 'i')
 763     i = F.gen()
 764     S = _real_symmetric_basis(n, field=F)
 765     T = []
 766     for s in S:
 767         T.append(s)
 768         T.append(i*s)
 769     embed_T = [ _embed_complex_matrix(t) for t in T ]
 770     Qs = _multiplication_table_from_matrix_basis(embed_T)
 771     return FiniteDimensionalEuclideanJordanAlgebra(field, Qs, rank=n)
 772
 773 def QuaternionHermitianSimpleEJA(n):
 774     """
 775     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
 776     matrices, the usual symmetric Jordan product, and the
 777     real-part-of-trace inner product. It has dimension `2n^2 - n` over
 778     the reals.
 779     """
 780     pass
 781
 782 def OctonionHermitianSimpleEJA(n):
 783     """
 784     This shit be crazy. It has dimension 27 over the reals.
 785     """
 786     n = 3
 787     pass
 788
 789 def JordanSpinSimpleEJA(n, field=QQ):
 790     """
 791     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
 792     with the usual inner product and jordan product ``x*y =
 793     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
 794     the reals.
 795
 796     EXAMPLES:
 797
 798     This multiplication table can be verified by hand::
 799
 800         sage: J = JordanSpinSimpleEJA(4)
 801         sage: e0,e1,e2,e3 = J.gens()
 802         sage: e0*e0
 803         e0
 804         sage: e0*e1
 805         e1
 806         sage: e0*e2
 807         e2
 808         sage: e0*e3
 809         e3
 810         sage: e1*e2
 811         0
 812         sage: e1*e3
 813         0
 814         sage: e2*e3
 815         0
 816
 817     In one dimension, this is the reals under multiplication::
 818
 819       sage: J1 = JordanSpinSimpleEJA(1)
 820       sage: J2 = eja_rn(1)
 821       sage: J1 == J2
 822       True
 823
 824     """
 825     Qs = []
 826     id_matrix = identity_matrix(field, n)
 827     for i in xrange(n):
 828         ei = id_matrix.column(i)
 829         Qi = zero_matrix(field, n)
 830         Qi.set_row(0, ei)
 831         Qi.set_column(0, ei)
 832         Qi += diagonal_matrix(n, [ei[0]]*n)
 833         # The addition of the diagonal matrix adds an extra ei[0] in the
 834         # upper-left corner of the matrix.
 835         Qi[0,0] = Qi[0,0] * ~field(2)
 836         Qs.append(Qi)
 837
 838     # The rank of the spin factor algebra is two, UNLESS we're in a
 839     # one-dimensional ambient space (the rank is bounded by the
 840     # ambient dimension).
 841     return FiniteDimensionalEuclideanJordanAlgebra(field, Qs, rank=min(n,2))