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