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 _complex_hermitian_basis(n, field=QQ):
 608     """
 609     Returns a basis for the space of complex Hermitian n-by-n matrices.
 610
 611     TESTS::
 612
 613         sage: set_random_seed()
 614         sage: n = ZZ.random_element(1,5).abs()
 615         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
 616         True
 617
 618     """
 619     F = QuadraticField(-1, 'I')
 620     I = F.gen()
 621
 622     # This is like the symmetric case, but we need to be careful:
 623     #
 624     #   * We want conjugate-symmetry, not just symmetry.
 625     #   * The diagonal will (as a result) be real.
 626     #
 627     S = []
 628     for i in xrange(n):
 629         for j in xrange(i+1):
 630             Eij = matrix(field, n, lambda k,l: k==i and l==j)
 631             if i == j:
 632                 Sij = _embed_complex_matrix(Eij)
 633                 S.append(Sij)
 634             else:
 635                 # Beware, orthogonal but not normalized! The second one
 636                 # has a minus because it's conjugated.
 637                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
 638                 S.append(Sij_real)
 639                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
 640                 S.append(Sij_imag)
 641     return S
 642
 643
 644 def _multiplication_table_from_matrix_basis(basis):
 645     """
 646     At least three of the five simple Euclidean Jordan algebras have the
 647     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
 648     multiplication on the right is matrix multiplication. Given a basis
 649     for the underlying matrix space, this function returns a
 650     multiplication table (obtained by looping through the basis
 651     elements) for an algebra of those matrices.
 652     """
 653     # In S^2, for example, we nominally have four coordinates even
 654     # though the space is of dimension three only. The vector space V
 655     # is supposed to hold the entire long vector, and the subspace W
 656     # of V will be spanned by the vectors that arise from symmetric
 657     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
 658     field = basis[0].base_ring()
 659     dimension = basis[0].nrows()
 660
 661     def mat2vec(m):
 662         return vector(field, m.list())
 663
 664     def vec2mat(v):
 665         return matrix(field, dimension, v.list())
 666
 667     V = VectorSpace(field, dimension**2)
 668     W = V.span( mat2vec(s) for s in basis )
 669
 670     # Taking the span above reorders our basis (thanks, jerk!) so we
 671     # need to put our "matrix basis" in the same order as the
 672     # (reordered) vector basis.
 673     S = [ vec2mat(b) for b in W.basis() ]
 674
 675     Qs = []
 676     for s in S:
 677         # Brute force the multiplication-by-s matrix by looping
 678         # through all elements of the basis and doing the computation
 679         # to find out what the corresponding row should be. BEWARE:
 680         # these multiplication tables won't be symmetric! It therefore
 681         # becomes REALLY IMPORTANT that the underlying algebra
 682         # constructor uses ROW vectors and not COLUMN vectors. That's
 683         # why we're computing rows here and not columns.
 684         Q_rows = []
 685         for t in S:
 686             this_row = mat2vec((s*t + t*s)/2)
 687             Q_rows.append(W.coordinates(this_row))
 688         Q = matrix(field, W.dimension(), Q_rows)
 689         Qs.append(Q)
 690
 691     return Qs
 692
 693
 694 def _embed_complex_matrix(M):
 695     """
 696     Embed the n-by-n complex matrix ``M`` into the space of real
 697     matrices of size 2n-by-2n via the map the sends each entry `z = a +
 698     bi` to the block matrix ``[[a,b],[-b,a]]``.
 699
 700     EXAMPLES::
 701
 702         sage: F = QuadraticField(-1,'i')
 703         sage: x1 = F(4 - 2*i)
 704         sage: x2 = F(1 + 2*i)
 705         sage: x3 = F(-i)
 706         sage: x4 = F(6)
 707         sage: M = matrix(F,2,[x1,x2,x3,x4])
 708         sage: _embed_complex_matrix(M)
 709         [ 4  2| 1 -2]
 710         [-2  4| 2  1]
 711         [-----+-----]
 712         [ 0  1| 6  0]
 713         [-1  0| 0  6]
 714
 715     """
 716     n = M.nrows()
 717     if M.ncols() != n:
 718         raise ArgumentError("the matrix 'M' must be square")
 719     field = M.base_ring()
 720     blocks = []
 721     for z in M.list():
 722         a = z.real()
 723         b = z.imag()
 724         blocks.append(matrix(field, 2, [[a,-b],[b,a]]))
 725
 726     # We can drop the imaginaries here.
 727     return block_matrix(field.base_ring(), n, blocks)
 728
 729
 730 def _unembed_complex_matrix(M):
 731     """
 732     The inverse of _embed_complex_matrix().
 733
 734     EXAMPLES::
 735
 736         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
 737         ....:                 [-2,  1,  -4,  3],
 738         ....:                 [ 9,  10, 11, 12],
 739         ....:                 [-10, 9, -12, 11] ])
 740         sage: _unembed_complex_matrix(A)
 741         [  -2*i + 1   -4*i + 3]
 742         [ -10*i + 9 -12*i + 11]
 743     """
 744     n = ZZ(M.nrows())
 745     if M.ncols() != n:
 746         raise ArgumentError("the matrix 'M' must be square")
 747     if not n.mod(2).is_zero():
 748         raise ArgumentError("the matrix 'M' must be a complex embedding")
 749
 750     F = QuadraticField(-1, 'i')
 751     i = F.gen()
 752
 753     # Go top-left to bottom-right (reading order), converting every
 754     # 2-by-2 block we see to a single complex element.
 755     elements = []
 756     for k in xrange(n/2):
 757         for j in xrange(n/2):
 758             submat = M[2*k:2*k+2,2*j:2*j+2]
 759             if submat[0,0] != submat[1,1]:
 760                 raise ArgumentError('bad real submatrix')
 761             if submat[0,1] != -submat[1,0]:
 762                 raise ArgumentError('bad imag submatrix')
 763             z = submat[0,0] + submat[1,0]*i
 764             elements.append(z)
 765
 766     return matrix(F, n/2, elements)
 767
 768
 769 def RealSymmetricSimpleEJA(n, field=QQ):
 770     """
 771     The rank-n simple EJA consisting of real symmetric n-by-n
 772     matrices, the usual symmetric Jordan product, and the trace inner
 773     product. It has dimension `(n^2 + n)/2` over the reals.
 774
 775     EXAMPLES::
 776
 777         sage: J = RealSymmetricSimpleEJA(2)
 778         sage: e0, e1, e2 = J.gens()
 779         sage: e0*e0
 780         e0
 781         sage: e1*e1
 782         e0 + e2
 783         sage: e2*e2
 784         e2
 785
 786     TESTS:
 787
 788     The degree of this algebra is `(n^2 + n) / 2`::
 789
 790         sage: set_random_seed()
 791         sage: n = ZZ.random_element(1,5).abs()
 792         sage: J = RealSymmetricSimpleEJA(n)
 793         sage: J.degree() == (n^2 + n)/2
 794         True
 795
 796     """
 797     S = _real_symmetric_basis(n, field=field)
 798     Qs = _multiplication_table_from_matrix_basis(S)
 799
 800     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=n)
 801
 802
 803 def ComplexHermitianSimpleEJA(n, field=QQ):
 804     """
 805     The rank-n simple EJA consisting of complex Hermitian n-by-n
 806     matrices over the real numbers, the usual symmetric Jordan product,
 807     and the real-part-of-trace inner product. It has dimension `n^2` over
 808     the reals.
 809
 810     TESTS:
 811
 812     The degree of this algebra is `n^2`::
 813
 814         sage: set_random_seed()
 815         sage: n = ZZ.random_element(1,5).abs()
 816         sage: J = ComplexHermitianSimpleEJA(n)
 817         sage: J.degree() == n^2
 818         True
 819
 820     """
 821     S = _complex_hermitian_basis(n)
 822     Qs = _multiplication_table_from_matrix_basis(S)
 823     return FiniteDimensionalEuclideanJordanAlgebra(field, Qs, rank=n)
 824
 825
 826 def QuaternionHermitianSimpleEJA(n):
 827     """
 828     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
 829     matrices, the usual symmetric Jordan product, and the
 830     real-part-of-trace inner product. It has dimension `2n^2 - n` over
 831     the reals.
 832     """
 833     pass
 834
 835 def OctonionHermitianSimpleEJA(n):
 836     """
 837     This shit be crazy. It has dimension 27 over the reals.
 838     """
 839     n = 3
 840     pass
 841
 842 def JordanSpinSimpleEJA(n, field=QQ):
 843     """
 844     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
 845     with the usual inner product and jordan product ``x*y =
 846     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
 847     the reals.
 848
 849     EXAMPLES:
 850
 851     This multiplication table can be verified by hand::
 852
 853         sage: J = JordanSpinSimpleEJA(4)
 854         sage: e0,e1,e2,e3 = J.gens()
 855         sage: e0*e0
 856         e0
 857         sage: e0*e1
 858         e1
 859         sage: e0*e2
 860         e2
 861         sage: e0*e3
 862         e3
 863         sage: e1*e2
 864         0
 865         sage: e1*e3
 866         0
 867         sage: e2*e3
 868         0
 869
 870     In one dimension, this is the reals under multiplication::
 871
 872       sage: J1 = JordanSpinSimpleEJA(1)
 873       sage: J2 = eja_rn(1)
 874       sage: J1 == J2
 875       True
 876
 877     """
 878     Qs = []
 879     id_matrix = identity_matrix(field, n)
 880     for i in xrange(n):
 881         ei = id_matrix.column(i)
 882         Qi = zero_matrix(field, n)
 883         Qi.set_row(0, ei)
 884         Qi.set_column(0, ei)
 885         Qi += diagonal_matrix(n, [ei[0]]*n)
 886         # The addition of the diagonal matrix adds an extra ei[0] in the
 887         # upper-left corner of the matrix.
 888         Qi[0,0] = Qi[0,0] * ~field(2)
 889         Qs.append(Qi)
 890
 891     # The rank of the spin factor algebra is two, UNLESS we're in a
 892     # one-dimensional ambient space (the rank is bounded by the
 893     # ambient dimension).
 894     return FiniteDimensionalEuclideanJordanAlgebra(field, Qs, rank=min(n,2))