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 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_morphism import FiniteDimensionalAlgebraMorphism
  15
  16
  17 class FiniteDimensionalEuclideanJordanAlgebraMorphism(FiniteDimensionalAlgebraMorphism):
  18     """
  19     A very thin wrapper around FiniteDimensionalAlgebraMorphism that
  20     does only two things:
  21
  22       1. Avoids the ``unitary`` and ``check`` arguments to the constructor
  23          that will always be ``False``. This is necessary because these
  24          are homomorphisms with respect to ADDITION, but the SageMath
  25          machinery wants to check that they're homomorphisms with respect
  26          to (Jordan) MULTIPLICATION. That obviously doesn't work.
  27
  28       2. Inputs and outputs the underlying matrix with respect to COLUMN
  29          vectors, unlike the parent class.
  30
  31     If this seems a bit heavyweight, it is. I would have been happy to
  32     use a the ring morphism that underlies the finite-dimensional
  33     algebra morphism, but they don't seem to be callable on elements of
  34     our EJA.
  35     """
  36     def __init__(self, parent, f):
  37         FiniteDimensionalAlgebraMorphism.__init__(self,
  38                                                   parent,
  39                                                   f.transpose(),
  40                                                   unitary=False,
  41                                                   check=False)
  42
  43
  44     def _repr_(self):
  45         """
  46         We override only the representation that is shown to the user,
  47         because we want the matrix to be with respect to COLUMN vectors.
  48
  49         TESTS:
  50
  51         Ensure that we see the transpose of the underlying matrix object:
  52
  53             sage: J = RealSymmetricEJA(3)
  54             sage: x = J.linear_combination(zip(range(len(J.gens())), J.gens()))
  55             sage: L = x.operator()
  56             sage: L
  57             Morphism from Euclidean Jordan algebra of degree 6 over Rational
  58             Field to Euclidean Jordan algebra of degree 6 over Rational Field
  59             given by matrix
  60             [  0   1   2   0   0   0]
  61             [1/2 3/2   2 1/2   1   0]
  62             [  1   2 5/2   0 1/2   1]
  63             [  0   1   0   3   4   0]
  64             [  0   1 1/2   2   4   2]
  65             [  0   0   2   0   4   5]
  66             sage: L._matrix
  67             [  0 1/2   1   0   0   0]
  68             [  1 3/2   2   1   1   0]
  69             [  2   2 5/2   0 1/2   2]
  70             [  0 1/2   0   3   2   0]
  71             [  0   1 1/2   4   4   4]
  72             [  0   0   1   0   2   5]
  73
  74         """
  75         return "Morphism from {} to {} given by matrix\n{}".format(
  76             self.domain(), self.codomain(), self.matrix())
  77
  78     def matrix(self):
  79         """
  80         Return the matrix of this morphism with respect to a left-action
  81         on column vectors.
  82         """
  83         return FiniteDimensionalAlgebraMorphism.matrix(self).transpose()
  84
  85
  86 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
  87     @staticmethod
  88     def __classcall_private__(cls,
  89                               field,
  90                               mult_table,
  91                               names='e',
  92                               assume_associative=False,
  93                               category=None,
  94                               rank=None,
  95                               natural_basis=None):
  96         n = len(mult_table)
  97         mult_table = [b.base_extend(field) for b in mult_table]
  98         for b in mult_table:
  99             b.set_immutable()
 100             if not (is_Matrix(b) and b.dimensions() == (n, n)):
 101                 raise ValueError("input is not a multiplication table")
 102         mult_table = tuple(mult_table)
 103
 104         cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
 105         cat.or_subcategory(category)
 106         if assume_associative:
 107             cat = cat.Associative()
 108
 109         names = normalize_names(n, names)
 110
 111         fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls)
 112         return fda.__classcall__(cls,
 113                                  field,
 114                                  mult_table,
 115                                  assume_associative=assume_associative,
 116                                  names=names,
 117                                  category=cat,
 118                                  rank=rank,
 119                                  natural_basis=natural_basis)
 120
 121
 122     def __init__(self,
 123                  field,
 124                  mult_table,
 125                  names='e',
 126                  assume_associative=False,
 127                  category=None,
 128                  rank=None,
 129                  natural_basis=None):
 130         """
 131         EXAMPLES:
 132
 133         By definition, Jordan multiplication commutes::
 134
 135             sage: set_random_seed()
 136             sage: J = random_eja()
 137             sage: x = J.random_element()
 138             sage: y = J.random_element()
 139             sage: x*y == y*x
 140             True
 141
 142         """
 143         self._rank = rank
 144         self._natural_basis = natural_basis
 145         self._multiplication_table = mult_table
 146         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
 147         fda.__init__(field,
 148                      mult_table,
 149                      names=names,
 150                      category=category)
 151
 152
 153     def _repr_(self):
 154         """
 155         Return a string representation of ``self``.
 156         """
 157         fmt = "Euclidean Jordan algebra of degree {} over {}"
 158         return fmt.format(self.degree(), self.base_ring())
 159
 160
 161     def _a_regular_element(self):
 162         """
 163         Guess a regular element. Needed to compute the basis for our
 164         characteristic polynomial coefficients.
 165         """
 166         gs = self.gens()
 167         z = self.sum( (i+1)*gs[i] for i in range(len(gs)) )
 168         if not z.is_regular():
 169             raise ValueError("don't know a regular element")
 170         return z
 171
 172
 173     @cached_method
 174     def _charpoly_basis_space(self):
 175         """
 176         Return the vector space spanned by the basis used in our
 177         characteristic polynomial coefficients. This is used not only to
 178         compute those coefficients, but also any time we need to
 179         evaluate the coefficients (like when we compute the trace or
 180         determinant).
 181         """
 182         z = self._a_regular_element()
 183         V = z.vector().parent().ambient_vector_space()
 184         V1 = V.span_of_basis( (z**k).vector() for k in range(self.rank()) )
 185         b =  (V1.basis() + V1.complement().basis())
 186         return V.span_of_basis(b)
 187
 188
 189     @cached_method
 190     def _charpoly_coeff(self, i):
 191         """
 192         Return the coefficient polynomial "a_{i}" of this algebra's
 193         general characteristic polynomial.
 194
 195         Having this be a separate cached method lets us compute and
 196         store the trace/determinant (a_{r-1} and a_{0} respectively)
 197         separate from the entire characteristic polynomial.
 198         """
 199         (A_of_x, x, xr, detA) = self._charpoly_matrix_system()
 200         R = A_of_x.base_ring()
 201         if i >= self.rank():
 202             # Guaranteed by theory
 203             return R.zero()
 204
 205         # Danger: the in-place modification is done for performance
 206         # reasons (reconstructing a matrix with huge polynomial
 207         # entries is slow), but I don't know how cached_method works,
 208         # so it's highly possible that we're modifying some global
 209         # list variable by reference, here. In other words, you
 210         # probably shouldn't call this method twice on the same
 211         # algebra, at the same time, in two threads
 212         Ai_orig = A_of_x.column(i)
 213         A_of_x.set_column(i,xr)
 214         numerator = A_of_x.det()
 215         A_of_x.set_column(i,Ai_orig)
 216
 217         # We're relying on the theory here to ensure that each a_i is
 218         # indeed back in R, and the added negative signs are to make
 219         # the whole charpoly expression sum to zero.
 220         return R(-numerator/detA)
 221
 222
 223     @cached_method
 224     def _charpoly_matrix_system(self):
 225         """
 226         Compute the matrix whose entries A_ij are polynomials in
 227         X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
 228         corresponding to `x^r` and the determinent of the matrix A =
 229         [A_ij]. In other words, all of the fixed (cachable) data needed
 230         to compute the coefficients of the characteristic polynomial.
 231         """
 232         r = self.rank()
 233         n = self.dimension()
 234
 235         # Construct a new algebra over a multivariate polynomial ring...
 236         names = ['X' + str(i) for i in range(1,n+1)]
 237         R = PolynomialRing(self.base_ring(), names)
 238         J = FiniteDimensionalEuclideanJordanAlgebra(R,
 239                                                     self._multiplication_table,
 240                                                     rank=r)
 241
 242         idmat = identity_matrix(J.base_ring(), n)
 243
 244         W = self._charpoly_basis_space()
 245         W = W.change_ring(R.fraction_field())
 246
 247         # Starting with the standard coordinates x = (X1,X2,...,Xn)
 248         # and then converting the entries to W-coordinates allows us
 249         # to pass in the standard coordinates to the charpoly and get
 250         # back the right answer. Specifically, with x = (X1,X2,...,Xn),
 251         # we have
 252         #
 253         #   W.coordinates(x^2) eval'd at (standard z-coords)
 254         #     =
 255         #   W-coords of (z^2)
 256         #     =
 257         #   W-coords of (standard coords of x^2 eval'd at std-coords of z)
 258         #
 259         # We want the middle equivalent thing in our matrix, but use
 260         # the first equivalent thing instead so that we can pass in
 261         # standard coordinates.
 262         x = J(vector(R, R.gens()))
 263         l1 = [column_matrix(W.coordinates((x**k).vector())) for k in range(r)]
 264         l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)]
 265         A_of_x = block_matrix(R, 1, n, (l1 + l2))
 266         xr = W.coordinates((x**r).vector())
 267         return (A_of_x, x, xr, A_of_x.det())
 268
 269
 270     @cached_method
 271     def characteristic_polynomial(self):
 272         """
 273
 274         .. WARNING::
 275
 276             This implementation doesn't guarantee that the polynomial
 277             denominator in the coefficients is not identically zero, so
 278             theoretically it could crash. The way that this is handled
 279             in e.g. Faraut and Koranyi is to use a basis that guarantees
 280             the denominator is non-zero. But, doing so requires knowledge
 281             of at least one regular element, and we don't even know how
 282             to do that. The trade-off is that, if we use the standard basis,
 283             the resulting polynomial will accept the "usual" coordinates. In
 284             other words, we don't have to do a change of basis before e.g.
 285             computing the trace or determinant.
 286
 287         EXAMPLES:
 288
 289         The characteristic polynomial in the spin algebra is given in
 290         Alizadeh, Example 11.11::
 291
 292             sage: J = JordanSpinEJA(3)
 293             sage: p = J.characteristic_polynomial(); p
 294             X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
 295             sage: xvec = J.one().vector()
 296             sage: p(*xvec)
 297             t^2 - 2*t + 1
 298
 299         """
 300         r = self.rank()
 301         n = self.dimension()
 302
 303         # The list of coefficient polynomials a_1, a_2, ..., a_n.
 304         a = [ self._charpoly_coeff(i) for i in range(n) ]
 305
 306         # We go to a bit of trouble here to reorder the
 307         # indeterminates, so that it's easier to evaluate the
 308         # characteristic polynomial at x's coordinates and get back
 309         # something in terms of t, which is what we want.
 310         R = a[0].parent()
 311         S = PolynomialRing(self.base_ring(),'t')
 312         t = S.gen(0)
 313         S = PolynomialRing(S, R.variable_names())
 314         t = S(t)
 315
 316         # Note: all entries past the rth should be zero. The
 317         # coefficient of the highest power (x^r) is 1, but it doesn't
 318         # appear in the solution vector which contains coefficients
 319         # for the other powers (to make them sum to x^r).
 320         if (r < n):
 321             a[r] = 1 # corresponds to x^r
 322         else:
 323             # When the rank is equal to the dimension, trying to
 324             # assign a[r] goes out-of-bounds.
 325             a.append(1) # corresponds to x^r
 326
 327         return sum( a[k]*(t**k) for k in range(len(a)) )
 328
 329
 330     def inner_product(self, x, y):
 331         """
 332         The inner product associated with this Euclidean Jordan algebra.
 333
 334         Defaults to the trace inner product, but can be overridden by
 335         subclasses if they are sure that the necessary properties are
 336         satisfied.
 337
 338         EXAMPLES:
 339
 340         The inner product must satisfy its axiom for this algebra to truly
 341         be a Euclidean Jordan Algebra::
 342
 343             sage: set_random_seed()
 344             sage: J = random_eja()
 345             sage: x = J.random_element()
 346             sage: y = J.random_element()
 347             sage: z = J.random_element()
 348             sage: (x*y).inner_product(z) == y.inner_product(x*z)
 349             True
 350
 351         """
 352         if (not x in self) or (not y in self):
 353             raise TypeError("arguments must live in this algebra")
 354         return x.trace_inner_product(y)
 355
 356
 357     def natural_basis(self):
 358         """
 359         Return a more-natural representation of this algebra's basis.
 360
 361         Every finite-dimensional Euclidean Jordan Algebra is a direct
 362         sum of five simple algebras, four of which comprise Hermitian
 363         matrices. This method returns the original "natural" basis
 364         for our underlying vector space. (Typically, the natural basis
 365         is used to construct the multiplication table in the first place.)
 366
 367         Note that this will always return a matrix. The standard basis
 368         in `R^n` will be returned as `n`-by-`1` column matrices.
 369
 370         EXAMPLES::
 371
 372             sage: J = RealSymmetricEJA(2)
 373             sage: J.basis()
 374             Family (e0, e1, e2)
 375             sage: J.natural_basis()
 376             (
 377             [1 0]  [0 1]  [0 0]
 378             [0 0], [1 0], [0 1]
 379             )
 380
 381         ::
 382
 383             sage: J = JordanSpinEJA(2)
 384             sage: J.basis()
 385             Family (e0, e1)
 386             sage: J.natural_basis()
 387             (
 388             [1]  [0]
 389             [0], [1]
 390             )
 391
 392         """
 393         if self._natural_basis is None:
 394             return tuple( b.vector().column() for b in self.basis() )
 395         else:
 396             return self._natural_basis
 397
 398
 399     def rank(self):
 400         """
 401         Return the rank of this EJA.
 402         """
 403         if self._rank is None:
 404             raise ValueError("no rank specified at genesis")
 405         else:
 406             return self._rank
 407
 408
 409     class Element(FiniteDimensionalAlgebraElement):
 410         """
 411         An element of a Euclidean Jordan algebra.
 412         """
 413
 414         def __dir__(self):
 415             """
 416             Oh man, I should not be doing this. This hides the "disabled"
 417             methods ``left_matrix`` and ``matrix`` from introspection;
 418             in particular it removes them from tab-completion.
 419             """
 420             return filter(lambda s: s not in ['left_matrix', 'matrix'],
 421                           dir(self.__class__) )
 422
 423
 424         def __init__(self, A, elt=None):
 425             """
 426             EXAMPLES:
 427
 428             The identity in `S^n` is converted to the identity in the EJA::
 429
 430                 sage: J = RealSymmetricEJA(3)
 431                 sage: I = identity_matrix(QQ,3)
 432                 sage: J(I) == J.one()
 433                 True
 434
 435             This skew-symmetric matrix can't be represented in the EJA::
 436
 437                 sage: J = RealSymmetricEJA(3)
 438                 sage: A = matrix(QQ,3, lambda i,j: i-j)
 439                 sage: J(A)
 440                 Traceback (most recent call last):
 441                 ...
 442                 ArithmeticError: vector is not in free module
 443
 444             """
 445             # Goal: if we're given a matrix, and if it lives in our
 446             # parent algebra's "natural ambient space," convert it
 447             # into an algebra element.
 448             #
 449             # The catch is, we make a recursive call after converting
 450             # the given matrix into a vector that lives in the algebra.
 451             # This we need to try the parent class initializer first,
 452             # to avoid recursing forever if we're given something that
 453             # already fits into the algebra, but also happens to live
 454             # in the parent's "natural ambient space" (this happens with
 455             # vectors in R^n).
 456             try:
 457                 FiniteDimensionalAlgebraElement.__init__(self, A, elt)
 458             except ValueError:
 459                 natural_basis = A.natural_basis()
 460                 if elt in natural_basis[0].matrix_space():
 461                     # Thanks for nothing! Matrix spaces aren't vector
 462                     # spaces in Sage, so we have to figure out its
 463                     # natural-basis coordinates ourselves.
 464                     V = VectorSpace(elt.base_ring(), elt.nrows()**2)
 465                     W = V.span( _mat2vec(s) for s in natural_basis )
 466                     coords =  W.coordinates(_mat2vec(elt))
 467                     FiniteDimensionalAlgebraElement.__init__(self, A, coords)
 468
 469         def __pow__(self, n):
 470             """
 471             Return ``self`` raised to the power ``n``.
 472
 473             Jordan algebras are always power-associative; see for
 474             example Faraut and Koranyi, Proposition II.1.2 (ii).
 475
 476             .. WARNING:
 477
 478                 We have to override this because our superclass uses row vectors
 479                 instead of column vectors! We, on the other hand, assume column
 480                 vectors everywhere.
 481
 482             EXAMPLES::
 483
 484                 sage: set_random_seed()
 485                 sage: x = random_eja().random_element()
 486                 sage: x.operator_matrix()*x.vector() == (x^2).vector()
 487                 True
 488
 489             A few examples of power-associativity::
 490
 491                 sage: set_random_seed()
 492                 sage: x = random_eja().random_element()
 493                 sage: x*(x*x)*(x*x) == x^5
 494                 True
 495                 sage: (x*x)*(x*x*x) == x^5
 496                 True
 497
 498             We also know that powers operator-commute (Koecher, Chapter
 499             III, Corollary 1)::
 500
 501                 sage: set_random_seed()
 502                 sage: x = random_eja().random_element()
 503                 sage: m = ZZ.random_element(0,10)
 504                 sage: n = ZZ.random_element(0,10)
 505                 sage: Lxm = (x^m).operator_matrix()
 506                 sage: Lxn = (x^n).operator_matrix()
 507                 sage: Lxm*Lxn == Lxn*Lxm
 508                 True
 509
 510             """
 511             A = self.parent()
 512             if n == 0:
 513                 return A.one()
 514             elif n == 1:
 515                 return self
 516             else:
 517                 return A( (self.operator_matrix()**(n-1))*self.vector() )
 518
 519
 520         def apply_univariate_polynomial(self, p):
 521             """
 522             Apply the univariate polynomial ``p`` to this element.
 523
 524             A priori, SageMath won't allow us to apply a univariate
 525             polynomial to an element of an EJA, because we don't know
 526             that EJAs are rings (they are usually not associative). Of
 527             course, we know that EJAs are power-associative, so the
 528             operation is ultimately kosher. This function sidesteps
 529             the CAS to get the answer we want and expect.
 530
 531             EXAMPLES::
 532
 533                 sage: R = PolynomialRing(QQ, 't')
 534                 sage: t = R.gen(0)
 535                 sage: p = t^4 - t^3 + 5*t - 2
 536                 sage: J = RealCartesianProductEJA(5)
 537                 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
 538                 True
 539
 540             TESTS:
 541
 542             We should always get back an element of the algebra::
 543
 544                 sage: set_random_seed()
 545                 sage: p = PolynomialRing(QQ, 't').random_element()
 546                 sage: J = random_eja()
 547                 sage: x = J.random_element()
 548                 sage: x.apply_univariate_polynomial(p) in J
 549                 True
 550
 551             """
 552             if len(p.variables()) > 1:
 553                 raise ValueError("not a univariate polynomial")
 554             P = self.parent()
 555             R = P.base_ring()
 556             # Convert the coeficcients to the parent's base ring,
 557             # because a priori they might live in an (unnecessarily)
 558             # larger ring for which P.sum() would fail below.
 559             cs = [ R(c) for c in p.coefficients(sparse=False) ]
 560             return P.sum( cs[k]*(self**k) for k in range(len(cs)) )
 561
 562
 563         def characteristic_polynomial(self):
 564             """
 565             Return the characteristic polynomial of this element.
 566
 567             EXAMPLES:
 568
 569             The rank of `R^3` is three, and the minimal polynomial of
 570             the identity element is `(t-1)` from which it follows that
 571             the characteristic polynomial should be `(t-1)^3`::
 572
 573                 sage: J = RealCartesianProductEJA(3)
 574                 sage: J.one().characteristic_polynomial()
 575                 t^3 - 3*t^2 + 3*t - 1
 576
 577             Likewise, the characteristic of the zero element in the
 578             rank-three algebra `R^{n}` should be `t^{3}`::
 579
 580                 sage: J = RealCartesianProductEJA(3)
 581                 sage: J.zero().characteristic_polynomial()
 582                 t^3
 583
 584             The characteristic polynomial of an element should evaluate
 585             to zero on that element::
 586
 587                 sage: set_random_seed()
 588                 sage: x = RealCartesianProductEJA(3).random_element()
 589                 sage: p = x.characteristic_polynomial()
 590                 sage: x.apply_univariate_polynomial(p)
 591                 0
 592
 593             """
 594             p = self.parent().characteristic_polynomial()
 595             return p(*self.vector())
 596
 597
 598         def inner_product(self, other):
 599             """
 600             Return the parent algebra's inner product of myself and ``other``.
 601
 602             EXAMPLES:
 603
 604             The inner product in the Jordan spin algebra is the usual
 605             inner product on `R^n` (this example only works because the
 606             basis for the Jordan algebra is the standard basis in `R^n`)::
 607
 608                 sage: J = JordanSpinEJA(3)
 609                 sage: x = vector(QQ,[1,2,3])
 610                 sage: y = vector(QQ,[4,5,6])
 611                 sage: x.inner_product(y)
 612                 32
 613                 sage: J(x).inner_product(J(y))
 614                 32
 615
 616             The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
 617             multiplication is the usual matrix multiplication in `S^n`,
 618             so the inner product of the identity matrix with itself
 619             should be the `n`::
 620
 621                 sage: J = RealSymmetricEJA(3)
 622                 sage: J.one().inner_product(J.one())
 623                 3
 624
 625             Likewise, the inner product on `C^n` is `<X,Y> =
 626             Re(trace(X*Y))`, where we must necessarily take the real
 627             part because the product of Hermitian matrices may not be
 628             Hermitian::
 629
 630                 sage: J = ComplexHermitianEJA(3)
 631                 sage: J.one().inner_product(J.one())
 632                 3
 633
 634             Ditto for the quaternions::
 635
 636                 sage: J = QuaternionHermitianEJA(3)
 637                 sage: J.one().inner_product(J.one())
 638                 3
 639
 640             TESTS:
 641
 642             Ensure that we can always compute an inner product, and that
 643             it gives us back a real number::
 644
 645                 sage: set_random_seed()
 646                 sage: J = random_eja()
 647                 sage: x = J.random_element()
 648                 sage: y = J.random_element()
 649                 sage: x.inner_product(y) in RR
 650                 True
 651
 652             """
 653             P = self.parent()
 654             if not other in P:
 655                 raise TypeError("'other' must live in the same algebra")
 656
 657             return P.inner_product(self, other)
 658
 659
 660         def operator_commutes_with(self, other):
 661             """
 662             Return whether or not this element operator-commutes
 663             with ``other``.
 664
 665             EXAMPLES:
 666
 667             The definition of a Jordan algebra says that any element
 668             operator-commutes with its square::
 669
 670                 sage: set_random_seed()
 671                 sage: x = random_eja().random_element()
 672                 sage: x.operator_commutes_with(x^2)
 673                 True
 674
 675             TESTS:
 676
 677             Test Lemma 1 from Chapter III of Koecher::
 678
 679                 sage: set_random_seed()
 680                 sage: J = random_eja()
 681                 sage: u = J.random_element()
 682                 sage: v = J.random_element()
 683                 sage: lhs = u.operator_commutes_with(u*v)
 684                 sage: rhs = v.operator_commutes_with(u^2)
 685                 sage: lhs == rhs
 686                 True
 687
 688             """
 689             if not other in self.parent():
 690                 raise TypeError("'other' must live in the same algebra")
 691
 692             A = self.operator_matrix()
 693             B = other.operator_matrix()
 694             return (A*B == B*A)
 695
 696
 697         def det(self):
 698             """
 699             Return my determinant, the product of my eigenvalues.
 700
 701             EXAMPLES::
 702
 703                 sage: J = JordanSpinEJA(2)
 704                 sage: e0,e1 = J.gens()
 705                 sage: x = sum( J.gens() )
 706                 sage: x.det()
 707                 0
 708
 709             ::
 710
 711                 sage: J = JordanSpinEJA(3)
 712                 sage: e0,e1,e2 = J.gens()
 713                 sage: x = sum( J.gens() )
 714                 sage: x.det()
 715                 -1
 716
 717             TESTS:
 718
 719             An element is invertible if and only if its determinant is
 720             non-zero::
 721
 722                 sage: set_random_seed()
 723                 sage: x = random_eja().random_element()
 724                 sage: x.is_invertible() == (x.det() != 0)
 725                 True
 726
 727             """
 728             P = self.parent()
 729             r = P.rank()
 730             p = P._charpoly_coeff(0)
 731             # The _charpoly_coeff function already adds the factor of
 732             # -1 to ensure that _charpoly_coeff(0) is really what
 733             # appears in front of t^{0} in the charpoly. However,
 734             # we want (-1)^r times THAT for the determinant.
 735             return ((-1)**r)*p(*self.vector())
 736
 737
 738         def inverse(self):
 739             """
 740             Return the Jordan-multiplicative inverse of this element.
 741
 742             ALGORITHM:
 743
 744             We appeal to the quadratic representation as in Koecher's
 745             Theorem 12 in Chapter III, Section 5.
 746
 747             EXAMPLES:
 748
 749             The inverse in the spin factor algebra is given in Alizadeh's
 750             Example 11.11::
 751
 752                 sage: set_random_seed()
 753                 sage: n = ZZ.random_element(1,10)
 754                 sage: J = JordanSpinEJA(n)
 755                 sage: x = J.random_element()
 756                 sage: while x.is_zero():
 757                 ....:     x = J.random_element()
 758                 sage: x_vec = x.vector()
 759                 sage: x0 = x_vec[0]
 760                 sage: x_bar = x_vec[1:]
 761                 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
 762                 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
 763                 sage: x_inverse = coeff*inv_vec
 764                 sage: x.inverse() == J(x_inverse)
 765                 True
 766
 767             TESTS:
 768
 769             The identity element is its own inverse::
 770
 771                 sage: set_random_seed()
 772                 sage: J = random_eja()
 773                 sage: J.one().inverse() == J.one()
 774                 True
 775
 776             If an element has an inverse, it acts like one::
 777
 778                 sage: set_random_seed()
 779                 sage: J = random_eja()
 780                 sage: x = J.random_element()
 781                 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
 782                 True
 783
 784             The inverse of the inverse is what we started with::
 785
 786                 sage: set_random_seed()
 787                 sage: J = random_eja()
 788                 sage: x = J.random_element()
 789                 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
 790                 True
 791
 792             The zero element is never invertible::
 793
 794                 sage: set_random_seed()
 795                 sage: J = random_eja().zero().inverse()
 796                 Traceback (most recent call last):
 797                 ...
 798                 ValueError: element is not invertible
 799
 800             """
 801             if not self.is_invertible():
 802                 raise ValueError("element is not invertible")
 803
 804             P = self.parent()
 805             return P(self.quadratic_representation().inverse()*self.vector())
 806
 807
 808         def is_invertible(self):
 809             """
 810             Return whether or not this element is invertible.
 811
 812             We can't use the superclass method because it relies on
 813             the algebra being associative.
 814
 815             ALGORITHM:
 816
 817             The usual way to do this is to check if the determinant is
 818             zero, but we need the characteristic polynomial for the
 819             determinant. The minimal polynomial is a lot easier to get,
 820             so we use Corollary 2 in Chapter V of Koecher to check
 821             whether or not the paren't algebra's zero element is a root
 822             of this element's minimal polynomial.
 823
 824             TESTS:
 825
 826             The identity element is always invertible::
 827
 828                 sage: set_random_seed()
 829                 sage: J = random_eja()
 830                 sage: J.one().is_invertible()
 831                 True
 832
 833             The zero element is never invertible::
 834
 835                 sage: set_random_seed()
 836                 sage: J = random_eja()
 837                 sage: J.zero().is_invertible()
 838                 False
 839
 840             """
 841             zero = self.parent().zero()
 842             p = self.minimal_polynomial()
 843             return not (p(zero) == zero)
 844
 845
 846         def is_nilpotent(self):
 847             """
 848             Return whether or not some power of this element is zero.
 849
 850             The superclass method won't work unless we're in an
 851             associative algebra, and we aren't. However, we generate
 852             an assocoative subalgebra and we're nilpotent there if and
 853             only if we're nilpotent here (probably).
 854
 855             TESTS:
 856
 857             The identity element is never nilpotent::
 858
 859                 sage: set_random_seed()
 860                 sage: random_eja().one().is_nilpotent()
 861                 False
 862
 863             The additive identity is always nilpotent::
 864
 865                 sage: set_random_seed()
 866                 sage: random_eja().zero().is_nilpotent()
 867                 True
 868
 869             """
 870             # The element we're going to call "is_nilpotent()" on.
 871             # Either myself, interpreted as an element of a finite-
 872             # dimensional algebra, or an element of an associative
 873             # subalgebra.
 874             elt = None
 875
 876             if self.parent().is_associative():
 877                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 878             else:
 879                 V = self.span_of_powers()
 880                 assoc_subalg = self.subalgebra_generated_by()
 881                 # Mis-design warning: the basis used for span_of_powers()
 882                 # and subalgebra_generated_by() must be the same, and in
 883                 # the same order!
 884                 elt = assoc_subalg(V.coordinates(self.vector()))
 885
 886             # Recursive call, but should work since elt lives in an
 887             # associative algebra.
 888             return elt.is_nilpotent()
 889
 890
 891         def is_regular(self):
 892             """
 893             Return whether or not this is a regular element.
 894
 895             EXAMPLES:
 896
 897             The identity element always has degree one, but any element
 898             linearly-independent from it is regular::
 899
 900                 sage: J = JordanSpinEJA(5)
 901                 sage: J.one().is_regular()
 902                 False
 903                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 904                 sage: for x in J.gens():
 905                 ....:     (J.one() + x).is_regular()
 906                 False
 907                 True
 908                 True
 909                 True
 910                 True
 911
 912             """
 913             return self.degree() == self.parent().rank()
 914
 915
 916         def degree(self):
 917             """
 918             Compute the degree of this element the straightforward way
 919             according to the definition; by appending powers to a list
 920             and figuring out its dimension (that is, whether or not
 921             they're linearly dependent).
 922
 923             EXAMPLES::
 924
 925                 sage: J = JordanSpinEJA(4)
 926                 sage: J.one().degree()
 927                 1
 928                 sage: e0,e1,e2,e3 = J.gens()
 929                 sage: (e0 - e1).degree()
 930                 2
 931
 932             In the spin factor algebra (of rank two), all elements that
 933             aren't multiples of the identity are regular::
 934
 935                 sage: set_random_seed()
 936                 sage: n = ZZ.random_element(1,10)
 937                 sage: J = JordanSpinEJA(n)
 938                 sage: x = J.random_element()
 939                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 940                 True
 941
 942             """
 943             return self.span_of_powers().dimension()
 944
 945
 946         def left_matrix(self):
 947             """
 948             Our parent class defines ``left_matrix`` and ``matrix``
 949             methods whose names are misleading. We don't want them.
 950             """
 951             raise NotImplementedError("use operator_matrix() instead")
 952
 953         matrix = left_matrix
 954
 955
 956         def minimal_polynomial(self):
 957             """
 958             Return the minimal polynomial of this element,
 959             as a function of the variable `t`.
 960
 961             ALGORITHM:
 962
 963             We restrict ourselves to the associative subalgebra
 964             generated by this element, and then return the minimal
 965             polynomial of this element's operator matrix (in that
 966             subalgebra). This works by Baes Proposition 2.3.16.
 967
 968             TESTS:
 969
 970             The minimal polynomial of the identity and zero elements are
 971             always the same::
 972
 973                 sage: set_random_seed()
 974                 sage: J = random_eja()
 975                 sage: J.one().minimal_polynomial()
 976                 t - 1
 977                 sage: J.zero().minimal_polynomial()
 978                 t
 979
 980             The degree of an element is (by one definition) the degree
 981             of its minimal polynomial::
 982
 983                 sage: set_random_seed()
 984                 sage: x = random_eja().random_element()
 985                 sage: x.degree() == x.minimal_polynomial().degree()
 986                 True
 987
 988             The minimal polynomial and the characteristic polynomial coincide
 989             and are known (see Alizadeh, Example 11.11) for all elements of
 990             the spin factor algebra that aren't scalar multiples of the
 991             identity::
 992
 993                 sage: set_random_seed()
 994                 sage: n = ZZ.random_element(2,10)
 995                 sage: J = JordanSpinEJA(n)
 996                 sage: y = J.random_element()
 997                 sage: while y == y.coefficient(0)*J.one():
 998                 ....:     y = J.random_element()
 999                 sage: y0 = y.vector()[0]
1000                 sage: y_bar = y.vector()[1:]
1001                 sage: actual = y.minimal_polynomial()
1002                 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
1003                 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
1004                 sage: bool(actual == expected)
1005                 True
1006
1007             The minimal polynomial should always kill its element::
1008
1009                 sage: set_random_seed()
1010                 sage: x = random_eja().random_element()
1011                 sage: p = x.minimal_polynomial()
1012                 sage: x.apply_univariate_polynomial(p)
1013                 0
1014
1015             """
1016             V = self.span_of_powers()
1017             assoc_subalg = self.subalgebra_generated_by()
1018             # Mis-design warning: the basis used for span_of_powers()
1019             # and subalgebra_generated_by() must be the same, and in
1020             # the same order!
1021             elt = assoc_subalg(V.coordinates(self.vector()))
1022
1023             # We get back a symbolic polynomial in 'x' but want a real
1024             # polynomial in 't'.
1025             p_of_x = elt.operator_matrix().minimal_polynomial()
1026             return p_of_x.change_variable_name('t')
1027
1028
1029         def natural_representation(self):
1030             """
1031             Return a more-natural representation of this element.
1032
1033             Every finite-dimensional Euclidean Jordan Algebra is a
1034             direct sum of five simple algebras, four of which comprise
1035             Hermitian matrices. This method returns the original
1036             "natural" representation of this element as a Hermitian
1037             matrix, if it has one. If not, you get the usual representation.
1038
1039             EXAMPLES::
1040
1041                 sage: J = ComplexHermitianEJA(3)
1042                 sage: J.one()
1043                 e0 + e5 + e8
1044                 sage: J.one().natural_representation()
1045                 [1 0 0 0 0 0]
1046                 [0 1 0 0 0 0]
1047                 [0 0 1 0 0 0]
1048                 [0 0 0 1 0 0]
1049                 [0 0 0 0 1 0]
1050                 [0 0 0 0 0 1]
1051
1052             ::
1053
1054                 sage: J = QuaternionHermitianEJA(3)
1055                 sage: J.one()
1056                 e0 + e9 + e14
1057                 sage: J.one().natural_representation()
1058                 [1 0 0 0 0 0 0 0 0 0 0 0]
1059                 [0 1 0 0 0 0 0 0 0 0 0 0]
1060                 [0 0 1 0 0 0 0 0 0 0 0 0]
1061                 [0 0 0 1 0 0 0 0 0 0 0 0]
1062                 [0 0 0 0 1 0 0 0 0 0 0 0]
1063                 [0 0 0 0 0 1 0 0 0 0 0 0]
1064                 [0 0 0 0 0 0 1 0 0 0 0 0]
1065                 [0 0 0 0 0 0 0 1 0 0 0 0]
1066                 [0 0 0 0 0 0 0 0 1 0 0 0]
1067                 [0 0 0 0 0 0 0 0 0 1 0 0]
1068                 [0 0 0 0 0 0 0 0 0 0 1 0]
1069                 [0 0 0 0 0 0 0 0 0 0 0 1]
1070
1071             """
1072             B = self.parent().natural_basis()
1073             W = B[0].matrix_space()
1074             return W.linear_combination(zip(self.vector(), B))
1075
1076
1077         def operator_matrix(self):
1078             """
1079             Return the matrix that represents left- (or right-)
1080             multiplication by this element in the parent algebra.
1081
1082             We implement this ourselves to work around the fact that
1083             our parent class represents everything with row vectors.
1084
1085             EXAMPLES:
1086
1087             Test the first polarization identity from my notes, Koecher Chapter
1088             III, or from Baes (2.3)::
1089
1090                 sage: set_random_seed()
1091                 sage: J = random_eja()
1092                 sage: x = J.random_element()
1093                 sage: y = J.random_element()
1094                 sage: Lx = x.operator_matrix()
1095                 sage: Ly = y.operator_matrix()
1096                 sage: Lxx = (x*x).operator_matrix()
1097                 sage: Lxy = (x*y).operator_matrix()
1098                 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
1099                 True
1100
1101             Test the second polarization identity from my notes or from
1102             Baes (2.4)::
1103
1104                 sage: set_random_seed()
1105                 sage: J = random_eja()
1106                 sage: x = J.random_element()
1107                 sage: y = J.random_element()
1108                 sage: z = J.random_element()
1109                 sage: Lx = x.operator_matrix()
1110                 sage: Ly = y.operator_matrix()
1111                 sage: Lz = z.operator_matrix()
1112                 sage: Lzy = (z*y).operator_matrix()
1113                 sage: Lxy = (x*y).operator_matrix()
1114                 sage: Lxz = (x*z).operator_matrix()
1115                 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
1116                 True
1117
1118             Test the third polarization identity from my notes or from
1119             Baes (2.5)::
1120
1121                 sage: set_random_seed()
1122                 sage: J = random_eja()
1123                 sage: u = J.random_element()
1124                 sage: y = J.random_element()
1125                 sage: z = J.random_element()
1126                 sage: Lu = u.operator_matrix()
1127                 sage: Ly = y.operator_matrix()
1128                 sage: Lz = z.operator_matrix()
1129                 sage: Lzy = (z*y).operator_matrix()
1130                 sage: Luy = (u*y).operator_matrix()
1131                 sage: Luz = (u*z).operator_matrix()
1132                 sage: Luyz = (u*(y*z)).operator_matrix()
1133                 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
1134                 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
1135                 sage: bool(lhs == rhs)
1136                 True
1137
1138             """
1139             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
1140             return fda_elt.matrix().transpose()
1141
1142
1143         def quadratic_representation(self, other=None):
1144             """
1145             Return the quadratic representation of this element.
1146
1147             EXAMPLES:
1148
1149             The explicit form in the spin factor algebra is given by
1150             Alizadeh's Example 11.12::
1151
1152                 sage: set_random_seed()
1153                 sage: n = ZZ.random_element(1,10)
1154                 sage: J = JordanSpinEJA(n)
1155                 sage: x = J.random_element()
1156                 sage: x_vec = x.vector()
1157                 sage: x0 = x_vec[0]
1158                 sage: x_bar = x_vec[1:]
1159                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1160                 sage: B = 2*x0*x_bar.row()
1161                 sage: C = 2*x0*x_bar.column()
1162                 sage: D = identity_matrix(QQ, n-1)
1163                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1164                 sage: D = D + 2*x_bar.tensor_product(x_bar)
1165                 sage: Q = block_matrix(2,2,[A,B,C,D])
1166                 sage: Q == x.quadratic_representation()
1167                 True
1168
1169             Test all of the properties from Theorem 11.2 in Alizadeh::
1170
1171                 sage: set_random_seed()
1172                 sage: J = random_eja()
1173                 sage: x = J.random_element()
1174                 sage: y = J.random_element()
1175                 sage: Lx = x.operator_matrix()
1176                 sage: Lxx = (x*x).operator_matrix()
1177                 sage: Qx = x.quadratic_representation()
1178                 sage: Qy = y.quadratic_representation()
1179                 sage: Qxy = x.quadratic_representation(y)
1180                 sage: Qex = J.one().quadratic_representation(x)
1181                 sage: n = ZZ.random_element(10)
1182                 sage: Qxn = (x^n).quadratic_representation()
1183
1184             Property 1:
1185
1186                 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1187                 True
1188
1189             Property 2:
1190
1191                 sage: alpha = QQ.random_element()
1192                 sage: (alpha*x).quadratic_representation() == (alpha^2)*Qx
1193                 True
1194
1195             Property 3:
1196
1197                 sage: not x.is_invertible() or (
1198                 ....:     Qx*x.inverse().vector() == x.vector() )
1199                 True
1200
1201                 sage: not x.is_invertible() or (
1202                 ....:   Qx.inverse()
1203                 ....:   ==
1204                 ....:   x.inverse().quadratic_representation() )
1205                 True
1206
1207                 sage: Qxy*(J.one().vector()) == (x*y).vector()
1208                 True
1209
1210             Property 4:
1211
1212                 sage: not x.is_invertible() or (
1213                 ....:   x.quadratic_representation(x.inverse())*Qx
1214                 ....:   == Qx*x.quadratic_representation(x.inverse()) )
1215                 True
1216
1217                 sage: not x.is_invertible() or (
1218                 ....:   x.quadratic_representation(x.inverse())*Qx
1219                 ....:   ==
1220                 ....:   2*x.operator_matrix()*Qex - Qx )
1221                 True
1222
1223                 sage: 2*x.operator_matrix()*Qex - Qx == Lxx
1224                 True
1225
1226             Property 5:
1227
1228                 sage: J(Qy*x.vector()).quadratic_representation() == Qy*Qx*Qy
1229                 True
1230
1231             Property 6:
1232
1233                 sage: Qxn == (Qx)^n
1234                 True
1235
1236             Property 7:
1237
1238                 sage: not x.is_invertible() or (
1239                 ....:   Qx*x.inverse().operator_matrix() == Lx )
1240                 True
1241
1242             Property 8:
1243
1244                 sage: not x.operator_commutes_with(y) or (
1245                 ....:   J(Qx*y.vector())^n == J(Qxn*(y^n).vector()) )
1246                 True
1247
1248             """
1249             if other is None:
1250                 other=self
1251             elif not other in self.parent():
1252                 raise TypeError("'other' must live in the same algebra")
1253
1254             L = self.operator_matrix()
1255             M = other.operator_matrix()
1256             return ( L*M + M*L - (self*other).operator_matrix() )
1257
1258
1259         def span_of_powers(self):
1260             """
1261             Return the vector space spanned by successive powers of
1262             this element.
1263             """
1264             # The dimension of the subalgebra can't be greater than
1265             # the big algebra, so just put everything into a list
1266             # and let span() get rid of the excess.
1267             #
1268             # We do the extra ambient_vector_space() in case we're messing
1269             # with polynomials and the direct parent is a module.
1270             V = self.vector().parent().ambient_vector_space()
1271             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
1272
1273
1274         def subalgebra_generated_by(self):
1275             """
1276             Return the associative subalgebra of the parent EJA generated
1277             by this element.
1278
1279             TESTS::
1280
1281                 sage: set_random_seed()
1282                 sage: x = random_eja().random_element()
1283                 sage: x.subalgebra_generated_by().is_associative()
1284                 True
1285
1286             Squaring in the subalgebra should be the same thing as
1287             squaring in the superalgebra::
1288
1289                 sage: set_random_seed()
1290                 sage: x = random_eja().random_element()
1291                 sage: u = x.subalgebra_generated_by().random_element()
1292                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
1293                 True
1294
1295             """
1296             # First get the subspace spanned by the powers of myself...
1297             V = self.span_of_powers()
1298             F = self.base_ring()
1299
1300             # Now figure out the entries of the right-multiplication
1301             # matrix for the successive basis elements b0, b1,... of
1302             # that subspace.
1303             mats = []
1304             for b_right in V.basis():
1305                 eja_b_right = self.parent()(b_right)
1306                 b_right_rows = []
1307                 # The first row of the right-multiplication matrix by
1308                 # b1 is what we get if we apply that matrix to b1. The
1309                 # second row of the right multiplication matrix by b1
1310                 # is what we get when we apply that matrix to b2...
1311                 #
1312                 # IMPORTANT: this assumes that all vectors are COLUMN
1313                 # vectors, unlike our superclass (which uses row vectors).
1314                 for b_left in V.basis():
1315                     eja_b_left = self.parent()(b_left)
1316                     # Multiply in the original EJA, but then get the
1317                     # coordinates from the subalgebra in terms of its
1318                     # basis.
1319                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
1320                     b_right_rows.append(this_row)
1321                 b_right_matrix = matrix(F, b_right_rows)
1322                 mats.append(b_right_matrix)
1323
1324             # It's an algebra of polynomials in one element, and EJAs
1325             # are power-associative.
1326             #
1327             # TODO: choose generator names intelligently.
1328             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
1329
1330
1331         def subalgebra_idempotent(self):
1332             """
1333             Find an idempotent in the associative subalgebra I generate
1334             using Proposition 2.3.5 in Baes.
1335
1336             TESTS::
1337
1338                 sage: set_random_seed()
1339                 sage: J = random_eja()
1340                 sage: x = J.random_element()
1341                 sage: while x.is_nilpotent():
1342                 ....:     x = J.random_element()
1343                 sage: c = x.subalgebra_idempotent()
1344                 sage: c^2 == c
1345                 True
1346
1347             """
1348             if self.is_nilpotent():
1349                 raise ValueError("this only works with non-nilpotent elements!")
1350
1351             V = self.span_of_powers()
1352             J = self.subalgebra_generated_by()
1353             # Mis-design warning: the basis used for span_of_powers()
1354             # and subalgebra_generated_by() must be the same, and in
1355             # the same order!
1356             u = J(V.coordinates(self.vector()))
1357
1358             # The image of the matrix of left-u^m-multiplication
1359             # will be minimal for some natural number s...
1360             s = 0
1361             minimal_dim = V.dimension()
1362             for i in xrange(1, V.dimension()):
1363                 this_dim = (u**i).operator_matrix().image().dimension()
1364                 if this_dim < minimal_dim:
1365                     minimal_dim = this_dim
1366                     s = i
1367
1368             # Now minimal_matrix should correspond to the smallest
1369             # non-zero subspace in Baes's (or really, Koecher's)
1370             # proposition.
1371             #
1372             # However, we need to restrict the matrix to work on the
1373             # subspace... or do we? Can't we just solve, knowing that
1374             # A(c) = u^(s+1) should have a solution in the big space,
1375             # too?
1376             #
1377             # Beware, solve_right() means that we're using COLUMN vectors.
1378             # Our FiniteDimensionalAlgebraElement superclass uses rows.
1379             u_next = u**(s+1)
1380             A = u_next.operator_matrix()
1381             c_coordinates = A.solve_right(u_next.vector())
1382
1383             # Now c_coordinates is the idempotent we want, but it's in
1384             # the coordinate system of the subalgebra.
1385             #
1386             # We need the basis for J, but as elements of the parent algebra.
1387             #
1388             basis = [self.parent(v) for v in V.basis()]
1389             return self.parent().linear_combination(zip(c_coordinates, basis))
1390
1391
1392         def trace(self):
1393             """
1394             Return my trace, the sum of my eigenvalues.
1395
1396             EXAMPLES::
1397
1398                 sage: J = JordanSpinEJA(3)
1399                 sage: x = sum(J.gens())
1400                 sage: x.trace()
1401                 2
1402
1403             ::
1404
1405                 sage: J = RealCartesianProductEJA(5)
1406                 sage: J.one().trace()
1407                 5
1408
1409             TESTS:
1410
1411             The trace of an element is a real number::
1412
1413                 sage: set_random_seed()
1414                 sage: J = random_eja()
1415                 sage: J.random_element().trace() in J.base_ring()
1416                 True
1417
1418             """
1419             P = self.parent()
1420             r = P.rank()
1421             p = P._charpoly_coeff(r-1)
1422             # The _charpoly_coeff function already adds the factor of
1423             # -1 to ensure that _charpoly_coeff(r-1) is really what
1424             # appears in front of t^{r-1} in the charpoly. However,
1425             # we want the negative of THAT for the trace.
1426             return -p(*self.vector())
1427
1428
1429         def trace_inner_product(self, other):
1430             """
1431             Return the trace inner product of myself and ``other``.
1432
1433             TESTS:
1434
1435             The trace inner product is commutative::
1436
1437                 sage: set_random_seed()
1438                 sage: J = random_eja()
1439                 sage: x = J.random_element(); y = J.random_element()
1440                 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1441                 True
1442
1443             The trace inner product is bilinear::
1444
1445                 sage: set_random_seed()
1446                 sage: J = random_eja()
1447                 sage: x = J.random_element()
1448                 sage: y = J.random_element()
1449                 sage: z = J.random_element()
1450                 sage: a = QQ.random_element();
1451                 sage: actual = (a*(x+z)).trace_inner_product(y)
1452                 sage: expected = ( a*x.trace_inner_product(y) +
1453                 ....:              a*z.trace_inner_product(y) )
1454                 sage: actual == expected
1455                 True
1456                 sage: actual = x.trace_inner_product(a*(y+z))
1457                 sage: expected = ( a*x.trace_inner_product(y) +
1458                 ....:              a*x.trace_inner_product(z) )
1459                 sage: actual == expected
1460                 True
1461
1462             The trace inner product satisfies the compatibility
1463             condition in the definition of a Euclidean Jordan algebra::
1464
1465                 sage: set_random_seed()
1466                 sage: J = random_eja()
1467                 sage: x = J.random_element()
1468                 sage: y = J.random_element()
1469                 sage: z = J.random_element()
1470                 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1471                 True
1472
1473             """
1474             if not other in self.parent():
1475                 raise TypeError("'other' must live in the same algebra")
1476
1477             return (self*other).trace()
1478
1479
1480 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
1481     """
1482     Return the Euclidean Jordan Algebra corresponding to the set
1483     `R^n` under the Hadamard product.
1484
1485     Note: this is nothing more than the Cartesian product of ``n``
1486     copies of the spin algebra. Once Cartesian product algebras
1487     are implemented, this can go.
1488
1489     EXAMPLES:
1490
1491     This multiplication table can be verified by hand::
1492
1493         sage: J = RealCartesianProductEJA(3)
1494         sage: e0,e1,e2 = J.gens()
1495         sage: e0*e0
1496         e0
1497         sage: e0*e1
1498         0
1499         sage: e0*e2
1500         0
1501         sage: e1*e1
1502         e1
1503         sage: e1*e2
1504         0
1505         sage: e2*e2
1506         e2
1507
1508     """
1509     @staticmethod
1510     def __classcall_private__(cls, n, field=QQ):
1511         # The FiniteDimensionalAlgebra constructor takes a list of
1512         # matrices, the ith representing right multiplication by the ith
1513         # basis element in the vector space. So if e_1 = (1,0,0), then
1514         # right (Hadamard) multiplication of x by e_1 picks out the first
1515         # component of x; and likewise for the ith basis element e_i.
1516         Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i))
1517                for i in xrange(n) ]
1518
1519         fdeja = super(RealCartesianProductEJA, cls)
1520         return fdeja.__classcall_private__(cls, field, Qs, rank=n)
1521
1522     def inner_product(self, x, y):
1523         return _usual_ip(x,y)
1524
1525
1526 def random_eja():
1527     """
1528     Return a "random" finite-dimensional Euclidean Jordan Algebra.
1529
1530     ALGORITHM:
1531
1532     For now, we choose a random natural number ``n`` (greater than zero)
1533     and then give you back one of the following:
1534
1535       * The cartesian product of the rational numbers ``n`` times; this is
1536         ``QQ^n`` with the Hadamard product.
1537
1538       * The Jordan spin algebra on ``QQ^n``.
1539
1540       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1541         product.
1542
1543       * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1544         in the space of ``2n``-by-``2n`` real symmetric matrices.
1545
1546       * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1547         in the space of ``4n``-by-``4n`` real symmetric matrices.
1548
1549     Later this might be extended to return Cartesian products of the
1550     EJAs above.
1551
1552     TESTS::
1553
1554         sage: random_eja()
1555         Euclidean Jordan algebra of degree...
1556
1557     """
1558
1559     # The max_n component lets us choose different upper bounds on the
1560     # value "n" that gets passed to the constructor. This is needed
1561     # because e.g. R^{10} is reasonable to test, while the Hermitian
1562     # 10-by-10 quaternion matrices are not.
1563     (constructor, max_n) = choice([(RealCartesianProductEJA, 6),
1564                                    (JordanSpinEJA, 6),
1565                                    (RealSymmetricEJA, 5),
1566                                    (ComplexHermitianEJA, 4),
1567                                    (QuaternionHermitianEJA, 3)])
1568     n = ZZ.random_element(1, max_n)
1569     return constructor(n, field=QQ)
1570
1571
1572
1573 def _real_symmetric_basis(n, field=QQ):
1574     """
1575     Return a basis for the space of real symmetric n-by-n matrices.
1576     """
1577     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1578     # coordinates.
1579     S = []
1580     for i in xrange(n):
1581         for j in xrange(i+1):
1582             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1583             if i == j:
1584                 Sij = Eij
1585             else:
1586                 # Beware, orthogonal but not normalized!
1587                 Sij = Eij + Eij.transpose()
1588             S.append(Sij)
1589     return tuple(S)
1590
1591
1592 def _complex_hermitian_basis(n, field=QQ):
1593     """
1594     Returns a basis for the space of complex Hermitian n-by-n matrices.
1595
1596     TESTS::
1597
1598         sage: set_random_seed()
1599         sage: n = ZZ.random_element(1,5)
1600         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1601         True
1602
1603     """
1604     F = QuadraticField(-1, 'I')
1605     I = F.gen()
1606
1607     # This is like the symmetric case, but we need to be careful:
1608     #
1609     #   * We want conjugate-symmetry, not just symmetry.
1610     #   * The diagonal will (as a result) be real.
1611     #
1612     S = []
1613     for i in xrange(n):
1614         for j in xrange(i+1):
1615             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1616             if i == j:
1617                 Sij = _embed_complex_matrix(Eij)
1618                 S.append(Sij)
1619             else:
1620                 # Beware, orthogonal but not normalized! The second one
1621                 # has a minus because it's conjugated.
1622                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1623                 S.append(Sij_real)
1624                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1625                 S.append(Sij_imag)
1626     return tuple(S)
1627
1628
1629 def _quaternion_hermitian_basis(n, field=QQ):
1630     """
1631     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1632
1633     TESTS::
1634
1635         sage: set_random_seed()
1636         sage: n = ZZ.random_element(1,5)
1637         sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1638         True
1639
1640     """
1641     Q = QuaternionAlgebra(QQ,-1,-1)
1642     I,J,K = Q.gens()
1643
1644     # This is like the symmetric case, but we need to be careful:
1645     #
1646     #   * We want conjugate-symmetry, not just symmetry.
1647     #   * The diagonal will (as a result) be real.
1648     #
1649     S = []
1650     for i in xrange(n):
1651         for j in xrange(i+1):
1652             Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1653             if i == j:
1654                 Sij = _embed_quaternion_matrix(Eij)
1655                 S.append(Sij)
1656             else:
1657                 # Beware, orthogonal but not normalized! The second,
1658                 # third, and fourth ones have a minus because they're
1659                 # conjugated.
1660                 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1661                 S.append(Sij_real)
1662                 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1663                 S.append(Sij_I)
1664                 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1665                 S.append(Sij_J)
1666                 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1667                 S.append(Sij_K)
1668     return tuple(S)
1669
1670
1671 def _mat2vec(m):
1672         return vector(m.base_ring(), m.list())
1673
1674 def _vec2mat(v):
1675         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1676
1677 def _multiplication_table_from_matrix_basis(basis):
1678     """
1679     At least three of the five simple Euclidean Jordan algebras have the
1680     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1681     multiplication on the right is matrix multiplication. Given a basis
1682     for the underlying matrix space, this function returns a
1683     multiplication table (obtained by looping through the basis
1684     elements) for an algebra of those matrices. A reordered copy
1685     of the basis is also returned to work around the fact that
1686     the ``span()`` in this function will change the order of the basis
1687     from what we think it is, to... something else.
1688     """
1689     # In S^2, for example, we nominally have four coordinates even
1690     # though the space is of dimension three only. The vector space V
1691     # is supposed to hold the entire long vector, and the subspace W
1692     # of V will be spanned by the vectors that arise from symmetric
1693     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1694     field = basis[0].base_ring()
1695     dimension = basis[0].nrows()
1696
1697     V = VectorSpace(field, dimension**2)
1698     W = V.span( _mat2vec(s) for s in basis )
1699
1700     # Taking the span above reorders our basis (thanks, jerk!) so we
1701     # need to put our "matrix basis" in the same order as the
1702     # (reordered) vector basis.
1703     S = tuple( _vec2mat(b) for b in W.basis() )
1704
1705     Qs = []
1706     for s in S:
1707         # Brute force the multiplication-by-s matrix by looping
1708         # through all elements of the basis and doing the computation
1709         # to find out what the corresponding row should be. BEWARE:
1710         # these multiplication tables won't be symmetric! It therefore
1711         # becomes REALLY IMPORTANT that the underlying algebra
1712         # constructor uses ROW vectors and not COLUMN vectors. That's
1713         # why we're computing rows here and not columns.
1714         Q_rows = []
1715         for t in S:
1716             this_row = _mat2vec((s*t + t*s)/2)
1717             Q_rows.append(W.coordinates(this_row))
1718         Q = matrix(field, W.dimension(), Q_rows)
1719         Qs.append(Q)
1720
1721     return (Qs, S)
1722
1723
1724 def _embed_complex_matrix(M):
1725     """
1726     Embed the n-by-n complex matrix ``M`` into the space of real
1727     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1728     bi` to the block matrix ``[[a,b],[-b,a]]``.
1729
1730     EXAMPLES::
1731
1732         sage: F = QuadraticField(-1,'i')
1733         sage: x1 = F(4 - 2*i)
1734         sage: x2 = F(1 + 2*i)
1735         sage: x3 = F(-i)
1736         sage: x4 = F(6)
1737         sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1738         sage: _embed_complex_matrix(M)
1739         [ 4 -2| 1  2]
1740         [ 2  4|-2  1]
1741         [-----+-----]
1742         [ 0 -1| 6  0]
1743         [ 1  0| 0  6]
1744
1745     TESTS:
1746
1747     Embedding is a homomorphism (isomorphism, in fact)::
1748
1749         sage: set_random_seed()
1750         sage: n = ZZ.random_element(5)
1751         sage: F = QuadraticField(-1, 'i')
1752         sage: X = random_matrix(F, n)
1753         sage: Y = random_matrix(F, n)
1754         sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1755         sage: expected = _embed_complex_matrix(X*Y)
1756         sage: actual == expected
1757         True
1758
1759     """
1760     n = M.nrows()
1761     if M.ncols() != n:
1762         raise ValueError("the matrix 'M' must be square")
1763     field = M.base_ring()
1764     blocks = []
1765     for z in M.list():
1766         a = z.real()
1767         b = z.imag()
1768         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1769
1770     # We can drop the imaginaries here.
1771     return block_matrix(field.base_ring(), n, blocks)
1772
1773
1774 def _unembed_complex_matrix(M):
1775     """
1776     The inverse of _embed_complex_matrix().
1777
1778     EXAMPLES::
1779
1780         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1781         ....:                 [-2,  1,  -4,  3],
1782         ....:                 [ 9,  10, 11, 12],
1783         ....:                 [-10, 9, -12, 11] ])
1784         sage: _unembed_complex_matrix(A)
1785         [  2*i + 1   4*i + 3]
1786         [ 10*i + 9 12*i + 11]
1787
1788     TESTS:
1789
1790     Unembedding is the inverse of embedding::
1791
1792         sage: set_random_seed()
1793         sage: F = QuadraticField(-1, 'i')
1794         sage: M = random_matrix(F, 3)
1795         sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1796         True
1797
1798     """
1799     n = ZZ(M.nrows())
1800     if M.ncols() != n:
1801         raise ValueError("the matrix 'M' must be square")
1802     if not n.mod(2).is_zero():
1803         raise ValueError("the matrix 'M' must be a complex embedding")
1804
1805     F = QuadraticField(-1, 'i')
1806     i = F.gen()
1807
1808     # Go top-left to bottom-right (reading order), converting every
1809     # 2-by-2 block we see to a single complex element.
1810     elements = []
1811     for k in xrange(n/2):
1812         for j in xrange(n/2):
1813             submat = M[2*k:2*k+2,2*j:2*j+2]
1814             if submat[0,0] != submat[1,1]:
1815                 raise ValueError('bad on-diagonal submatrix')
1816             if submat[0,1] != -submat[1,0]:
1817                 raise ValueError('bad off-diagonal submatrix')
1818             z = submat[0,0] + submat[0,1]*i
1819             elements.append(z)
1820
1821     return matrix(F, n/2, elements)
1822
1823
1824 def _embed_quaternion_matrix(M):
1825     """
1826     Embed the n-by-n quaternion matrix ``M`` into the space of real
1827     matrices of size 4n-by-4n by first sending each quaternion entry
1828     `z = a + bi + cj + dk` to the block-complex matrix
1829     ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1830     a real matrix.
1831
1832     EXAMPLES::
1833
1834         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1835         sage: i,j,k = Q.gens()
1836         sage: x = 1 + 2*i + 3*j + 4*k
1837         sage: M = matrix(Q, 1, [[x]])
1838         sage: _embed_quaternion_matrix(M)
1839         [ 1  2  3  4]
1840         [-2  1 -4  3]
1841         [-3  4  1 -2]
1842         [-4 -3  2  1]
1843
1844     Embedding is a homomorphism (isomorphism, in fact)::
1845
1846         sage: set_random_seed()
1847         sage: n = ZZ.random_element(5)
1848         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1849         sage: X = random_matrix(Q, n)
1850         sage: Y = random_matrix(Q, n)
1851         sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1852         sage: expected = _embed_quaternion_matrix(X*Y)
1853         sage: actual == expected
1854         True
1855
1856     """
1857     quaternions = M.base_ring()
1858     n = M.nrows()
1859     if M.ncols() != n:
1860         raise ValueError("the matrix 'M' must be square")
1861
1862     F = QuadraticField(-1, 'i')
1863     i = F.gen()
1864
1865     blocks = []
1866     for z in M.list():
1867         t = z.coefficient_tuple()
1868         a = t[0]
1869         b = t[1]
1870         c = t[2]
1871         d = t[3]
1872         cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1873                                     [-c + d*i, a - b*i]])
1874         blocks.append(_embed_complex_matrix(cplx_matrix))
1875
1876     # We should have real entries by now, so use the realest field
1877     # we've got for the return value.
1878     return block_matrix(quaternions.base_ring(), n, blocks)
1879
1880
1881 def _unembed_quaternion_matrix(M):
1882     """
1883     The inverse of _embed_quaternion_matrix().
1884
1885     EXAMPLES::
1886
1887         sage: M = matrix(QQ, [[ 1,  2,  3,  4],
1888         ....:                 [-2,  1, -4,  3],
1889         ....:                 [-3,  4,  1, -2],
1890         ....:                 [-4, -3,  2,  1]])
1891         sage: _unembed_quaternion_matrix(M)
1892         [1 + 2*i + 3*j + 4*k]
1893
1894     TESTS:
1895
1896     Unembedding is the inverse of embedding::
1897
1898         sage: set_random_seed()
1899         sage: Q = QuaternionAlgebra(QQ, -1, -1)
1900         sage: M = random_matrix(Q, 3)
1901         sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1902         True
1903
1904     """
1905     n = ZZ(M.nrows())
1906     if M.ncols() != n:
1907         raise ValueError("the matrix 'M' must be square")
1908     if not n.mod(4).is_zero():
1909         raise ValueError("the matrix 'M' must be a complex embedding")
1910
1911     Q = QuaternionAlgebra(QQ,-1,-1)
1912     i,j,k = Q.gens()
1913
1914     # Go top-left to bottom-right (reading order), converting every
1915     # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1916     # quaternion block.
1917     elements = []
1918     for l in xrange(n/4):
1919         for m in xrange(n/4):
1920             submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1921             if submat[0,0] != submat[1,1].conjugate():
1922                 raise ValueError('bad on-diagonal submatrix')
1923             if submat[0,1] != -submat[1,0].conjugate():
1924                 raise ValueError('bad off-diagonal submatrix')
1925             z  = submat[0,0].real() + submat[0,0].imag()*i
1926             z += submat[0,1].real()*j + submat[0,1].imag()*k
1927             elements.append(z)
1928
1929     return matrix(Q, n/4, elements)
1930
1931
1932 # The usual inner product on R^n.
1933 def _usual_ip(x,y):
1934     return x.vector().inner_product(y.vector())
1935
1936 # The inner product used for the real symmetric simple EJA.
1937 # We keep it as a separate function because e.g. the complex
1938 # algebra uses the same inner product, except divided by 2.
1939 def _matrix_ip(X,Y):
1940     X_mat = X.natural_representation()
1941     Y_mat = Y.natural_representation()
1942     return (X_mat*Y_mat).trace()
1943
1944
1945 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
1946     """
1947     The rank-n simple EJA consisting of real symmetric n-by-n
1948     matrices, the usual symmetric Jordan product, and the trace inner
1949     product. It has dimension `(n^2 + n)/2` over the reals.
1950
1951     EXAMPLES::
1952
1953         sage: J = RealSymmetricEJA(2)
1954         sage: e0, e1, e2 = J.gens()
1955         sage: e0*e0
1956         e0
1957         sage: e1*e1
1958         e0 + e2
1959         sage: e2*e2
1960         e2
1961
1962     TESTS:
1963
1964     The degree of this algebra is `(n^2 + n) / 2`::
1965
1966         sage: set_random_seed()
1967         sage: n = ZZ.random_element(1,5)
1968         sage: J = RealSymmetricEJA(n)
1969         sage: J.degree() == (n^2 + n)/2
1970         True
1971
1972     The Jordan multiplication is what we think it is::
1973
1974         sage: set_random_seed()
1975         sage: n = ZZ.random_element(1,5)
1976         sage: J = RealSymmetricEJA(n)
1977         sage: x = J.random_element()
1978         sage: y = J.random_element()
1979         sage: actual = (x*y).natural_representation()
1980         sage: X = x.natural_representation()
1981         sage: Y = y.natural_representation()
1982         sage: expected = (X*Y + Y*X)/2
1983         sage: actual == expected
1984         True
1985         sage: J(expected) == x*y
1986         True
1987
1988     """
1989     @staticmethod
1990     def __classcall_private__(cls, n, field=QQ):
1991         S = _real_symmetric_basis(n, field=field)
1992         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1993
1994         fdeja = super(RealSymmetricEJA, cls)
1995         return fdeja.__classcall_private__(cls,
1996                                            field,
1997                                            Qs,
1998                                            rank=n,
1999                                            natural_basis=T)
2000
2001     def inner_product(self, x, y):
2002         return _matrix_ip(x,y)
2003
2004
2005 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
2006     """
2007     The rank-n simple EJA consisting of complex Hermitian n-by-n
2008     matrices over the real numbers, the usual symmetric Jordan product,
2009     and the real-part-of-trace inner product. It has dimension `n^2` over
2010     the reals.
2011
2012     TESTS:
2013
2014     The degree of this algebra is `n^2`::
2015
2016         sage: set_random_seed()
2017         sage: n = ZZ.random_element(1,5)
2018         sage: J = ComplexHermitianEJA(n)
2019         sage: J.degree() == n^2
2020         True
2021
2022     The Jordan multiplication is what we think it is::
2023
2024         sage: set_random_seed()
2025         sage: n = ZZ.random_element(1,5)
2026         sage: J = ComplexHermitianEJA(n)
2027         sage: x = J.random_element()
2028         sage: y = J.random_element()
2029         sage: actual = (x*y).natural_representation()
2030         sage: X = x.natural_representation()
2031         sage: Y = y.natural_representation()
2032         sage: expected = (X*Y + Y*X)/2
2033         sage: actual == expected
2034         True
2035         sage: J(expected) == x*y
2036         True
2037
2038     """
2039     @staticmethod
2040     def __classcall_private__(cls, n, field=QQ):
2041         S = _complex_hermitian_basis(n)
2042         (Qs, T) = _multiplication_table_from_matrix_basis(S)
2043
2044         fdeja = super(ComplexHermitianEJA, cls)
2045         return fdeja.__classcall_private__(cls,
2046                                            field,
2047                                            Qs,
2048                                            rank=n,
2049                                            natural_basis=T)
2050
2051     def inner_product(self, x, y):
2052         # Since a+bi on the diagonal is represented as
2053         #
2054         #   a + bi  = [  a  b  ]
2055         #             [ -b  a  ],
2056         #
2057         # we'll double-count the "a" entries if we take the trace of
2058         # the embedding.
2059         return _matrix_ip(x,y)/2
2060
2061
2062 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
2063     """
2064     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2065     matrices, the usual symmetric Jordan product, and the
2066     real-part-of-trace inner product. It has dimension `2n^2 - n` over
2067     the reals.
2068
2069     TESTS:
2070
2071     The degree of this algebra is `n^2`::
2072
2073         sage: set_random_seed()
2074         sage: n = ZZ.random_element(1,5)
2075         sage: J = QuaternionHermitianEJA(n)
2076         sage: J.degree() == 2*(n^2) - n
2077         True
2078
2079     The Jordan multiplication is what we think it is::
2080
2081         sage: set_random_seed()
2082         sage: n = ZZ.random_element(1,5)
2083         sage: J = QuaternionHermitianEJA(n)
2084         sage: x = J.random_element()
2085         sage: y = J.random_element()
2086         sage: actual = (x*y).natural_representation()
2087         sage: X = x.natural_representation()
2088         sage: Y = y.natural_representation()
2089         sage: expected = (X*Y + Y*X)/2
2090         sage: actual == expected
2091         True
2092         sage: J(expected) == x*y
2093         True
2094
2095     """
2096     @staticmethod
2097     def __classcall_private__(cls, n, field=QQ):
2098         S = _quaternion_hermitian_basis(n)
2099         (Qs, T) = _multiplication_table_from_matrix_basis(S)
2100
2101         fdeja = super(QuaternionHermitianEJA, cls)
2102         return fdeja.__classcall_private__(cls,
2103                                            field,
2104                                            Qs,
2105                                            rank=n,
2106                                            natural_basis=T)
2107
2108     def inner_product(self, x, y):
2109         # Since a+bi+cj+dk on the diagonal is represented as
2110         #
2111         #   a + bi +cj + dk = [  a  b  c  d]
2112         #                     [ -b  a -d  c]
2113         #                     [ -c  d  a -b]
2114         #                     [ -d -c  b  a],
2115         #
2116         # we'll quadruple-count the "a" entries if we take the trace of
2117         # the embedding.
2118         return _matrix_ip(x,y)/4
2119
2120
2121 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
2122     """
2123     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2124     with the usual inner product and jordan product ``x*y =
2125     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2126     the reals.
2127
2128     EXAMPLES:
2129
2130     This multiplication table can be verified by hand::
2131
2132         sage: J = JordanSpinEJA(4)
2133         sage: e0,e1,e2,e3 = J.gens()
2134         sage: e0*e0
2135         e0
2136         sage: e0*e1
2137         e1
2138         sage: e0*e2
2139         e2
2140         sage: e0*e3
2141         e3
2142         sage: e1*e2
2143         0
2144         sage: e1*e3
2145         0
2146         sage: e2*e3
2147         0
2148
2149     """
2150     @staticmethod
2151     def __classcall_private__(cls, n, field=QQ):
2152         Qs = []
2153         id_matrix = identity_matrix(field, n)
2154         for i in xrange(n):
2155             ei = id_matrix.column(i)
2156             Qi = zero_matrix(field, n)
2157             Qi.set_row(0, ei)
2158             Qi.set_column(0, ei)
2159             Qi += diagonal_matrix(n, [ei[0]]*n)
2160             # The addition of the diagonal matrix adds an extra ei[0] in the
2161             # upper-left corner of the matrix.
2162             Qi[0,0] = Qi[0,0] * ~field(2)
2163             Qs.append(Qi)
2164
2165         # The rank of the spin algebra is two, unless we're in a
2166         # one-dimensional ambient space (because the rank is bounded by
2167         # the ambient dimension).
2168         fdeja = super(JordanSpinEJA, cls)
2169         return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2))
2170
2171     def inner_product(self, x, y):
2172         return _usual_ip(x,y)