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(self):
1078             """
1079             Return the left-multiplication-by-this-element
1080             operator on the ambient algebra.
1081
1082             TESTS::
1083
1084                 sage: set_random_seed()
1085                 sage: J = random_eja()
1086                 sage: x = J.random_element()
1087                 sage: y = J.random_element()
1088                 sage: x.operator()(y) == x*y
1089                 True
1090                 sage: y.operator()(x) == x*y
1091                 True
1092
1093             """
1094             P = self.parent()
1095             return FiniteDimensionalEuclideanJordanAlgebraMorphism(
1096                      Hom(P,P),
1097                      self.operator_matrix() )
1098
1099
1100
1101         def operator_matrix(self):
1102             """
1103             Return the matrix that represents left- (or right-)
1104             multiplication by this element in the parent algebra.
1105
1106             We implement this ourselves to work around the fact that
1107             our parent class represents everything with row vectors.
1108
1109             EXAMPLES:
1110
1111             Test the first polarization identity from my notes, Koecher Chapter
1112             III, or from Baes (2.3)::
1113
1114                 sage: set_random_seed()
1115                 sage: J = random_eja()
1116                 sage: x = J.random_element()
1117                 sage: y = J.random_element()
1118                 sage: Lx = x.operator_matrix()
1119                 sage: Ly = y.operator_matrix()
1120                 sage: Lxx = (x*x).operator_matrix()
1121                 sage: Lxy = (x*y).operator_matrix()
1122                 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
1123                 True
1124
1125             Test the second polarization identity from my notes or from
1126             Baes (2.4)::
1127
1128                 sage: set_random_seed()
1129                 sage: J = random_eja()
1130                 sage: x = J.random_element()
1131                 sage: y = J.random_element()
1132                 sage: z = J.random_element()
1133                 sage: Lx = x.operator_matrix()
1134                 sage: Ly = y.operator_matrix()
1135                 sage: Lz = z.operator_matrix()
1136                 sage: Lzy = (z*y).operator_matrix()
1137                 sage: Lxy = (x*y).operator_matrix()
1138                 sage: Lxz = (x*z).operator_matrix()
1139                 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
1140                 True
1141
1142             Test the third polarization identity from my notes or from
1143             Baes (2.5)::
1144
1145                 sage: set_random_seed()
1146                 sage: J = random_eja()
1147                 sage: u = J.random_element()
1148                 sage: y = J.random_element()
1149                 sage: z = J.random_element()
1150                 sage: Lu = u.operator_matrix()
1151                 sage: Ly = y.operator_matrix()
1152                 sage: Lz = z.operator_matrix()
1153                 sage: Lzy = (z*y).operator_matrix()
1154                 sage: Luy = (u*y).operator_matrix()
1155                 sage: Luz = (u*z).operator_matrix()
1156                 sage: Luyz = (u*(y*z)).operator_matrix()
1157                 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
1158                 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
1159                 sage: bool(lhs == rhs)
1160                 True
1161
1162             """
1163             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
1164             return fda_elt.matrix().transpose()
1165
1166
1167         def quadratic_representation(self, other=None):
1168             """
1169             Return the quadratic representation of this element.
1170
1171             EXAMPLES:
1172
1173             The explicit form in the spin factor algebra is given by
1174             Alizadeh's Example 11.12::
1175
1176                 sage: set_random_seed()
1177                 sage: n = ZZ.random_element(1,10)
1178                 sage: J = JordanSpinEJA(n)
1179                 sage: x = J.random_element()
1180                 sage: x_vec = x.vector()
1181                 sage: x0 = x_vec[0]
1182                 sage: x_bar = x_vec[1:]
1183                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1184                 sage: B = 2*x0*x_bar.row()
1185                 sage: C = 2*x0*x_bar.column()
1186                 sage: D = identity_matrix(QQ, n-1)
1187                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1188                 sage: D = D + 2*x_bar.tensor_product(x_bar)
1189                 sage: Q = block_matrix(2,2,[A,B,C,D])
1190                 sage: Q == x.quadratic_representation()
1191                 True
1192
1193             Test all of the properties from Theorem 11.2 in Alizadeh::
1194
1195                 sage: set_random_seed()
1196                 sage: J = random_eja()
1197                 sage: x = J.random_element()
1198                 sage: y = J.random_element()
1199                 sage: Lx = x.operator_matrix()
1200                 sage: Lxx = (x*x).operator_matrix()
1201                 sage: Qx = x.quadratic_representation()
1202                 sage: Qy = y.quadratic_representation()
1203                 sage: Qxy = x.quadratic_representation(y)
1204                 sage: Qex = J.one().quadratic_representation(x)
1205                 sage: n = ZZ.random_element(10)
1206                 sage: Qxn = (x^n).quadratic_representation()
1207
1208             Property 1:
1209
1210                 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1211                 True
1212
1213             Property 2:
1214
1215                 sage: alpha = QQ.random_element()
1216                 sage: (alpha*x).quadratic_representation() == (alpha^2)*Qx
1217                 True
1218
1219             Property 3:
1220
1221                 sage: not x.is_invertible() or (
1222                 ....:     Qx*x.inverse().vector() == x.vector() )
1223                 True
1224
1225                 sage: not x.is_invertible() or (
1226                 ....:   Qx.inverse()
1227                 ....:   ==
1228                 ....:   x.inverse().quadratic_representation() )
1229                 True
1230
1231                 sage: Qxy*(J.one().vector()) == (x*y).vector()
1232                 True
1233
1234             Property 4:
1235
1236                 sage: not x.is_invertible() or (
1237                 ....:   x.quadratic_representation(x.inverse())*Qx
1238                 ....:   == Qx*x.quadratic_representation(x.inverse()) )
1239                 True
1240
1241                 sage: not x.is_invertible() or (
1242                 ....:   x.quadratic_representation(x.inverse())*Qx
1243                 ....:   ==
1244                 ....:   2*x.operator_matrix()*Qex - Qx )
1245                 True
1246
1247                 sage: 2*x.operator_matrix()*Qex - Qx == Lxx
1248                 True
1249
1250             Property 5:
1251
1252                 sage: J(Qy*x.vector()).quadratic_representation() == Qy*Qx*Qy
1253                 True
1254
1255             Property 6:
1256
1257                 sage: Qxn == (Qx)^n
1258                 True
1259
1260             Property 7:
1261
1262                 sage: not x.is_invertible() or (
1263                 ....:   Qx*x.inverse().operator_matrix() == Lx )
1264                 True
1265
1266             Property 8:
1267
1268                 sage: not x.operator_commutes_with(y) or (
1269                 ....:   J(Qx*y.vector())^n == J(Qxn*(y^n).vector()) )
1270                 True
1271
1272             """
1273             if other is None:
1274                 other=self
1275             elif not other in self.parent():
1276                 raise TypeError("'other' must live in the same algebra")
1277
1278             L = self.operator_matrix()
1279             M = other.operator_matrix()
1280             return ( L*M + M*L - (self*other).operator_matrix() )
1281
1282
1283         def span_of_powers(self):
1284             """
1285             Return the vector space spanned by successive powers of
1286             this element.
1287             """
1288             # The dimension of the subalgebra can't be greater than
1289             # the big algebra, so just put everything into a list
1290             # and let span() get rid of the excess.
1291             #
1292             # We do the extra ambient_vector_space() in case we're messing
1293             # with polynomials and the direct parent is a module.
1294             V = self.vector().parent().ambient_vector_space()
1295             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
1296
1297
1298         def subalgebra_generated_by(self):
1299             """
1300             Return the associative subalgebra of the parent EJA generated
1301             by this element.
1302
1303             TESTS::
1304
1305                 sage: set_random_seed()
1306                 sage: x = random_eja().random_element()
1307                 sage: x.subalgebra_generated_by().is_associative()
1308                 True
1309
1310             Squaring in the subalgebra should be the same thing as
1311             squaring in the superalgebra::
1312
1313                 sage: set_random_seed()
1314                 sage: x = random_eja().random_element()
1315                 sage: u = x.subalgebra_generated_by().random_element()
1316                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
1317                 True
1318
1319             """
1320             # First get the subspace spanned by the powers of myself...
1321             V = self.span_of_powers()
1322             F = self.base_ring()
1323
1324             # Now figure out the entries of the right-multiplication
1325             # matrix for the successive basis elements b0, b1,... of
1326             # that subspace.
1327             mats = []
1328             for b_right in V.basis():
1329                 eja_b_right = self.parent()(b_right)
1330                 b_right_rows = []
1331                 # The first row of the right-multiplication matrix by
1332                 # b1 is what we get if we apply that matrix to b1. The
1333                 # second row of the right multiplication matrix by b1
1334                 # is what we get when we apply that matrix to b2...
1335                 #
1336                 # IMPORTANT: this assumes that all vectors are COLUMN
1337                 # vectors, unlike our superclass (which uses row vectors).
1338                 for b_left in V.basis():
1339                     eja_b_left = self.parent()(b_left)
1340                     # Multiply in the original EJA, but then get the
1341                     # coordinates from the subalgebra in terms of its
1342                     # basis.
1343                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
1344                     b_right_rows.append(this_row)
1345                 b_right_matrix = matrix(F, b_right_rows)
1346                 mats.append(b_right_matrix)
1347
1348             # It's an algebra of polynomials in one element, and EJAs
1349             # are power-associative.
1350             #
1351             # TODO: choose generator names intelligently.
1352             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
1353
1354
1355         def subalgebra_idempotent(self):
1356             """
1357             Find an idempotent in the associative subalgebra I generate
1358             using Proposition 2.3.5 in Baes.
1359
1360             TESTS::
1361
1362                 sage: set_random_seed()
1363                 sage: J = random_eja()
1364                 sage: x = J.random_element()
1365                 sage: while x.is_nilpotent():
1366                 ....:     x = J.random_element()
1367                 sage: c = x.subalgebra_idempotent()
1368                 sage: c^2 == c
1369                 True
1370
1371             """
1372             if self.is_nilpotent():
1373                 raise ValueError("this only works with non-nilpotent elements!")
1374
1375             V = self.span_of_powers()
1376             J = self.subalgebra_generated_by()
1377             # Mis-design warning: the basis used for span_of_powers()
1378             # and subalgebra_generated_by() must be the same, and in
1379             # the same order!
1380             u = J(V.coordinates(self.vector()))
1381
1382             # The image of the matrix of left-u^m-multiplication
1383             # will be minimal for some natural number s...
1384             s = 0
1385             minimal_dim = V.dimension()
1386             for i in xrange(1, V.dimension()):
1387                 this_dim = (u**i).operator_matrix().image().dimension()
1388                 if this_dim < minimal_dim:
1389                     minimal_dim = this_dim
1390                     s = i
1391
1392             # Now minimal_matrix should correspond to the smallest
1393             # non-zero subspace in Baes's (or really, Koecher's)
1394             # proposition.
1395             #
1396             # However, we need to restrict the matrix to work on the
1397             # subspace... or do we? Can't we just solve, knowing that
1398             # A(c) = u^(s+1) should have a solution in the big space,
1399             # too?
1400             #
1401             # Beware, solve_right() means that we're using COLUMN vectors.
1402             # Our FiniteDimensionalAlgebraElement superclass uses rows.
1403             u_next = u**(s+1)
1404             A = u_next.operator_matrix()
1405             c_coordinates = A.solve_right(u_next.vector())
1406
1407             # Now c_coordinates is the idempotent we want, but it's in
1408             # the coordinate system of the subalgebra.
1409             #
1410             # We need the basis for J, but as elements of the parent algebra.
1411             #
1412             basis = [self.parent(v) for v in V.basis()]
1413             return self.parent().linear_combination(zip(c_coordinates, basis))
1414
1415
1416         def trace(self):
1417             """
1418             Return my trace, the sum of my eigenvalues.
1419
1420             EXAMPLES::
1421
1422                 sage: J = JordanSpinEJA(3)
1423                 sage: x = sum(J.gens())
1424                 sage: x.trace()
1425                 2
1426
1427             ::
1428
1429                 sage: J = RealCartesianProductEJA(5)
1430                 sage: J.one().trace()
1431                 5
1432
1433             TESTS:
1434
1435             The trace of an element is a real number::
1436
1437                 sage: set_random_seed()
1438                 sage: J = random_eja()
1439                 sage: J.random_element().trace() in J.base_ring()
1440                 True
1441
1442             """
1443             P = self.parent()
1444             r = P.rank()
1445             p = P._charpoly_coeff(r-1)
1446             # The _charpoly_coeff function already adds the factor of
1447             # -1 to ensure that _charpoly_coeff(r-1) is really what
1448             # appears in front of t^{r-1} in the charpoly. However,
1449             # we want the negative of THAT for the trace.
1450             return -p(*self.vector())
1451
1452
1453         def trace_inner_product(self, other):
1454             """
1455             Return the trace inner product of myself and ``other``.
1456
1457             TESTS:
1458
1459             The trace inner product is commutative::
1460
1461                 sage: set_random_seed()
1462                 sage: J = random_eja()
1463                 sage: x = J.random_element(); y = J.random_element()
1464                 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1465                 True
1466
1467             The trace inner product is bilinear::
1468
1469                 sage: set_random_seed()
1470                 sage: J = random_eja()
1471                 sage: x = J.random_element()
1472                 sage: y = J.random_element()
1473                 sage: z = J.random_element()
1474                 sage: a = QQ.random_element();
1475                 sage: actual = (a*(x+z)).trace_inner_product(y)
1476                 sage: expected = ( a*x.trace_inner_product(y) +
1477                 ....:              a*z.trace_inner_product(y) )
1478                 sage: actual == expected
1479                 True
1480                 sage: actual = x.trace_inner_product(a*(y+z))
1481                 sage: expected = ( a*x.trace_inner_product(y) +
1482                 ....:              a*x.trace_inner_product(z) )
1483                 sage: actual == expected
1484                 True
1485
1486             The trace inner product satisfies the compatibility
1487             condition in the definition of a Euclidean Jordan algebra::
1488
1489                 sage: set_random_seed()
1490                 sage: J = random_eja()
1491                 sage: x = J.random_element()
1492                 sage: y = J.random_element()
1493                 sage: z = J.random_element()
1494                 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1495                 True
1496
1497             """
1498             if not other in self.parent():
1499                 raise TypeError("'other' must live in the same algebra")
1500
1501             return (self*other).trace()
1502
1503
1504 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
1505     """
1506     Return the Euclidean Jordan Algebra corresponding to the set
1507     `R^n` under the Hadamard product.
1508
1509     Note: this is nothing more than the Cartesian product of ``n``
1510     copies of the spin algebra. Once Cartesian product algebras
1511     are implemented, this can go.
1512
1513     EXAMPLES:
1514
1515     This multiplication table can be verified by hand::
1516
1517         sage: J = RealCartesianProductEJA(3)
1518         sage: e0,e1,e2 = J.gens()
1519         sage: e0*e0
1520         e0
1521         sage: e0*e1
1522         0
1523         sage: e0*e2
1524         0
1525         sage: e1*e1
1526         e1
1527         sage: e1*e2
1528         0
1529         sage: e2*e2
1530         e2
1531
1532     """
1533     @staticmethod
1534     def __classcall_private__(cls, n, field=QQ):
1535         # The FiniteDimensionalAlgebra constructor takes a list of
1536         # matrices, the ith representing right multiplication by the ith
1537         # basis element in the vector space. So if e_1 = (1,0,0), then
1538         # right (Hadamard) multiplication of x by e_1 picks out the first
1539         # component of x; and likewise for the ith basis element e_i.
1540         Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i))
1541                for i in xrange(n) ]
1542
1543         fdeja = super(RealCartesianProductEJA, cls)
1544         return fdeja.__classcall_private__(cls, field, Qs, rank=n)
1545
1546     def inner_product(self, x, y):
1547         return _usual_ip(x,y)
1548
1549
1550 def random_eja():
1551     """
1552     Return a "random" finite-dimensional Euclidean Jordan Algebra.
1553
1554     ALGORITHM:
1555
1556     For now, we choose a random natural number ``n`` (greater than zero)
1557     and then give you back one of the following:
1558
1559       * The cartesian product of the rational numbers ``n`` times; this is
1560         ``QQ^n`` with the Hadamard product.
1561
1562       * The Jordan spin algebra on ``QQ^n``.
1563
1564       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1565         product.
1566
1567       * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1568         in the space of ``2n``-by-``2n`` real symmetric matrices.
1569
1570       * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1571         in the space of ``4n``-by-``4n`` real symmetric matrices.
1572
1573     Later this might be extended to return Cartesian products of the
1574     EJAs above.
1575
1576     TESTS::
1577
1578         sage: random_eja()
1579         Euclidean Jordan algebra of degree...
1580
1581     """
1582
1583     # The max_n component lets us choose different upper bounds on the
1584     # value "n" that gets passed to the constructor. This is needed
1585     # because e.g. R^{10} is reasonable to test, while the Hermitian
1586     # 10-by-10 quaternion matrices are not.
1587     (constructor, max_n) = choice([(RealCartesianProductEJA, 6),
1588                                    (JordanSpinEJA, 6),
1589                                    (RealSymmetricEJA, 5),
1590                                    (ComplexHermitianEJA, 4),
1591                                    (QuaternionHermitianEJA, 3)])
1592     n = ZZ.random_element(1, max_n)
1593     return constructor(n, field=QQ)
1594
1595
1596
1597 def _real_symmetric_basis(n, field=QQ):
1598     """
1599     Return a basis for the space of real symmetric n-by-n matrices.
1600     """
1601     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1602     # coordinates.
1603     S = []
1604     for i in xrange(n):
1605         for j in xrange(i+1):
1606             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1607             if i == j:
1608                 Sij = Eij
1609             else:
1610                 # Beware, orthogonal but not normalized!
1611                 Sij = Eij + Eij.transpose()
1612             S.append(Sij)
1613     return tuple(S)
1614
1615
1616 def _complex_hermitian_basis(n, field=QQ):
1617     """
1618     Returns a basis for the space of complex Hermitian n-by-n matrices.
1619
1620     TESTS::
1621
1622         sage: set_random_seed()
1623         sage: n = ZZ.random_element(1,5)
1624         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1625         True
1626
1627     """
1628     F = QuadraticField(-1, 'I')
1629     I = F.gen()
1630
1631     # This is like the symmetric case, but we need to be careful:
1632     #
1633     #   * We want conjugate-symmetry, not just symmetry.
1634     #   * The diagonal will (as a result) be real.
1635     #
1636     S = []
1637     for i in xrange(n):
1638         for j in xrange(i+1):
1639             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1640             if i == j:
1641                 Sij = _embed_complex_matrix(Eij)
1642                 S.append(Sij)
1643             else:
1644                 # Beware, orthogonal but not normalized! The second one
1645                 # has a minus because it's conjugated.
1646                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1647                 S.append(Sij_real)
1648                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1649                 S.append(Sij_imag)
1650     return tuple(S)
1651
1652
1653 def _quaternion_hermitian_basis(n, field=QQ):
1654     """
1655     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1656
1657     TESTS::
1658
1659         sage: set_random_seed()
1660         sage: n = ZZ.random_element(1,5)
1661         sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1662         True
1663
1664     """
1665     Q = QuaternionAlgebra(QQ,-1,-1)
1666     I,J,K = Q.gens()
1667
1668     # This is like the symmetric case, but we need to be careful:
1669     #
1670     #   * We want conjugate-symmetry, not just symmetry.
1671     #   * The diagonal will (as a result) be real.
1672     #
1673     S = []
1674     for i in xrange(n):
1675         for j in xrange(i+1):
1676             Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1677             if i == j:
1678                 Sij = _embed_quaternion_matrix(Eij)
1679                 S.append(Sij)
1680             else:
1681                 # Beware, orthogonal but not normalized! The second,
1682                 # third, and fourth ones have a minus because they're
1683                 # conjugated.
1684                 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1685                 S.append(Sij_real)
1686                 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1687                 S.append(Sij_I)
1688                 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1689                 S.append(Sij_J)
1690                 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1691                 S.append(Sij_K)
1692     return tuple(S)
1693
1694
1695 def _mat2vec(m):
1696         return vector(m.base_ring(), m.list())
1697
1698 def _vec2mat(v):
1699         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1700
1701 def _multiplication_table_from_matrix_basis(basis):
1702     """
1703     At least three of the five simple Euclidean Jordan algebras have the
1704     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1705     multiplication on the right is matrix multiplication. Given a basis
1706     for the underlying matrix space, this function returns a
1707     multiplication table (obtained by looping through the basis
1708     elements) for an algebra of those matrices. A reordered copy
1709     of the basis is also returned to work around the fact that
1710     the ``span()`` in this function will change the order of the basis
1711     from what we think it is, to... something else.
1712     """
1713     # In S^2, for example, we nominally have four coordinates even
1714     # though the space is of dimension three only. The vector space V
1715     # is supposed to hold the entire long vector, and the subspace W
1716     # of V will be spanned by the vectors that arise from symmetric
1717     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1718     field = basis[0].base_ring()
1719     dimension = basis[0].nrows()
1720
1721     V = VectorSpace(field, dimension**2)
1722     W = V.span( _mat2vec(s) for s in basis )
1723
1724     # Taking the span above reorders our basis (thanks, jerk!) so we
1725     # need to put our "matrix basis" in the same order as the
1726     # (reordered) vector basis.
1727     S = tuple( _vec2mat(b) for b in W.basis() )
1728
1729     Qs = []
1730     for s in S:
1731         # Brute force the multiplication-by-s matrix by looping
1732         # through all elements of the basis and doing the computation
1733         # to find out what the corresponding row should be. BEWARE:
1734         # these multiplication tables won't be symmetric! It therefore
1735         # becomes REALLY IMPORTANT that the underlying algebra
1736         # constructor uses ROW vectors and not COLUMN vectors. That's
1737         # why we're computing rows here and not columns.
1738         Q_rows = []
1739         for t in S:
1740             this_row = _mat2vec((s*t + t*s)/2)
1741             Q_rows.append(W.coordinates(this_row))
1742         Q = matrix(field, W.dimension(), Q_rows)
1743         Qs.append(Q)
1744
1745     return (Qs, S)
1746
1747
1748 def _embed_complex_matrix(M):
1749     """
1750     Embed the n-by-n complex matrix ``M`` into the space of real
1751     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1752     bi` to the block matrix ``[[a,b],[-b,a]]``.
1753
1754     EXAMPLES::
1755
1756         sage: F = QuadraticField(-1,'i')
1757         sage: x1 = F(4 - 2*i)
1758         sage: x2 = F(1 + 2*i)
1759         sage: x3 = F(-i)
1760         sage: x4 = F(6)
1761         sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1762         sage: _embed_complex_matrix(M)
1763         [ 4 -2| 1  2]
1764         [ 2  4|-2  1]
1765         [-----+-----]
1766         [ 0 -1| 6  0]
1767         [ 1  0| 0  6]
1768
1769     TESTS:
1770
1771     Embedding is a homomorphism (isomorphism, in fact)::
1772
1773         sage: set_random_seed()
1774         sage: n = ZZ.random_element(5)
1775         sage: F = QuadraticField(-1, 'i')
1776         sage: X = random_matrix(F, n)
1777         sage: Y = random_matrix(F, n)
1778         sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1779         sage: expected = _embed_complex_matrix(X*Y)
1780         sage: actual == expected
1781         True
1782
1783     """
1784     n = M.nrows()
1785     if M.ncols() != n:
1786         raise ValueError("the matrix 'M' must be square")
1787     field = M.base_ring()
1788     blocks = []
1789     for z in M.list():
1790         a = z.real()
1791         b = z.imag()
1792         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1793
1794     # We can drop the imaginaries here.
1795     return block_matrix(field.base_ring(), n, blocks)
1796
1797
1798 def _unembed_complex_matrix(M):
1799     """
1800     The inverse of _embed_complex_matrix().
1801
1802     EXAMPLES::
1803
1804         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1805         ....:                 [-2,  1,  -4,  3],
1806         ....:                 [ 9,  10, 11, 12],
1807         ....:                 [-10, 9, -12, 11] ])
1808         sage: _unembed_complex_matrix(A)
1809         [  2*i + 1   4*i + 3]
1810         [ 10*i + 9 12*i + 11]
1811
1812     TESTS:
1813
1814     Unembedding is the inverse of embedding::
1815
1816         sage: set_random_seed()
1817         sage: F = QuadraticField(-1, 'i')
1818         sage: M = random_matrix(F, 3)
1819         sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1820         True
1821
1822     """
1823     n = ZZ(M.nrows())
1824     if M.ncols() != n:
1825         raise ValueError("the matrix 'M' must be square")
1826     if not n.mod(2).is_zero():
1827         raise ValueError("the matrix 'M' must be a complex embedding")
1828
1829     F = QuadraticField(-1, 'i')
1830     i = F.gen()
1831
1832     # Go top-left to bottom-right (reading order), converting every
1833     # 2-by-2 block we see to a single complex element.
1834     elements = []
1835     for k in xrange(n/2):
1836         for j in xrange(n/2):
1837             submat = M[2*k:2*k+2,2*j:2*j+2]
1838             if submat[0,0] != submat[1,1]:
1839                 raise ValueError('bad on-diagonal submatrix')
1840             if submat[0,1] != -submat[1,0]:
1841                 raise ValueError('bad off-diagonal submatrix')
1842             z = submat[0,0] + submat[0,1]*i
1843             elements.append(z)
1844
1845     return matrix(F, n/2, elements)
1846
1847
1848 def _embed_quaternion_matrix(M):
1849     """
1850     Embed the n-by-n quaternion matrix ``M`` into the space of real
1851     matrices of size 4n-by-4n by first sending each quaternion entry
1852     `z = a + bi + cj + dk` to the block-complex matrix
1853     ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1854     a real matrix.
1855
1856     EXAMPLES::
1857
1858         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1859         sage: i,j,k = Q.gens()
1860         sage: x = 1 + 2*i + 3*j + 4*k
1861         sage: M = matrix(Q, 1, [[x]])
1862         sage: _embed_quaternion_matrix(M)
1863         [ 1  2  3  4]
1864         [-2  1 -4  3]
1865         [-3  4  1 -2]
1866         [-4 -3  2  1]
1867
1868     Embedding is a homomorphism (isomorphism, in fact)::
1869
1870         sage: set_random_seed()
1871         sage: n = ZZ.random_element(5)
1872         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1873         sage: X = random_matrix(Q, n)
1874         sage: Y = random_matrix(Q, n)
1875         sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1876         sage: expected = _embed_quaternion_matrix(X*Y)
1877         sage: actual == expected
1878         True
1879
1880     """
1881     quaternions = M.base_ring()
1882     n = M.nrows()
1883     if M.ncols() != n:
1884         raise ValueError("the matrix 'M' must be square")
1885
1886     F = QuadraticField(-1, 'i')
1887     i = F.gen()
1888
1889     blocks = []
1890     for z in M.list():
1891         t = z.coefficient_tuple()
1892         a = t[0]
1893         b = t[1]
1894         c = t[2]
1895         d = t[3]
1896         cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1897                                     [-c + d*i, a - b*i]])
1898         blocks.append(_embed_complex_matrix(cplx_matrix))
1899
1900     # We should have real entries by now, so use the realest field
1901     # we've got for the return value.
1902     return block_matrix(quaternions.base_ring(), n, blocks)
1903
1904
1905 def _unembed_quaternion_matrix(M):
1906     """
1907     The inverse of _embed_quaternion_matrix().
1908
1909     EXAMPLES::
1910
1911         sage: M = matrix(QQ, [[ 1,  2,  3,  4],
1912         ....:                 [-2,  1, -4,  3],
1913         ....:                 [-3,  4,  1, -2],
1914         ....:                 [-4, -3,  2,  1]])
1915         sage: _unembed_quaternion_matrix(M)
1916         [1 + 2*i + 3*j + 4*k]
1917
1918     TESTS:
1919
1920     Unembedding is the inverse of embedding::
1921
1922         sage: set_random_seed()
1923         sage: Q = QuaternionAlgebra(QQ, -1, -1)
1924         sage: M = random_matrix(Q, 3)
1925         sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1926         True
1927
1928     """
1929     n = ZZ(M.nrows())
1930     if M.ncols() != n:
1931         raise ValueError("the matrix 'M' must be square")
1932     if not n.mod(4).is_zero():
1933         raise ValueError("the matrix 'M' must be a complex embedding")
1934
1935     Q = QuaternionAlgebra(QQ,-1,-1)
1936     i,j,k = Q.gens()
1937
1938     # Go top-left to bottom-right (reading order), converting every
1939     # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1940     # quaternion block.
1941     elements = []
1942     for l in xrange(n/4):
1943         for m in xrange(n/4):
1944             submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1945             if submat[0,0] != submat[1,1].conjugate():
1946                 raise ValueError('bad on-diagonal submatrix')
1947             if submat[0,1] != -submat[1,0].conjugate():
1948                 raise ValueError('bad off-diagonal submatrix')
1949             z  = submat[0,0].real() + submat[0,0].imag()*i
1950             z += submat[0,1].real()*j + submat[0,1].imag()*k
1951             elements.append(z)
1952
1953     return matrix(Q, n/4, elements)
1954
1955
1956 # The usual inner product on R^n.
1957 def _usual_ip(x,y):
1958     return x.vector().inner_product(y.vector())
1959
1960 # The inner product used for the real symmetric simple EJA.
1961 # We keep it as a separate function because e.g. the complex
1962 # algebra uses the same inner product, except divided by 2.
1963 def _matrix_ip(X,Y):
1964     X_mat = X.natural_representation()
1965     Y_mat = Y.natural_representation()
1966     return (X_mat*Y_mat).trace()
1967
1968
1969 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
1970     """
1971     The rank-n simple EJA consisting of real symmetric n-by-n
1972     matrices, the usual symmetric Jordan product, and the trace inner
1973     product. It has dimension `(n^2 + n)/2` over the reals.
1974
1975     EXAMPLES::
1976
1977         sage: J = RealSymmetricEJA(2)
1978         sage: e0, e1, e2 = J.gens()
1979         sage: e0*e0
1980         e0
1981         sage: e1*e1
1982         e0 + e2
1983         sage: e2*e2
1984         e2
1985
1986     TESTS:
1987
1988     The degree of this algebra is `(n^2 + n) / 2`::
1989
1990         sage: set_random_seed()
1991         sage: n = ZZ.random_element(1,5)
1992         sage: J = RealSymmetricEJA(n)
1993         sage: J.degree() == (n^2 + n)/2
1994         True
1995
1996     The Jordan multiplication is what we think it is::
1997
1998         sage: set_random_seed()
1999         sage: n = ZZ.random_element(1,5)
2000         sage: J = RealSymmetricEJA(n)
2001         sage: x = J.random_element()
2002         sage: y = J.random_element()
2003         sage: actual = (x*y).natural_representation()
2004         sage: X = x.natural_representation()
2005         sage: Y = y.natural_representation()
2006         sage: expected = (X*Y + Y*X)/2
2007         sage: actual == expected
2008         True
2009         sage: J(expected) == x*y
2010         True
2011
2012     """
2013     @staticmethod
2014     def __classcall_private__(cls, n, field=QQ):
2015         S = _real_symmetric_basis(n, field=field)
2016         (Qs, T) = _multiplication_table_from_matrix_basis(S)
2017
2018         fdeja = super(RealSymmetricEJA, cls)
2019         return fdeja.__classcall_private__(cls,
2020                                            field,
2021                                            Qs,
2022                                            rank=n,
2023                                            natural_basis=T)
2024
2025     def inner_product(self, x, y):
2026         return _matrix_ip(x,y)
2027
2028
2029 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
2030     """
2031     The rank-n simple EJA consisting of complex Hermitian n-by-n
2032     matrices over the real numbers, the usual symmetric Jordan product,
2033     and the real-part-of-trace inner product. It has dimension `n^2` over
2034     the reals.
2035
2036     TESTS:
2037
2038     The degree of this algebra is `n^2`::
2039
2040         sage: set_random_seed()
2041         sage: n = ZZ.random_element(1,5)
2042         sage: J = ComplexHermitianEJA(n)
2043         sage: J.degree() == n^2
2044         True
2045
2046     The Jordan multiplication is what we think it is::
2047
2048         sage: set_random_seed()
2049         sage: n = ZZ.random_element(1,5)
2050         sage: J = ComplexHermitianEJA(n)
2051         sage: x = J.random_element()
2052         sage: y = J.random_element()
2053         sage: actual = (x*y).natural_representation()
2054         sage: X = x.natural_representation()
2055         sage: Y = y.natural_representation()
2056         sage: expected = (X*Y + Y*X)/2
2057         sage: actual == expected
2058         True
2059         sage: J(expected) == x*y
2060         True
2061
2062     """
2063     @staticmethod
2064     def __classcall_private__(cls, n, field=QQ):
2065         S = _complex_hermitian_basis(n)
2066         (Qs, T) = _multiplication_table_from_matrix_basis(S)
2067
2068         fdeja = super(ComplexHermitianEJA, cls)
2069         return fdeja.__classcall_private__(cls,
2070                                            field,
2071                                            Qs,
2072                                            rank=n,
2073                                            natural_basis=T)
2074
2075     def inner_product(self, x, y):
2076         # Since a+bi on the diagonal is represented as
2077         #
2078         #   a + bi  = [  a  b  ]
2079         #             [ -b  a  ],
2080         #
2081         # we'll double-count the "a" entries if we take the trace of
2082         # the embedding.
2083         return _matrix_ip(x,y)/2
2084
2085
2086 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
2087     """
2088     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2089     matrices, the usual symmetric Jordan product, and the
2090     real-part-of-trace inner product. It has dimension `2n^2 - n` over
2091     the reals.
2092
2093     TESTS:
2094
2095     The degree of this algebra is `n^2`::
2096
2097         sage: set_random_seed()
2098         sage: n = ZZ.random_element(1,5)
2099         sage: J = QuaternionHermitianEJA(n)
2100         sage: J.degree() == 2*(n^2) - n
2101         True
2102
2103     The Jordan multiplication is what we think it is::
2104
2105         sage: set_random_seed()
2106         sage: n = ZZ.random_element(1,5)
2107         sage: J = QuaternionHermitianEJA(n)
2108         sage: x = J.random_element()
2109         sage: y = J.random_element()
2110         sage: actual = (x*y).natural_representation()
2111         sage: X = x.natural_representation()
2112         sage: Y = y.natural_representation()
2113         sage: expected = (X*Y + Y*X)/2
2114         sage: actual == expected
2115         True
2116         sage: J(expected) == x*y
2117         True
2118
2119     """
2120     @staticmethod
2121     def __classcall_private__(cls, n, field=QQ):
2122         S = _quaternion_hermitian_basis(n)
2123         (Qs, T) = _multiplication_table_from_matrix_basis(S)
2124
2125         fdeja = super(QuaternionHermitianEJA, cls)
2126         return fdeja.__classcall_private__(cls,
2127                                            field,
2128                                            Qs,
2129                                            rank=n,
2130                                            natural_basis=T)
2131
2132     def inner_product(self, x, y):
2133         # Since a+bi+cj+dk on the diagonal is represented as
2134         #
2135         #   a + bi +cj + dk = [  a  b  c  d]
2136         #                     [ -b  a -d  c]
2137         #                     [ -c  d  a -b]
2138         #                     [ -d -c  b  a],
2139         #
2140         # we'll quadruple-count the "a" entries if we take the trace of
2141         # the embedding.
2142         return _matrix_ip(x,y)/4
2143
2144
2145 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
2146     """
2147     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2148     with the usual inner product and jordan product ``x*y =
2149     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2150     the reals.
2151
2152     EXAMPLES:
2153
2154     This multiplication table can be verified by hand::
2155
2156         sage: J = JordanSpinEJA(4)
2157         sage: e0,e1,e2,e3 = J.gens()
2158         sage: e0*e0
2159         e0
2160         sage: e0*e1
2161         e1
2162         sage: e0*e2
2163         e2
2164         sage: e0*e3
2165         e3
2166         sage: e1*e2
2167         0
2168         sage: e1*e3
2169         0
2170         sage: e2*e3
2171         0
2172
2173     """
2174     @staticmethod
2175     def __classcall_private__(cls, n, field=QQ):
2176         Qs = []
2177         id_matrix = identity_matrix(field, n)
2178         for i in xrange(n):
2179             ei = id_matrix.column(i)
2180             Qi = zero_matrix(field, n)
2181             Qi.set_row(0, ei)
2182             Qi.set_column(0, ei)
2183             Qi += diagonal_matrix(n, [ei[0]]*n)
2184             # The addition of the diagonal matrix adds an extra ei[0] in the
2185             # upper-left corner of the matrix.
2186             Qi[0,0] = Qi[0,0] * ~field(2)
2187             Qs.append(Qi)
2188
2189         # The rank of the spin algebra is two, unless we're in a
2190         # one-dimensional ambient space (because the rank is bounded by
2191         # the ambient dimension).
2192         fdeja = super(JordanSpinEJA, cls)
2193         return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2))
2194
2195     def inner_product(self, x, y):
2196         return _usual_ip(x,y)