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