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                               natural_basis=None,
  25                               inner_product=None):
  26         n = len(mult_table)
  27         mult_table = [b.base_extend(field) for b in mult_table]
  28         for b in mult_table:
  29             b.set_immutable()
  30             if not (is_Matrix(b) and b.dimensions() == (n, n)):
  31                 raise ValueError("input is not a multiplication table")
  32         mult_table = tuple(mult_table)
  33
  34         cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
  35         cat.or_subcategory(category)
  36         if assume_associative:
  37             cat = cat.Associative()
  38
  39         names = normalize_names(n, names)
  40
  41         fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls)
  42         return fda.__classcall__(cls,
  43                                  field,
  44                                  mult_table,
  45                                  assume_associative=assume_associative,
  46                                  names=names,
  47                                  category=cat,
  48                                  rank=rank,
  49                                  natural_basis=natural_basis,
  50                                  inner_product=inner_product)
  51
  52
  53     def __init__(self, field,
  54                  mult_table,
  55                  names='e',
  56                  assume_associative=False,
  57                  category=None,
  58                  rank=None,
  59                  natural_basis=None,
  60                  inner_product=None):
  61         """
  62         EXAMPLES:
  63
  64         By definition, Jordan multiplication commutes::
  65
  66             sage: set_random_seed()
  67             sage: J = random_eja()
  68             sage: x = J.random_element()
  69             sage: y = J.random_element()
  70             sage: x*y == y*x
  71             True
  72
  73         """
  74         self._rank = rank
  75         self._natural_basis = natural_basis
  76         self._inner_product = inner_product
  77         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
  78         fda.__init__(field,
  79                      mult_table,
  80                      names=names,
  81                      category=category)
  82
  83
  84     def _repr_(self):
  85         """
  86         Return a string representation of ``self``.
  87         """
  88         fmt = "Euclidean Jordan algebra of degree {} over {}"
  89         return fmt.format(self.degree(), self.base_ring())
  90
  91
  92     def inner_product(self, x, y):
  93         """
  94         The inner product associated with this Euclidean Jordan algebra.
  95
  96         Will default to the trace inner product if nothing else.
  97
  98         EXAMPLES:
  99
 100         The inner product must satisfy its axiom for this algebra to truly
 101         be a Euclidean Jordan Algebra::
 102
 103             sage: set_random_seed()
 104             sage: J = random_eja()
 105             sage: x = J.random_element()
 106             sage: y = J.random_element()
 107             sage: z = J.random_element()
 108             sage: (x*y).inner_product(z) == y.inner_product(x*z)
 109             True
 110
 111         """
 112         if (not x in self) or (not y in self):
 113             raise TypeError("arguments must live in this algebra")
 114         if self._inner_product is None:
 115             return x.trace_inner_product(y)
 116         else:
 117             return self._inner_product(x,y)
 118
 119
 120     def natural_basis(self):
 121         """
 122         Return a more-natural representation of this algebra's basis.
 123
 124         Every finite-dimensional Euclidean Jordan Algebra is a direct
 125         sum of five simple algebras, four of which comprise Hermitian
 126         matrices. This method returns the original "natural" basis
 127         for our underlying vector space. (Typically, the natural basis
 128         is used to construct the multiplication table in the first place.)
 129
 130         Note that this will always return a matrix. The standard basis
 131         in `R^n` will be returned as `n`-by-`1` column matrices.
 132
 133         EXAMPLES::
 134
 135             sage: J = RealSymmetricSimpleEJA(2)
 136             sage: J.basis()
 137             Family (e0, e1, e2)
 138             sage: J.natural_basis()
 139             (
 140             [1 0]  [0 1]  [0 0]
 141             [0 0], [1 0], [0 1]
 142             )
 143
 144         ::
 145
 146             sage: J = JordanSpinSimpleEJA(2)
 147             sage: J.basis()
 148             Family (e0, e1)
 149             sage: J.natural_basis()
 150             (
 151             [1]  [0]
 152             [0], [1]
 153             )
 154
 155         """
 156         if self._natural_basis is None:
 157             return tuple( b.vector().column() for b in self.basis() )
 158         else:
 159             return self._natural_basis
 160
 161
 162     def rank(self):
 163         """
 164         Return the rank of this EJA.
 165         """
 166         if self._rank is None:
 167             raise ValueError("no rank specified at genesis")
 168         else:
 169             return self._rank
 170
 171
 172     class Element(FiniteDimensionalAlgebraElement):
 173         """
 174         An element of a Euclidean Jordan algebra.
 175         """
 176
 177         def __init__(self, A, elt=None):
 178             """
 179             EXAMPLES:
 180
 181             The identity in `S^n` is converted to the identity in the EJA::
 182
 183                 sage: J = RealSymmetricSimpleEJA(3)
 184                 sage: I = identity_matrix(QQ,3)
 185                 sage: J(I) == J.one()
 186                 True
 187
 188             This skew-symmetric matrix can't be represented in the EJA::
 189
 190                 sage: J = RealSymmetricSimpleEJA(3)
 191                 sage: A = matrix(QQ,3, lambda i,j: i-j)
 192                 sage: J(A)
 193                 Traceback (most recent call last):
 194                 ...
 195                 ArithmeticError: vector is not in free module
 196
 197             """
 198             # Goal: if we're given a matrix, and if it lives in our
 199             # parent algebra's "natural ambient space," convert it
 200             # into an algebra element.
 201             #
 202             # The catch is, we make a recursive call after converting
 203             # the given matrix into a vector that lives in the algebra.
 204             # This we need to try the parent class initializer first,
 205             # to avoid recursing forever if we're given something that
 206             # already fits into the algebra, but also happens to live
 207             # in the parent's "natural ambient space" (this happens with
 208             # vectors in R^n).
 209             try:
 210                 FiniteDimensionalAlgebraElement.__init__(self, A, elt)
 211             except ValueError:
 212                 natural_basis = A.natural_basis()
 213                 if elt in natural_basis[0].matrix_space():
 214                     # Thanks for nothing! Matrix spaces aren't vector
 215                     # spaces in Sage, so we have to figure out its
 216                     # natural-basis coordinates ourselves.
 217                     V = VectorSpace(elt.base_ring(), elt.nrows()**2)
 218                     W = V.span( _mat2vec(s) for s in natural_basis )
 219                     coords =  W.coordinates(_mat2vec(elt))
 220                     FiniteDimensionalAlgebraElement.__init__(self, A, coords)
 221
 222         def __pow__(self, n):
 223             """
 224             Return ``self`` raised to the power ``n``.
 225
 226             Jordan algebras are always power-associative; see for
 227             example Faraut and Koranyi, Proposition II.1.2 (ii).
 228
 229             .. WARNING:
 230
 231                 We have to override this because our superclass uses row vectors
 232                 instead of column vectors! We, on the other hand, assume column
 233                 vectors everywhere.
 234
 235             EXAMPLES::
 236
 237                 sage: set_random_seed()
 238                 sage: x = random_eja().random_element()
 239                 sage: x.operator_matrix()*x.vector() == (x^2).vector()
 240                 True
 241
 242             A few examples of power-associativity::
 243
 244                 sage: set_random_seed()
 245                 sage: x = random_eja().random_element()
 246                 sage: x*(x*x)*(x*x) == x^5
 247                 True
 248                 sage: (x*x)*(x*x*x) == x^5
 249                 True
 250
 251             We also know that powers operator-commute (Koecher, Chapter
 252             III, Corollary 1)::
 253
 254                 sage: set_random_seed()
 255                 sage: x = random_eja().random_element()
 256                 sage: m = ZZ.random_element(0,10)
 257                 sage: n = ZZ.random_element(0,10)
 258                 sage: Lxm = (x^m).operator_matrix()
 259                 sage: Lxn = (x^n).operator_matrix()
 260                 sage: Lxm*Lxn == Lxn*Lxm
 261                 True
 262
 263             """
 264             A = self.parent()
 265             if n == 0:
 266                 return A.one()
 267             elif n == 1:
 268                 return self
 269             else:
 270                 return A( (self.operator_matrix()**(n-1))*self.vector() )
 271
 272
 273         def characteristic_polynomial(self):
 274             """
 275             Return my characteristic polynomial (if I'm a regular
 276             element).
 277
 278             Eventually this should be implemented in terms of the parent
 279             algebra's characteristic polynomial that works for ALL
 280             elements.
 281             """
 282             if self.is_regular():
 283                 return self.minimal_polynomial()
 284             else:
 285                 raise NotImplementedError('irregular element')
 286
 287
 288         def inner_product(self, other):
 289             """
 290             Return the parent algebra's inner product of myself and ``other``.
 291
 292             EXAMPLES:
 293
 294             The inner product in the Jordan spin algebra is the usual
 295             inner product on `R^n` (this example only works because the
 296             basis for the Jordan algebra is the standard basis in `R^n`)::
 297
 298                 sage: J = JordanSpinSimpleEJA(3)
 299                 sage: x = vector(QQ,[1,2,3])
 300                 sage: y = vector(QQ,[4,5,6])
 301                 sage: x.inner_product(y)
 302                 32
 303                 sage: J(x).inner_product(J(y))
 304                 32
 305
 306             The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
 307             multiplication is the usual matrix multiplication in `S^n`,
 308             so the inner product of the identity matrix with itself
 309             should be the `n`::
 310
 311                 sage: J = RealSymmetricSimpleEJA(3)
 312                 sage: J.one().inner_product(J.one())
 313                 3
 314
 315             Likewise, the inner product on `C^n` is `<X,Y> =
 316             Re(trace(X*Y))`, where we must necessarily take the real
 317             part because the product of Hermitian matrices may not be
 318             Hermitian::
 319
 320                 sage: J = ComplexHermitianSimpleEJA(3)
 321                 sage: J.one().inner_product(J.one())
 322                 3
 323
 324             Ditto for the quaternions::
 325
 326                 sage: J = QuaternionHermitianSimpleEJA(3)
 327                 sage: J.one().inner_product(J.one())
 328                 3
 329
 330             TESTS:
 331
 332             Ensure that we can always compute an inner product, and that
 333             it gives us back a real number::
 334
 335                 sage: set_random_seed()
 336                 sage: J = random_eja()
 337                 sage: x = J.random_element()
 338                 sage: y = J.random_element()
 339                 sage: x.inner_product(y) in RR
 340                 True
 341
 342             """
 343             P = self.parent()
 344             if not other in P:
 345                 raise TypeError("'other' must live in the same algebra")
 346
 347             return P.inner_product(self, other)
 348
 349
 350         def operator_commutes_with(self, other):
 351             """
 352             Return whether or not this element operator-commutes
 353             with ``other``.
 354
 355             EXAMPLES:
 356
 357             The definition of a Jordan algebra says that any element
 358             operator-commutes with its square::
 359
 360                 sage: set_random_seed()
 361                 sage: x = random_eja().random_element()
 362                 sage: x.operator_commutes_with(x^2)
 363                 True
 364
 365             TESTS:
 366
 367             Test Lemma 1 from Chapter III of Koecher::
 368
 369                 sage: set_random_seed()
 370                 sage: J = random_eja()
 371                 sage: u = J.random_element()
 372                 sage: v = J.random_element()
 373                 sage: lhs = u.operator_commutes_with(u*v)
 374                 sage: rhs = v.operator_commutes_with(u^2)
 375                 sage: lhs == rhs
 376                 True
 377
 378             """
 379             if not other in self.parent():
 380                 raise TypeError("'other' must live in the same algebra")
 381
 382             A = self.operator_matrix()
 383             B = other.operator_matrix()
 384             return (A*B == B*A)
 385
 386
 387         def det(self):
 388             """
 389             Return my determinant, the product of my eigenvalues.
 390
 391             EXAMPLES::
 392
 393                 sage: J = JordanSpinSimpleEJA(2)
 394                 sage: e0,e1 = J.gens()
 395                 sage: x = e0 + e1
 396                 sage: x.det()
 397                 0
 398                 sage: J = JordanSpinSimpleEJA(3)
 399                 sage: e0,e1,e2 = J.gens()
 400                 sage: x = e0 + e1 + e2
 401                 sage: x.det()
 402                 -1
 403
 404             """
 405             cs = self.characteristic_polynomial().coefficients(sparse=False)
 406             r = len(cs) - 1
 407             if r >= 0:
 408                 return cs[0] * (-1)**r
 409             else:
 410                 raise ValueError('charpoly had no coefficients')
 411
 412
 413         def inverse(self):
 414             """
 415             Return the Jordan-multiplicative inverse of this element.
 416
 417             We can't use the superclass method because it relies on the
 418             algebra being associative.
 419
 420             EXAMPLES:
 421
 422             The inverse in the spin factor algebra is given in Alizadeh's
 423             Example 11.11::
 424
 425                 sage: set_random_seed()
 426                 sage: n = ZZ.random_element(1,10)
 427                 sage: J = JordanSpinSimpleEJA(n)
 428                 sage: x = J.random_element()
 429                 sage: while x.is_zero():
 430                 ....:     x = J.random_element()
 431                 sage: x_vec = x.vector()
 432                 sage: x0 = x_vec[0]
 433                 sage: x_bar = x_vec[1:]
 434                 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
 435                 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
 436                 sage: x_inverse = coeff*inv_vec
 437                 sage: x.inverse() == J(x_inverse)
 438                 True
 439
 440             TESTS:
 441
 442             The identity element is its own inverse::
 443
 444                 sage: set_random_seed()
 445                 sage: J = random_eja()
 446                 sage: J.one().inverse() == J.one()
 447                 True
 448
 449             If an element has an inverse, it acts like one. TODO: this
 450             can be a lot less ugly once ``is_invertible`` doesn't crash
 451             on irregular elements::
 452
 453                 sage: set_random_seed()
 454                 sage: J = random_eja()
 455                 sage: x = J.random_element()
 456                 sage: try:
 457                 ....:     x.inverse()*x == J.one()
 458                 ....: except:
 459                 ....:     True
 460                 True
 461
 462             """
 463             if self.parent().is_associative():
 464                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 465                 return elt.inverse()
 466
 467             # TODO: we can do better once the call to is_invertible()
 468             # doesn't crash on irregular elements.
 469             #if not self.is_invertible():
 470             #    raise ValueError('element is not invertible')
 471
 472             # We do this a little different than the usual recursive
 473             # call to a finite-dimensional algebra element, because we
 474             # wind up with an inverse that lives in the subalgebra and
 475             # we need information about the parent to convert it back.
 476             V = self.span_of_powers()
 477             assoc_subalg = self.subalgebra_generated_by()
 478             # Mis-design warning: the basis used for span_of_powers()
 479             # and subalgebra_generated_by() must be the same, and in
 480             # the same order!
 481             elt = assoc_subalg(V.coordinates(self.vector()))
 482
 483             # This will be in the subalgebra's coordinates...
 484             fda_elt = FiniteDimensionalAlgebraElement(assoc_subalg, elt)
 485             subalg_inverse = fda_elt.inverse()
 486
 487             # So we have to convert back...
 488             basis = [ self.parent(v) for v in V.basis() ]
 489             pairs = zip(subalg_inverse.vector(), basis)
 490             return self.parent().linear_combination(pairs)
 491
 492
 493         def is_invertible(self):
 494             """
 495             Return whether or not this element is invertible.
 496
 497             We can't use the superclass method because it relies on
 498             the algebra being associative.
 499             """
 500             return not self.det().is_zero()
 501
 502
 503         def is_nilpotent(self):
 504             """
 505             Return whether or not some power of this element is zero.
 506
 507             The superclass method won't work unless we're in an
 508             associative algebra, and we aren't. However, we generate
 509             an assocoative subalgebra and we're nilpotent there if and
 510             only if we're nilpotent here (probably).
 511
 512             TESTS:
 513
 514             The identity element is never nilpotent::
 515
 516                 sage: set_random_seed()
 517                 sage: random_eja().one().is_nilpotent()
 518                 False
 519
 520             The additive identity is always nilpotent::
 521
 522                 sage: set_random_seed()
 523                 sage: random_eja().zero().is_nilpotent()
 524                 True
 525
 526             """
 527             # The element we're going to call "is_nilpotent()" on.
 528             # Either myself, interpreted as an element of a finite-
 529             # dimensional algebra, or an element of an associative
 530             # subalgebra.
 531             elt = None
 532
 533             if self.parent().is_associative():
 534                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 535             else:
 536                 V = self.span_of_powers()
 537                 assoc_subalg = self.subalgebra_generated_by()
 538                 # Mis-design warning: the basis used for span_of_powers()
 539                 # and subalgebra_generated_by() must be the same, and in
 540                 # the same order!
 541                 elt = assoc_subalg(V.coordinates(self.vector()))
 542
 543             # Recursive call, but should work since elt lives in an
 544             # associative algebra.
 545             return elt.is_nilpotent()
 546
 547
 548         def is_regular(self):
 549             """
 550             Return whether or not this is a regular element.
 551
 552             EXAMPLES:
 553
 554             The identity element always has degree one, but any element
 555             linearly-independent from it is regular::
 556
 557                 sage: J = JordanSpinSimpleEJA(5)
 558                 sage: J.one().is_regular()
 559                 False
 560                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 561                 sage: for x in J.gens():
 562                 ....:     (J.one() + x).is_regular()
 563                 False
 564                 True
 565                 True
 566                 True
 567                 True
 568
 569             """
 570             return self.degree() == self.parent().rank()
 571
 572
 573         def degree(self):
 574             """
 575             Compute the degree of this element the straightforward way
 576             according to the definition; by appending powers to a list
 577             and figuring out its dimension (that is, whether or not
 578             they're linearly dependent).
 579
 580             EXAMPLES::
 581
 582                 sage: J = JordanSpinSimpleEJA(4)
 583                 sage: J.one().degree()
 584                 1
 585                 sage: e0,e1,e2,e3 = J.gens()
 586                 sage: (e0 - e1).degree()
 587                 2
 588
 589             In the spin factor algebra (of rank two), all elements that
 590             aren't multiples of the identity are regular::
 591
 592                 sage: set_random_seed()
 593                 sage: n = ZZ.random_element(1,10)
 594                 sage: J = JordanSpinSimpleEJA(n)
 595                 sage: x = J.random_element()
 596                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 597                 True
 598
 599             """
 600             return self.span_of_powers().dimension()
 601
 602
 603         def minimal_polynomial(self):
 604             """
 605             EXAMPLES::
 606
 607                 sage: set_random_seed()
 608                 sage: x = random_eja().random_element()
 609                 sage: x.degree() == x.minimal_polynomial().degree()
 610                 True
 611
 612             ::
 613
 614                 sage: set_random_seed()
 615                 sage: x = random_eja().random_element()
 616                 sage: x.degree() == x.minimal_polynomial().degree()
 617                 True
 618
 619             The minimal polynomial and the characteristic polynomial coincide
 620             and are known (see Alizadeh, Example 11.11) for all elements of
 621             the spin factor algebra that aren't scalar multiples of the
 622             identity::
 623
 624                 sage: set_random_seed()
 625                 sage: n = ZZ.random_element(2,10)
 626                 sage: J = JordanSpinSimpleEJA(n)
 627                 sage: y = J.random_element()
 628                 sage: while y == y.coefficient(0)*J.one():
 629                 ....:     y = J.random_element()
 630                 sage: y0 = y.vector()[0]
 631                 sage: y_bar = y.vector()[1:]
 632                 sage: actual = y.minimal_polynomial()
 633                 sage: x = SR.symbol('x', domain='real')
 634                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 635                 sage: bool(actual == expected)
 636                 True
 637
 638             """
 639             # The element we're going to call "minimal_polynomial()" on.
 640             # Either myself, interpreted as an element of a finite-
 641             # dimensional algebra, or an element of an associative
 642             # subalgebra.
 643             elt = None
 644
 645             if self.parent().is_associative():
 646                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 647             else:
 648                 V = self.span_of_powers()
 649                 assoc_subalg = self.subalgebra_generated_by()
 650                 # Mis-design warning: the basis used for span_of_powers()
 651                 # and subalgebra_generated_by() must be the same, and in
 652                 # the same order!
 653                 elt = assoc_subalg(V.coordinates(self.vector()))
 654
 655             # Recursive call, but should work since elt lives in an
 656             # associative algebra.
 657             return elt.minimal_polynomial()
 658
 659
 660         def natural_representation(self):
 661             """
 662             Return a more-natural representation of this element.
 663
 664             Every finite-dimensional Euclidean Jordan Algebra is a
 665             direct sum of five simple algebras, four of which comprise
 666             Hermitian matrices. This method returns the original
 667             "natural" representation of this element as a Hermitian
 668             matrix, if it has one. If not, you get the usual representation.
 669
 670             EXAMPLES::
 671
 672                 sage: J = ComplexHermitianSimpleEJA(3)
 673                 sage: J.one()
 674                 e0 + e5 + e8
 675                 sage: J.one().natural_representation()
 676                 [1 0 0 0 0 0]
 677                 [0 1 0 0 0 0]
 678                 [0 0 1 0 0 0]
 679                 [0 0 0 1 0 0]
 680                 [0 0 0 0 1 0]
 681                 [0 0 0 0 0 1]
 682
 683             ::
 684
 685                 sage: J = QuaternionHermitianSimpleEJA(3)
 686                 sage: J.one()
 687                 e0 + e9 + e14
 688                 sage: J.one().natural_representation()
 689                 [1 0 0 0 0 0 0 0 0 0 0 0]
 690                 [0 1 0 0 0 0 0 0 0 0 0 0]
 691                 [0 0 1 0 0 0 0 0 0 0 0 0]
 692                 [0 0 0 1 0 0 0 0 0 0 0 0]
 693                 [0 0 0 0 1 0 0 0 0 0 0 0]
 694                 [0 0 0 0 0 1 0 0 0 0 0 0]
 695                 [0 0 0 0 0 0 1 0 0 0 0 0]
 696                 [0 0 0 0 0 0 0 1 0 0 0 0]
 697                 [0 0 0 0 0 0 0 0 1 0 0 0]
 698                 [0 0 0 0 0 0 0 0 0 1 0 0]
 699                 [0 0 0 0 0 0 0 0 0 0 1 0]
 700                 [0 0 0 0 0 0 0 0 0 0 0 1]
 701
 702             """
 703             B = self.parent().natural_basis()
 704             W = B[0].matrix_space()
 705             return W.linear_combination(zip(self.vector(), B))
 706
 707
 708         def operator_matrix(self):
 709             """
 710             Return the matrix that represents left- (or right-)
 711             multiplication by this element in the parent algebra.
 712
 713             We have to override this because the superclass method
 714             returns a matrix that acts on row vectors (that is, on
 715             the right).
 716
 717             EXAMPLES:
 718
 719             Test the first polarization identity from my notes, Koecher Chapter
 720             III, or from Baes (2.3)::
 721
 722                 sage: set_random_seed()
 723                 sage: J = random_eja()
 724                 sage: x = J.random_element()
 725                 sage: y = J.random_element()
 726                 sage: Lx = x.operator_matrix()
 727                 sage: Ly = y.operator_matrix()
 728                 sage: Lxx = (x*x).operator_matrix()
 729                 sage: Lxy = (x*y).operator_matrix()
 730                 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
 731                 True
 732
 733             Test the second polarization identity from my notes or from
 734             Baes (2.4)::
 735
 736                 sage: set_random_seed()
 737                 sage: J = random_eja()
 738                 sage: x = J.random_element()
 739                 sage: y = J.random_element()
 740                 sage: z = J.random_element()
 741                 sage: Lx = x.operator_matrix()
 742                 sage: Ly = y.operator_matrix()
 743                 sage: Lz = z.operator_matrix()
 744                 sage: Lzy = (z*y).operator_matrix()
 745                 sage: Lxy = (x*y).operator_matrix()
 746                 sage: Lxz = (x*z).operator_matrix()
 747                 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
 748                 True
 749
 750             Test the third polarization identity from my notes or from
 751             Baes (2.5)::
 752
 753                 sage: set_random_seed()
 754                 sage: J = random_eja()
 755                 sage: u = J.random_element()
 756                 sage: y = J.random_element()
 757                 sage: z = J.random_element()
 758                 sage: Lu = u.operator_matrix()
 759                 sage: Ly = y.operator_matrix()
 760                 sage: Lz = z.operator_matrix()
 761                 sage: Lzy = (z*y).operator_matrix()
 762                 sage: Luy = (u*y).operator_matrix()
 763                 sage: Luz = (u*z).operator_matrix()
 764                 sage: Luyz = (u*(y*z)).operator_matrix()
 765                 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
 766                 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
 767                 sage: bool(lhs == rhs)
 768                 True
 769
 770             """
 771             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 772             return fda_elt.matrix().transpose()
 773
 774
 775         def quadratic_representation(self, other=None):
 776             """
 777             Return the quadratic representation of this element.
 778
 779             EXAMPLES:
 780
 781             The explicit form in the spin factor algebra is given by
 782             Alizadeh's Example 11.12::
 783
 784                 sage: set_random_seed()
 785                 sage: n = ZZ.random_element(1,10)
 786                 sage: J = JordanSpinSimpleEJA(n)
 787                 sage: x = J.random_element()
 788                 sage: x_vec = x.vector()
 789                 sage: x0 = x_vec[0]
 790                 sage: x_bar = x_vec[1:]
 791                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 792                 sage: B = 2*x0*x_bar.row()
 793                 sage: C = 2*x0*x_bar.column()
 794                 sage: D = identity_matrix(QQ, n-1)
 795                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 796                 sage: D = D + 2*x_bar.tensor_product(x_bar)
 797                 sage: Q = block_matrix(2,2,[A,B,C,D])
 798                 sage: Q == x.quadratic_representation()
 799                 True
 800
 801             Test all of the properties from Theorem 11.2 in Alizadeh::
 802
 803                 sage: set_random_seed()
 804                 sage: J = random_eja()
 805                 sage: x = J.random_element()
 806                 sage: y = J.random_element()
 807
 808             Property 1:
 809
 810                 sage: actual = x.quadratic_representation(y)
 811                 sage: expected = ( (x+y).quadratic_representation()
 812                 ....:              -x.quadratic_representation()
 813                 ....:              -y.quadratic_representation() ) / 2
 814                 sage: actual == expected
 815                 True
 816
 817             Property 2:
 818
 819                 sage: alpha = QQ.random_element()
 820                 sage: actual = (alpha*x).quadratic_representation()
 821                 sage: expected = (alpha^2)*x.quadratic_representation()
 822                 sage: actual == expected
 823                 True
 824
 825             Property 5:
 826
 827                 sage: Qy = y.quadratic_representation()
 828                 sage: actual = J(Qy*x.vector()).quadratic_representation()
 829                 sage: expected = Qy*x.quadratic_representation()*Qy
 830                 sage: actual == expected
 831                 True
 832
 833             Property 6:
 834
 835                 sage: k = ZZ.random_element(1,10)
 836                 sage: actual = (x^k).quadratic_representation()
 837                 sage: expected = (x.quadratic_representation())^k
 838                 sage: actual == expected
 839                 True
 840
 841             """
 842             if other is None:
 843                 other=self
 844             elif not other in self.parent():
 845                 raise TypeError("'other' must live in the same algebra")
 846
 847             L = self.operator_matrix()
 848             M = other.operator_matrix()
 849             return ( L*M + M*L - (self*other).operator_matrix() )
 850
 851
 852         def span_of_powers(self):
 853             """
 854             Return the vector space spanned by successive powers of
 855             this element.
 856             """
 857             # The dimension of the subalgebra can't be greater than
 858             # the big algebra, so just put everything into a list
 859             # and let span() get rid of the excess.
 860             V = self.vector().parent()
 861             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 862
 863
 864         def subalgebra_generated_by(self):
 865             """
 866             Return the associative subalgebra of the parent EJA generated
 867             by this element.
 868
 869             TESTS::
 870
 871                 sage: set_random_seed()
 872                 sage: x = random_eja().random_element()
 873                 sage: x.subalgebra_generated_by().is_associative()
 874                 True
 875
 876             Squaring in the subalgebra should be the same thing as
 877             squaring in the superalgebra::
 878
 879                 sage: set_random_seed()
 880                 sage: x = random_eja().random_element()
 881                 sage: u = x.subalgebra_generated_by().random_element()
 882                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
 883                 True
 884
 885             """
 886             # First get the subspace spanned by the powers of myself...
 887             V = self.span_of_powers()
 888             F = self.base_ring()
 889
 890             # Now figure out the entries of the right-multiplication
 891             # matrix for the successive basis elements b0, b1,... of
 892             # that subspace.
 893             mats = []
 894             for b_right in V.basis():
 895                 eja_b_right = self.parent()(b_right)
 896                 b_right_rows = []
 897                 # The first row of the right-multiplication matrix by
 898                 # b1 is what we get if we apply that matrix to b1. The
 899                 # second row of the right multiplication matrix by b1
 900                 # is what we get when we apply that matrix to b2...
 901                 #
 902                 # IMPORTANT: this assumes that all vectors are COLUMN
 903                 # vectors, unlike our superclass (which uses row vectors).
 904                 for b_left in V.basis():
 905                     eja_b_left = self.parent()(b_left)
 906                     # Multiply in the original EJA, but then get the
 907                     # coordinates from the subalgebra in terms of its
 908                     # basis.
 909                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 910                     b_right_rows.append(this_row)
 911                 b_right_matrix = matrix(F, b_right_rows)
 912                 mats.append(b_right_matrix)
 913
 914             # It's an algebra of polynomials in one element, and EJAs
 915             # are power-associative.
 916             #
 917             # TODO: choose generator names intelligently.
 918             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 919
 920
 921         def subalgebra_idempotent(self):
 922             """
 923             Find an idempotent in the associative subalgebra I generate
 924             using Proposition 2.3.5 in Baes.
 925
 926             TESTS::
 927
 928                 sage: set_random_seed()
 929                 sage: J = eja_rn(5)
 930                 sage: c = J.random_element().subalgebra_idempotent()
 931                 sage: c^2 == c
 932                 True
 933                 sage: J = JordanSpinSimpleEJA(5)
 934                 sage: c = J.random_element().subalgebra_idempotent()
 935                 sage: c^2 == c
 936                 True
 937
 938             """
 939             if self.is_nilpotent():
 940                 raise ValueError("this only works with non-nilpotent elements!")
 941
 942             V = self.span_of_powers()
 943             J = self.subalgebra_generated_by()
 944             # Mis-design warning: the basis used for span_of_powers()
 945             # and subalgebra_generated_by() must be the same, and in
 946             # the same order!
 947             u = J(V.coordinates(self.vector()))
 948
 949             # The image of the matrix of left-u^m-multiplication
 950             # will be minimal for some natural number s...
 951             s = 0
 952             minimal_dim = V.dimension()
 953             for i in xrange(1, V.dimension()):
 954                 this_dim = (u**i).operator_matrix().image().dimension()
 955                 if this_dim < minimal_dim:
 956                     minimal_dim = this_dim
 957                     s = i
 958
 959             # Now minimal_matrix should correspond to the smallest
 960             # non-zero subspace in Baes's (or really, Koecher's)
 961             # proposition.
 962             #
 963             # However, we need to restrict the matrix to work on the
 964             # subspace... or do we? Can't we just solve, knowing that
 965             # A(c) = u^(s+1) should have a solution in the big space,
 966             # too?
 967             #
 968             # Beware, solve_right() means that we're using COLUMN vectors.
 969             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 970             u_next = u**(s+1)
 971             A = u_next.operator_matrix()
 972             c_coordinates = A.solve_right(u_next.vector())
 973
 974             # Now c_coordinates is the idempotent we want, but it's in
 975             # the coordinate system of the subalgebra.
 976             #
 977             # We need the basis for J, but as elements of the parent algebra.
 978             #
 979             basis = [self.parent(v) for v in V.basis()]
 980             return self.parent().linear_combination(zip(c_coordinates, basis))
 981
 982
 983         def trace(self):
 984             """
 985             Return my trace, the sum of my eigenvalues.
 986
 987             EXAMPLES::
 988
 989                 sage: J = JordanSpinSimpleEJA(3)
 990                 sage: e0,e1,e2 = J.gens()
 991                 sage: x = e0 + e1 + e2
 992                 sage: x.trace()
 993                 2
 994
 995             """
 996             cs = self.characteristic_polynomial().coefficients(sparse=False)
 997             if len(cs) >= 2:
 998                 return -1*cs[-2]
 999             else:
1000                 raise ValueError('charpoly had fewer than 2 coefficients')
1001
1002
1003         def trace_inner_product(self, other):
1004             """
1005             Return the trace inner product of myself and ``other``.
1006             """
1007             if not other in self.parent():
1008                 raise TypeError("'other' must live in the same algebra")
1009
1010             return (self*other).trace()
1011
1012
1013 def eja_rn(dimension, field=QQ):
1014     """
1015     Return the Euclidean Jordan Algebra corresponding to the set
1016     `R^n` under the Hadamard product.
1017
1018     EXAMPLES:
1019
1020     This multiplication table can be verified by hand::
1021
1022         sage: J = eja_rn(3)
1023         sage: e0,e1,e2 = J.gens()
1024         sage: e0*e0
1025         e0
1026         sage: e0*e1
1027         0
1028         sage: e0*e2
1029         0
1030         sage: e1*e1
1031         e1
1032         sage: e1*e2
1033         0
1034         sage: e2*e2
1035         e2
1036
1037     """
1038     # The FiniteDimensionalAlgebra constructor takes a list of
1039     # matrices, the ith representing right multiplication by the ith
1040     # basis element in the vector space. So if e_1 = (1,0,0), then
1041     # right (Hadamard) multiplication of x by e_1 picks out the first
1042     # component of x; and likewise for the ith basis element e_i.
1043     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
1044            for i in xrange(dimension) ]
1045
1046     return FiniteDimensionalEuclideanJordanAlgebra(field,
1047                                                    Qs,
1048                                                    rank=dimension,
1049                                                    inner_product=_usual_ip)
1050
1051
1052
1053 def random_eja():
1054     """
1055     Return a "random" finite-dimensional Euclidean Jordan Algebra.
1056
1057     ALGORITHM:
1058
1059     For now, we choose a random natural number ``n`` (greater than zero)
1060     and then give you back one of the following:
1061
1062       * The cartesian product of the rational numbers ``n`` times; this is
1063         ``QQ^n`` with the Hadamard product.
1064
1065       * The Jordan spin algebra on ``QQ^n``.
1066
1067       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1068         product.
1069
1070       * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1071         in the space of ``2n``-by-``2n`` real symmetric matrices.
1072
1073       * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1074         in the space of ``4n``-by-``4n`` real symmetric matrices.
1075
1076     Later this might be extended to return Cartesian products of the
1077     EJAs above.
1078
1079     TESTS::
1080
1081         sage: random_eja()
1082         Euclidean Jordan algebra of degree...
1083
1084     """
1085     n = ZZ.random_element(1,5)
1086     constructor = choice([eja_rn,
1087                           JordanSpinSimpleEJA,
1088                           RealSymmetricSimpleEJA,
1089                           ComplexHermitianSimpleEJA,
1090                           QuaternionHermitianSimpleEJA])
1091     return constructor(n, field=QQ)
1092
1093
1094
1095 def _real_symmetric_basis(n, field=QQ):
1096     """
1097     Return a basis for the space of real symmetric n-by-n matrices.
1098     """
1099     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1100     # coordinates.
1101     S = []
1102     for i in xrange(n):
1103         for j in xrange(i+1):
1104             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1105             if i == j:
1106                 Sij = Eij
1107             else:
1108                 # Beware, orthogonal but not normalized!
1109                 Sij = Eij + Eij.transpose()
1110             S.append(Sij)
1111     return tuple(S)
1112
1113
1114 def _complex_hermitian_basis(n, field=QQ):
1115     """
1116     Returns a basis for the space of complex Hermitian n-by-n matrices.
1117
1118     TESTS::
1119
1120         sage: set_random_seed()
1121         sage: n = ZZ.random_element(1,5)
1122         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1123         True
1124
1125     """
1126     F = QuadraticField(-1, 'I')
1127     I = F.gen()
1128
1129     # This is like the symmetric case, but we need to be careful:
1130     #
1131     #   * We want conjugate-symmetry, not just symmetry.
1132     #   * The diagonal will (as a result) be real.
1133     #
1134     S = []
1135     for i in xrange(n):
1136         for j in xrange(i+1):
1137             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1138             if i == j:
1139                 Sij = _embed_complex_matrix(Eij)
1140                 S.append(Sij)
1141             else:
1142                 # Beware, orthogonal but not normalized! The second one
1143                 # has a minus because it's conjugated.
1144                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1145                 S.append(Sij_real)
1146                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1147                 S.append(Sij_imag)
1148     return tuple(S)
1149
1150
1151 def _quaternion_hermitian_basis(n, field=QQ):
1152     """
1153     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1154
1155     TESTS::
1156
1157         sage: set_random_seed()
1158         sage: n = ZZ.random_element(1,5)
1159         sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1160         True
1161
1162     """
1163     Q = QuaternionAlgebra(QQ,-1,-1)
1164     I,J,K = Q.gens()
1165
1166     # This is like the symmetric case, but we need to be careful:
1167     #
1168     #   * We want conjugate-symmetry, not just symmetry.
1169     #   * The diagonal will (as a result) be real.
1170     #
1171     S = []
1172     for i in xrange(n):
1173         for j in xrange(i+1):
1174             Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1175             if i == j:
1176                 Sij = _embed_quaternion_matrix(Eij)
1177                 S.append(Sij)
1178             else:
1179                 # Beware, orthogonal but not normalized! The second,
1180                 # third, and fourth ones have a minus because they're
1181                 # conjugated.
1182                 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1183                 S.append(Sij_real)
1184                 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1185                 S.append(Sij_I)
1186                 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1187                 S.append(Sij_J)
1188                 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1189                 S.append(Sij_K)
1190     return tuple(S)
1191
1192
1193 def _mat2vec(m):
1194         return vector(m.base_ring(), m.list())
1195
1196 def _vec2mat(v):
1197         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1198
1199 def _multiplication_table_from_matrix_basis(basis):
1200     """
1201     At least three of the five simple Euclidean Jordan algebras have the
1202     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1203     multiplication on the right is matrix multiplication. Given a basis
1204     for the underlying matrix space, this function returns a
1205     multiplication table (obtained by looping through the basis
1206     elements) for an algebra of those matrices. A reordered copy
1207     of the basis is also returned to work around the fact that
1208     the ``span()`` in this function will change the order of the basis
1209     from what we think it is, to... something else.
1210     """
1211     # In S^2, for example, we nominally have four coordinates even
1212     # though the space is of dimension three only. The vector space V
1213     # is supposed to hold the entire long vector, and the subspace W
1214     # of V will be spanned by the vectors that arise from symmetric
1215     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1216     field = basis[0].base_ring()
1217     dimension = basis[0].nrows()
1218
1219     V = VectorSpace(field, dimension**2)
1220     W = V.span( _mat2vec(s) for s in basis )
1221
1222     # Taking the span above reorders our basis (thanks, jerk!) so we
1223     # need to put our "matrix basis" in the same order as the
1224     # (reordered) vector basis.
1225     S = tuple( _vec2mat(b) for b in W.basis() )
1226
1227     Qs = []
1228     for s in S:
1229         # Brute force the multiplication-by-s matrix by looping
1230         # through all elements of the basis and doing the computation
1231         # to find out what the corresponding row should be. BEWARE:
1232         # these multiplication tables won't be symmetric! It therefore
1233         # becomes REALLY IMPORTANT that the underlying algebra
1234         # constructor uses ROW vectors and not COLUMN vectors. That's
1235         # why we're computing rows here and not columns.
1236         Q_rows = []
1237         for t in S:
1238             this_row = _mat2vec((s*t + t*s)/2)
1239             Q_rows.append(W.coordinates(this_row))
1240         Q = matrix(field, W.dimension(), Q_rows)
1241         Qs.append(Q)
1242
1243     return (Qs, S)
1244
1245
1246 def _embed_complex_matrix(M):
1247     """
1248     Embed the n-by-n complex matrix ``M`` into the space of real
1249     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1250     bi` to the block matrix ``[[a,b],[-b,a]]``.
1251
1252     EXAMPLES::
1253
1254         sage: F = QuadraticField(-1,'i')
1255         sage: x1 = F(4 - 2*i)
1256         sage: x2 = F(1 + 2*i)
1257         sage: x3 = F(-i)
1258         sage: x4 = F(6)
1259         sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1260         sage: _embed_complex_matrix(M)
1261         [ 4 -2| 1  2]
1262         [ 2  4|-2  1]
1263         [-----+-----]
1264         [ 0 -1| 6  0]
1265         [ 1  0| 0  6]
1266
1267     TESTS:
1268
1269     Embedding is a homomorphism (isomorphism, in fact)::
1270
1271         sage: set_random_seed()
1272         sage: n = ZZ.random_element(5)
1273         sage: F = QuadraticField(-1, 'i')
1274         sage: X = random_matrix(F, n)
1275         sage: Y = random_matrix(F, n)
1276         sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1277         sage: expected = _embed_complex_matrix(X*Y)
1278         sage: actual == expected
1279         True
1280
1281     """
1282     n = M.nrows()
1283     if M.ncols() != n:
1284         raise ValueError("the matrix 'M' must be square")
1285     field = M.base_ring()
1286     blocks = []
1287     for z in M.list():
1288         a = z.real()
1289         b = z.imag()
1290         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1291
1292     # We can drop the imaginaries here.
1293     return block_matrix(field.base_ring(), n, blocks)
1294
1295
1296 def _unembed_complex_matrix(M):
1297     """
1298     The inverse of _embed_complex_matrix().
1299
1300     EXAMPLES::
1301
1302         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1303         ....:                 [-2,  1,  -4,  3],
1304         ....:                 [ 9,  10, 11, 12],
1305         ....:                 [-10, 9, -12, 11] ])
1306         sage: _unembed_complex_matrix(A)
1307         [  2*i + 1   4*i + 3]
1308         [ 10*i + 9 12*i + 11]
1309
1310     TESTS:
1311
1312     Unembedding is the inverse of embedding::
1313
1314         sage: set_random_seed()
1315         sage: F = QuadraticField(-1, 'i')
1316         sage: M = random_matrix(F, 3)
1317         sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1318         True
1319
1320     """
1321     n = ZZ(M.nrows())
1322     if M.ncols() != n:
1323         raise ValueError("the matrix 'M' must be square")
1324     if not n.mod(2).is_zero():
1325         raise ValueError("the matrix 'M' must be a complex embedding")
1326
1327     F = QuadraticField(-1, 'i')
1328     i = F.gen()
1329
1330     # Go top-left to bottom-right (reading order), converting every
1331     # 2-by-2 block we see to a single complex element.
1332     elements = []
1333     for k in xrange(n/2):
1334         for j in xrange(n/2):
1335             submat = M[2*k:2*k+2,2*j:2*j+2]
1336             if submat[0,0] != submat[1,1]:
1337                 raise ValueError('bad on-diagonal submatrix')
1338             if submat[0,1] != -submat[1,0]:
1339                 raise ValueError('bad off-diagonal submatrix')
1340             z = submat[0,0] + submat[0,1]*i
1341             elements.append(z)
1342
1343     return matrix(F, n/2, elements)
1344
1345
1346 def _embed_quaternion_matrix(M):
1347     """
1348     Embed the n-by-n quaternion matrix ``M`` into the space of real
1349     matrices of size 4n-by-4n by first sending each quaternion entry
1350     `z = a + bi + cj + dk` to the block-complex matrix
1351     ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1352     a real matrix.
1353
1354     EXAMPLES::
1355
1356         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1357         sage: i,j,k = Q.gens()
1358         sage: x = 1 + 2*i + 3*j + 4*k
1359         sage: M = matrix(Q, 1, [[x]])
1360         sage: _embed_quaternion_matrix(M)
1361         [ 1  2  3  4]
1362         [-2  1 -4  3]
1363         [-3  4  1 -2]
1364         [-4 -3  2  1]
1365
1366     Embedding is a homomorphism (isomorphism, in fact)::
1367
1368         sage: set_random_seed()
1369         sage: n = ZZ.random_element(5)
1370         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1371         sage: X = random_matrix(Q, n)
1372         sage: Y = random_matrix(Q, n)
1373         sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1374         sage: expected = _embed_quaternion_matrix(X*Y)
1375         sage: actual == expected
1376         True
1377
1378     """
1379     quaternions = M.base_ring()
1380     n = M.nrows()
1381     if M.ncols() != n:
1382         raise ValueError("the matrix 'M' must be square")
1383
1384     F = QuadraticField(-1, 'i')
1385     i = F.gen()
1386
1387     blocks = []
1388     for z in M.list():
1389         t = z.coefficient_tuple()
1390         a = t[0]
1391         b = t[1]
1392         c = t[2]
1393         d = t[3]
1394         cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1395                                     [-c + d*i, a - b*i]])
1396         blocks.append(_embed_complex_matrix(cplx_matrix))
1397
1398     # We should have real entries by now, so use the realest field
1399     # we've got for the return value.
1400     return block_matrix(quaternions.base_ring(), n, blocks)
1401
1402
1403 def _unembed_quaternion_matrix(M):
1404     """
1405     The inverse of _embed_quaternion_matrix().
1406
1407     EXAMPLES::
1408
1409         sage: M = matrix(QQ, [[ 1,  2,  3,  4],
1410         ....:                 [-2,  1, -4,  3],
1411         ....:                 [-3,  4,  1, -2],
1412         ....:                 [-4, -3,  2,  1]])
1413         sage: _unembed_quaternion_matrix(M)
1414         [1 + 2*i + 3*j + 4*k]
1415
1416     TESTS:
1417
1418     Unembedding is the inverse of embedding::
1419
1420         sage: set_random_seed()
1421         sage: Q = QuaternionAlgebra(QQ, -1, -1)
1422         sage: M = random_matrix(Q, 3)
1423         sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1424         True
1425
1426     """
1427     n = ZZ(M.nrows())
1428     if M.ncols() != n:
1429         raise ValueError("the matrix 'M' must be square")
1430     if not n.mod(4).is_zero():
1431         raise ValueError("the matrix 'M' must be a complex embedding")
1432
1433     Q = QuaternionAlgebra(QQ,-1,-1)
1434     i,j,k = Q.gens()
1435
1436     # Go top-left to bottom-right (reading order), converting every
1437     # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1438     # quaternion block.
1439     elements = []
1440     for l in xrange(n/4):
1441         for m in xrange(n/4):
1442             submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1443             if submat[0,0] != submat[1,1].conjugate():
1444                 raise ValueError('bad on-diagonal submatrix')
1445             if submat[0,1] != -submat[1,0].conjugate():
1446                 raise ValueError('bad off-diagonal submatrix')
1447             z  = submat[0,0].real() + submat[0,0].imag()*i
1448             z += submat[0,1].real()*j + submat[0,1].imag()*k
1449             elements.append(z)
1450
1451     return matrix(Q, n/4, elements)
1452
1453
1454 # The usual inner product on R^n.
1455 def _usual_ip(x,y):
1456     return x.vector().inner_product(y.vector())
1457
1458 # The inner product used for the real symmetric simple EJA.
1459 # We keep it as a separate function because e.g. the complex
1460 # algebra uses the same inner product, except divided by 2.
1461 def _matrix_ip(X,Y):
1462     X_mat = X.natural_representation()
1463     Y_mat = Y.natural_representation()
1464     return (X_mat*Y_mat).trace()
1465
1466
1467 def RealSymmetricSimpleEJA(n, field=QQ):
1468     """
1469     The rank-n simple EJA consisting of real symmetric n-by-n
1470     matrices, the usual symmetric Jordan product, and the trace inner
1471     product. It has dimension `(n^2 + n)/2` over the reals.
1472
1473     EXAMPLES::
1474
1475         sage: J = RealSymmetricSimpleEJA(2)
1476         sage: e0, e1, e2 = J.gens()
1477         sage: e0*e0
1478         e0
1479         sage: e1*e1
1480         e0 + e2
1481         sage: e2*e2
1482         e2
1483
1484     TESTS:
1485
1486     The degree of this algebra is `(n^2 + n) / 2`::
1487
1488         sage: set_random_seed()
1489         sage: n = ZZ.random_element(1,5)
1490         sage: J = RealSymmetricSimpleEJA(n)
1491         sage: J.degree() == (n^2 + n)/2
1492         True
1493
1494     The Jordan multiplication is what we think it is::
1495
1496         sage: set_random_seed()
1497         sage: n = ZZ.random_element(1,5)
1498         sage: J = RealSymmetricSimpleEJA(n)
1499         sage: x = J.random_element()
1500         sage: y = J.random_element()
1501         sage: actual = (x*y).natural_representation()
1502         sage: X = x.natural_representation()
1503         sage: Y = y.natural_representation()
1504         sage: expected = (X*Y + Y*X)/2
1505         sage: actual == expected
1506         True
1507         sage: J(expected) == x*y
1508         True
1509
1510     """
1511     S = _real_symmetric_basis(n, field=field)
1512     (Qs, T) = _multiplication_table_from_matrix_basis(S)
1513
1514     return FiniteDimensionalEuclideanJordanAlgebra(field,
1515                                                    Qs,
1516                                                    rank=n,
1517                                                    natural_basis=T,
1518                                                    inner_product=_matrix_ip)
1519
1520
1521 def ComplexHermitianSimpleEJA(n, field=QQ):
1522     """
1523     The rank-n simple EJA consisting of complex Hermitian n-by-n
1524     matrices over the real numbers, the usual symmetric Jordan product,
1525     and the real-part-of-trace inner product. It has dimension `n^2` over
1526     the reals.
1527
1528     TESTS:
1529
1530     The degree of this algebra is `n^2`::
1531
1532         sage: set_random_seed()
1533         sage: n = ZZ.random_element(1,5)
1534         sage: J = ComplexHermitianSimpleEJA(n)
1535         sage: J.degree() == n^2
1536         True
1537
1538     The Jordan multiplication is what we think it is::
1539
1540         sage: set_random_seed()
1541         sage: n = ZZ.random_element(1,5)
1542         sage: J = ComplexHermitianSimpleEJA(n)
1543         sage: x = J.random_element()
1544         sage: y = J.random_element()
1545         sage: actual = (x*y).natural_representation()
1546         sage: X = x.natural_representation()
1547         sage: Y = y.natural_representation()
1548         sage: expected = (X*Y + Y*X)/2
1549         sage: actual == expected
1550         True
1551         sage: J(expected) == x*y
1552         True
1553
1554     """
1555     S = _complex_hermitian_basis(n)
1556     (Qs, T) = _multiplication_table_from_matrix_basis(S)
1557
1558     # Since a+bi on the diagonal is represented as
1559     #
1560     #   a + bi  = [  a  b  ]
1561     #             [ -b  a  ],
1562     #
1563     # we'll double-count the "a" entries if we take the trace of
1564     # the embedding.
1565     ip = lambda X,Y: _matrix_ip(X,Y)/2
1566
1567     return FiniteDimensionalEuclideanJordanAlgebra(field,
1568                                                    Qs,
1569                                                    rank=n,
1570                                                    natural_basis=T,
1571                                                    inner_product=ip)
1572
1573
1574 def QuaternionHermitianSimpleEJA(n, field=QQ):
1575     """
1576     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1577     matrices, the usual symmetric Jordan product, and the
1578     real-part-of-trace inner product. It has dimension `2n^2 - n` over
1579     the reals.
1580
1581     TESTS:
1582
1583     The degree of this algebra is `n^2`::
1584
1585         sage: set_random_seed()
1586         sage: n = ZZ.random_element(1,5)
1587         sage: J = QuaternionHermitianSimpleEJA(n)
1588         sage: J.degree() == 2*(n^2) - n
1589         True
1590
1591     The Jordan multiplication is what we think it is::
1592
1593         sage: set_random_seed()
1594         sage: n = ZZ.random_element(1,5)
1595         sage: J = QuaternionHermitianSimpleEJA(n)
1596         sage: x = J.random_element()
1597         sage: y = J.random_element()
1598         sage: actual = (x*y).natural_representation()
1599         sage: X = x.natural_representation()
1600         sage: Y = y.natural_representation()
1601         sage: expected = (X*Y + Y*X)/2
1602         sage: actual == expected
1603         True
1604         sage: J(expected) == x*y
1605         True
1606
1607     """
1608     S = _quaternion_hermitian_basis(n)
1609     (Qs, T) = _multiplication_table_from_matrix_basis(S)
1610
1611     # Since a+bi+cj+dk on the diagonal is represented as
1612     #
1613     #   a + bi +cj + dk = [  a  b  c  d]
1614     #                     [ -b  a -d  c]
1615     #                     [ -c  d  a -b]
1616     #                     [ -d -c  b  a],
1617     #
1618     # we'll quadruple-count the "a" entries if we take the trace of
1619     # the embedding.
1620     ip = lambda X,Y: _matrix_ip(X,Y)/4
1621
1622     return FiniteDimensionalEuclideanJordanAlgebra(field,
1623                                                    Qs,
1624                                                    rank=n,
1625                                                    natural_basis=T,
1626                                                    inner_product=ip)
1627
1628
1629 def OctonionHermitianSimpleEJA(n):
1630     """
1631     This shit be crazy. It has dimension 27 over the reals.
1632     """
1633     n = 3
1634     pass
1635
1636 def JordanSpinSimpleEJA(n, field=QQ):
1637     """
1638     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1639     with the usual inner product and jordan product ``x*y =
1640     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1641     the reals.
1642
1643     EXAMPLES:
1644
1645     This multiplication table can be verified by hand::
1646
1647         sage: J = JordanSpinSimpleEJA(4)
1648         sage: e0,e1,e2,e3 = J.gens()
1649         sage: e0*e0
1650         e0
1651         sage: e0*e1
1652         e1
1653         sage: e0*e2
1654         e2
1655         sage: e0*e3
1656         e3
1657         sage: e1*e2
1658         0
1659         sage: e1*e3
1660         0
1661         sage: e2*e3
1662         0
1663
1664     In one dimension, this is the reals under multiplication::
1665
1666         sage: J1 = JordanSpinSimpleEJA(1)
1667         sage: J2 = eja_rn(1)
1668         sage: J1 == J2
1669         True
1670
1671     """
1672     Qs = []
1673     id_matrix = identity_matrix(field, n)
1674     for i in xrange(n):
1675         ei = id_matrix.column(i)
1676         Qi = zero_matrix(field, n)
1677         Qi.set_row(0, ei)
1678         Qi.set_column(0, ei)
1679         Qi += diagonal_matrix(n, [ei[0]]*n)
1680         # The addition of the diagonal matrix adds an extra ei[0] in the
1681         # upper-left corner of the matrix.
1682         Qi[0,0] = Qi[0,0] * ~field(2)
1683         Qs.append(Qi)
1684
1685     # The rank of the spin factor algebra is two, UNLESS we're in a
1686     # one-dimensional ambient space (the rank is bounded by the
1687     # ambient dimension).
1688     return FiniteDimensionalEuclideanJordanAlgebra(field,
1689                                                    Qs,
1690                                                    rank=min(n,2),
1691                                                    inner_product=_usual_ip)