mjo/eja/euclidean_jordan_algebra.py

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