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             ALGORITHM:
 662
 663             We appeal to the quadratic representation as in Koecher's
 664             Theorem 12 in Chapter III, Section 5.
 665
 666             EXAMPLES:
 667
 668             The inverse in the spin factor algebra is given in Alizadeh's
 669             Example 11.11::
 670
 671                 sage: set_random_seed()
 672                 sage: n = ZZ.random_element(1,10)
 673                 sage: J = JordanSpinEJA(n)
 674                 sage: x = J.random_element()
 675                 sage: while x.is_zero():
 676                 ....:     x = J.random_element()
 677                 sage: x_vec = x.vector()
 678                 sage: x0 = x_vec[0]
 679                 sage: x_bar = x_vec[1:]
 680                 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
 681                 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
 682                 sage: x_inverse = coeff*inv_vec
 683                 sage: x.inverse() == J(x_inverse)
 684                 True
 685
 686             TESTS:
 687
 688             The identity element is its own inverse::
 689
 690                 sage: set_random_seed()
 691                 sage: J = random_eja()
 692                 sage: J.one().inverse() == J.one()
 693                 True
 694
 695             If an element has an inverse, it acts like one::
 696
 697                 sage: set_random_seed()
 698                 sage: J = random_eja()
 699                 sage: x = J.random_element()
 700                 sage: (not x.is_invertible()) or (x.inverse()*x == J.one())
 701                 True
 702
 703             The inverse of the inverse is what we started with::
 704
 705                 sage: set_random_seed()
 706                 sage: J = random_eja()
 707                 sage: x = J.random_element()
 708                 sage: (not x.is_invertible()) or (x.inverse().inverse() == x)
 709                 True
 710
 711             The zero element is never invertible::
 712
 713                 sage: set_random_seed()
 714                 sage: J = random_eja().zero().inverse()
 715                 Traceback (most recent call last):
 716                 ...
 717                 ValueError: element is not invertible
 718
 719             """
 720             if not self.is_invertible():
 721                 raise ValueError("element is not invertible")
 722
 723             P = self.parent()
 724             return P(self.quadratic_representation().inverse()*self.vector())
 725
 726
 727         def is_invertible(self):
 728             """
 729             Return whether or not this element is invertible.
 730
 731             We can't use the superclass method because it relies on
 732             the algebra being associative.
 733
 734             ALGORITHM:
 735
 736             The usual way to do this is to check if the determinant is
 737             zero, but we need the characteristic polynomial for the
 738             determinant. The minimal polynomial is a lot easier to get,
 739             so we use Corollary 2 in Chapter V of Koecher to check
 740             whether or not the paren't algebra's zero element is a root
 741             of this element's minimal polynomial.
 742
 743             TESTS:
 744
 745             The identity element is always invertible::
 746
 747                 sage: set_random_seed()
 748                 sage: J = random_eja()
 749                 sage: J.one().is_invertible()
 750                 True
 751
 752             The zero element is never invertible::
 753
 754                 sage: set_random_seed()
 755                 sage: J = random_eja()
 756                 sage: J.zero().is_invertible()
 757                 False
 758
 759             """
 760             zero = self.parent().zero()
 761             p = self.minimal_polynomial()
 762             return not (p(zero) == zero)
 763
 764
 765         def is_nilpotent(self):
 766             """
 767             Return whether or not some power of this element is zero.
 768
 769             The superclass method won't work unless we're in an
 770             associative algebra, and we aren't. However, we generate
 771             an assocoative subalgebra and we're nilpotent there if and
 772             only if we're nilpotent here (probably).
 773
 774             TESTS:
 775
 776             The identity element is never nilpotent::
 777
 778                 sage: set_random_seed()
 779                 sage: random_eja().one().is_nilpotent()
 780                 False
 781
 782             The additive identity is always nilpotent::
 783
 784                 sage: set_random_seed()
 785                 sage: random_eja().zero().is_nilpotent()
 786                 True
 787
 788             """
 789             # The element we're going to call "is_nilpotent()" on.
 790             # Either myself, interpreted as an element of a finite-
 791             # dimensional algebra, or an element of an associative
 792             # subalgebra.
 793             elt = None
 794
 795             if self.parent().is_associative():
 796                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 797             else:
 798                 V = self.span_of_powers()
 799                 assoc_subalg = self.subalgebra_generated_by()
 800                 # Mis-design warning: the basis used for span_of_powers()
 801                 # and subalgebra_generated_by() must be the same, and in
 802                 # the same order!
 803                 elt = assoc_subalg(V.coordinates(self.vector()))
 804
 805             # Recursive call, but should work since elt lives in an
 806             # associative algebra.
 807             return elt.is_nilpotent()
 808
 809
 810         def is_regular(self):
 811             """
 812             Return whether or not this is a regular element.
 813
 814             EXAMPLES:
 815
 816             The identity element always has degree one, but any element
 817             linearly-independent from it is regular::
 818
 819                 sage: J = JordanSpinEJA(5)
 820                 sage: J.one().is_regular()
 821                 False
 822                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 823                 sage: for x in J.gens():
 824                 ....:     (J.one() + x).is_regular()
 825                 False
 826                 True
 827                 True
 828                 True
 829                 True
 830
 831             """
 832             return self.degree() == self.parent().rank()
 833
 834
 835         def degree(self):
 836             """
 837             Compute the degree of this element the straightforward way
 838             according to the definition; by appending powers to a list
 839             and figuring out its dimension (that is, whether or not
 840             they're linearly dependent).
 841
 842             EXAMPLES::
 843
 844                 sage: J = JordanSpinEJA(4)
 845                 sage: J.one().degree()
 846                 1
 847                 sage: e0,e1,e2,e3 = J.gens()
 848                 sage: (e0 - e1).degree()
 849                 2
 850
 851             In the spin factor algebra (of rank two), all elements that
 852             aren't multiples of the identity are regular::
 853
 854                 sage: set_random_seed()
 855                 sage: n = ZZ.random_element(1,10)
 856                 sage: J = JordanSpinEJA(n)
 857                 sage: x = J.random_element()
 858                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 859                 True
 860
 861             """
 862             return self.span_of_powers().dimension()
 863
 864
 865         def minimal_polynomial(self):
 866             """
 867             Return the minimal polynomial of this element,
 868             as a function of the variable `t`.
 869
 870             ALGORITHM:
 871
 872             We restrict ourselves to the associative subalgebra
 873             generated by this element, and then return the minimal
 874             polynomial of this element's operator matrix (in that
 875             subalgebra). This works by Baes Proposition 2.3.16.
 876
 877             TESTS:
 878
 879             The minimal polynomial of the identity and zero elements are
 880             always the same::
 881
 882                 sage: set_random_seed()
 883                 sage: J = random_eja()
 884                 sage: J.one().minimal_polynomial()
 885                 t - 1
 886                 sage: J.zero().minimal_polynomial()
 887                 t
 888
 889             The degree of an element is (by one definition) the degree
 890             of its minimal polynomial::
 891
 892                 sage: set_random_seed()
 893                 sage: x = random_eja().random_element()
 894                 sage: x.degree() == x.minimal_polynomial().degree()
 895                 True
 896
 897             The minimal polynomial and the characteristic polynomial coincide
 898             and are known (see Alizadeh, Example 11.11) for all elements of
 899             the spin factor algebra that aren't scalar multiples of the
 900             identity::
 901
 902                 sage: set_random_seed()
 903                 sage: n = ZZ.random_element(2,10)
 904                 sage: J = JordanSpinEJA(n)
 905                 sage: y = J.random_element()
 906                 sage: while y == y.coefficient(0)*J.one():
 907                 ....:     y = J.random_element()
 908                 sage: y0 = y.vector()[0]
 909                 sage: y_bar = y.vector()[1:]
 910                 sage: actual = y.minimal_polynomial()
 911                 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
 912                 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
 913                 sage: bool(actual == expected)
 914                 True
 915
 916             The minimal polynomial should always kill its element::
 917
 918                 sage: set_random_seed()
 919                 sage: x = random_eja().random_element()
 920                 sage: p = x.minimal_polynomial()
 921                 sage: x.apply_univariate_polynomial(p)
 922                 0
 923
 924             """
 925             V = self.span_of_powers()
 926             assoc_subalg = self.subalgebra_generated_by()
 927             # Mis-design warning: the basis used for span_of_powers()
 928             # and subalgebra_generated_by() must be the same, and in
 929             # the same order!
 930             elt = assoc_subalg(V.coordinates(self.vector()))
 931
 932             # We get back a symbolic polynomial in 'x' but want a real
 933             # polynomial in 't'.
 934             p_of_x = elt.operator_matrix().minimal_polynomial()
 935             return p_of_x.change_variable_name('t')
 936
 937
 938         def natural_representation(self):
 939             """
 940             Return a more-natural representation of this element.
 941
 942             Every finite-dimensional Euclidean Jordan Algebra is a
 943             direct sum of five simple algebras, four of which comprise
 944             Hermitian matrices. This method returns the original
 945             "natural" representation of this element as a Hermitian
 946             matrix, if it has one. If not, you get the usual representation.
 947
 948             EXAMPLES::
 949
 950                 sage: J = ComplexHermitianEJA(3)
 951                 sage: J.one()
 952                 e0 + e5 + e8
 953                 sage: J.one().natural_representation()
 954                 [1 0 0 0 0 0]
 955                 [0 1 0 0 0 0]
 956                 [0 0 1 0 0 0]
 957                 [0 0 0 1 0 0]
 958                 [0 0 0 0 1 0]
 959                 [0 0 0 0 0 1]
 960
 961             ::
 962
 963                 sage: J = QuaternionHermitianEJA(3)
 964                 sage: J.one()
 965                 e0 + e9 + e14
 966                 sage: J.one().natural_representation()
 967                 [1 0 0 0 0 0 0 0 0 0 0 0]
 968                 [0 1 0 0 0 0 0 0 0 0 0 0]
 969                 [0 0 1 0 0 0 0 0 0 0 0 0]
 970                 [0 0 0 1 0 0 0 0 0 0 0 0]
 971                 [0 0 0 0 1 0 0 0 0 0 0 0]
 972                 [0 0 0 0 0 1 0 0 0 0 0 0]
 973                 [0 0 0 0 0 0 1 0 0 0 0 0]
 974                 [0 0 0 0 0 0 0 1 0 0 0 0]
 975                 [0 0 0 0 0 0 0 0 1 0 0 0]
 976                 [0 0 0 0 0 0 0 0 0 1 0 0]
 977                 [0 0 0 0 0 0 0 0 0 0 1 0]
 978                 [0 0 0 0 0 0 0 0 0 0 0 1]
 979
 980             """
 981             B = self.parent().natural_basis()
 982             W = B[0].matrix_space()
 983             return W.linear_combination(zip(self.vector(), B))
 984
 985
 986         def operator_matrix(self):
 987             """
 988             Return the matrix that represents left- (or right-)
 989             multiplication by this element in the parent algebra.
 990
 991             We have to override this because the superclass method
 992             returns a matrix that acts on row vectors (that is, on
 993             the right).
 994
 995             EXAMPLES:
 996
 997             Test the first polarization identity from my notes, Koecher Chapter
 998             III, or from Baes (2.3)::
 999
1000                 sage: set_random_seed()
1001                 sage: J = random_eja()
1002                 sage: x = J.random_element()
1003                 sage: y = J.random_element()
1004                 sage: Lx = x.operator_matrix()
1005                 sage: Ly = y.operator_matrix()
1006                 sage: Lxx = (x*x).operator_matrix()
1007                 sage: Lxy = (x*y).operator_matrix()
1008                 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
1009                 True
1010
1011             Test the second polarization identity from my notes or from
1012             Baes (2.4)::
1013
1014                 sage: set_random_seed()
1015                 sage: J = random_eja()
1016                 sage: x = J.random_element()
1017                 sage: y = J.random_element()
1018                 sage: z = J.random_element()
1019                 sage: Lx = x.operator_matrix()
1020                 sage: Ly = y.operator_matrix()
1021                 sage: Lz = z.operator_matrix()
1022                 sage: Lzy = (z*y).operator_matrix()
1023                 sage: Lxy = (x*y).operator_matrix()
1024                 sage: Lxz = (x*z).operator_matrix()
1025                 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
1026                 True
1027
1028             Test the third polarization identity from my notes or from
1029             Baes (2.5)::
1030
1031                 sage: set_random_seed()
1032                 sage: J = random_eja()
1033                 sage: u = J.random_element()
1034                 sage: y = J.random_element()
1035                 sage: z = J.random_element()
1036                 sage: Lu = u.operator_matrix()
1037                 sage: Ly = y.operator_matrix()
1038                 sage: Lz = z.operator_matrix()
1039                 sage: Lzy = (z*y).operator_matrix()
1040                 sage: Luy = (u*y).operator_matrix()
1041                 sage: Luz = (u*z).operator_matrix()
1042                 sage: Luyz = (u*(y*z)).operator_matrix()
1043                 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
1044                 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
1045                 sage: bool(lhs == rhs)
1046                 True
1047
1048             """
1049             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
1050             return fda_elt.matrix().transpose()
1051
1052
1053         def quadratic_representation(self, other=None):
1054             """
1055             Return the quadratic representation of this element.
1056
1057             EXAMPLES:
1058
1059             The explicit form in the spin factor algebra is given by
1060             Alizadeh's Example 11.12::
1061
1062                 sage: set_random_seed()
1063                 sage: n = ZZ.random_element(1,10)
1064                 sage: J = JordanSpinEJA(n)
1065                 sage: x = J.random_element()
1066                 sage: x_vec = x.vector()
1067                 sage: x0 = x_vec[0]
1068                 sage: x_bar = x_vec[1:]
1069                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
1070                 sage: B = 2*x0*x_bar.row()
1071                 sage: C = 2*x0*x_bar.column()
1072                 sage: D = identity_matrix(QQ, n-1)
1073                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
1074                 sage: D = D + 2*x_bar.tensor_product(x_bar)
1075                 sage: Q = block_matrix(2,2,[A,B,C,D])
1076                 sage: Q == x.quadratic_representation()
1077                 True
1078
1079             Test all of the properties from Theorem 11.2 in Alizadeh::
1080
1081                 sage: set_random_seed()
1082                 sage: J = random_eja()
1083                 sage: x = J.random_element()
1084                 sage: y = J.random_element()
1085                 sage: Lx = x.operator_matrix()
1086                 sage: Lxx = (x*x).operator_matrix()
1087                 sage: Qx = x.quadratic_representation()
1088                 sage: Qy = y.quadratic_representation()
1089                 sage: Qxy = x.quadratic_representation(y)
1090                 sage: Qex = J.one().quadratic_representation(x)
1091                 sage: n = ZZ.random_element(10)
1092                 sage: Qxn = (x^n).quadratic_representation()
1093
1094             Property 1:
1095
1096                 sage: 2*Qxy == (x+y).quadratic_representation() - Qx - Qy
1097                 True
1098
1099             Property 2:
1100
1101                 sage: alpha = QQ.random_element()
1102                 sage: (alpha*x).quadratic_representation() == (alpha^2)*Qx
1103                 True
1104
1105             Property 3:
1106
1107                 sage: not x.is_invertible() or (
1108                 ....:     Qx*x.inverse().vector() == x.vector() )
1109                 True
1110
1111                 sage: not x.is_invertible() or (
1112                 ....:   Qx.inverse()
1113                 ....:   ==
1114                 ....:   x.inverse().quadratic_representation() )
1115                 True
1116
1117                 sage: Qxy*(J.one().vector()) == (x*y).vector()
1118                 True
1119
1120             Property 4:
1121
1122                 sage: not x.is_invertible() or (
1123                 ....:   x.quadratic_representation(x.inverse())*Qx
1124                 ....:   == Qx*x.quadratic_representation(x.inverse()) )
1125                 True
1126
1127                 sage: not x.is_invertible() or (
1128                 ....:   x.quadratic_representation(x.inverse())*Qx
1129                 ....:   ==
1130                 ....:   2*x.operator_matrix()*Qex - Qx )
1131                 True
1132
1133                 sage: 2*x.operator_matrix()*Qex - Qx == Lxx
1134                 True
1135
1136             Property 5:
1137
1138                 sage: J(Qy*x.vector()).quadratic_representation() == Qy*Qx*Qy
1139                 True
1140
1141             Property 6:
1142
1143                 sage: Qxn == (Qx)^n
1144                 True
1145
1146             Property 7:
1147
1148                 sage: not x.is_invertible() or (
1149                 ....:   Qx*x.inverse().operator_matrix() == Lx )
1150                 True
1151
1152             Property 8:
1153
1154                 sage: not x.operator_commutes_with(y) or (
1155                 ....:   J(Qx*y.vector())^n == J(Qxn*(y^n).vector()) )
1156                 True
1157
1158             """
1159             if other is None:
1160                 other=self
1161             elif not other in self.parent():
1162                 raise TypeError("'other' must live in the same algebra")
1163
1164             L = self.operator_matrix()
1165             M = other.operator_matrix()
1166             return ( L*M + M*L - (self*other).operator_matrix() )
1167
1168
1169         def span_of_powers(self):
1170             """
1171             Return the vector space spanned by successive powers of
1172             this element.
1173             """
1174             # The dimension of the subalgebra can't be greater than
1175             # the big algebra, so just put everything into a list
1176             # and let span() get rid of the excess.
1177             #
1178             # We do the extra ambient_vector_space() in case we're messing
1179             # with polynomials and the direct parent is a module.
1180             V = self.vector().parent().ambient_vector_space()
1181             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
1182
1183
1184         def subalgebra_generated_by(self):
1185             """
1186             Return the associative subalgebra of the parent EJA generated
1187             by this element.
1188
1189             TESTS::
1190
1191                 sage: set_random_seed()
1192                 sage: x = random_eja().random_element()
1193                 sage: x.subalgebra_generated_by().is_associative()
1194                 True
1195
1196             Squaring in the subalgebra should be the same thing as
1197             squaring in the superalgebra::
1198
1199                 sage: set_random_seed()
1200                 sage: x = random_eja().random_element()
1201                 sage: u = x.subalgebra_generated_by().random_element()
1202                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
1203                 True
1204
1205             """
1206             # First get the subspace spanned by the powers of myself...
1207             V = self.span_of_powers()
1208             F = self.base_ring()
1209
1210             # Now figure out the entries of the right-multiplication
1211             # matrix for the successive basis elements b0, b1,... of
1212             # that subspace.
1213             mats = []
1214             for b_right in V.basis():
1215                 eja_b_right = self.parent()(b_right)
1216                 b_right_rows = []
1217                 # The first row of the right-multiplication matrix by
1218                 # b1 is what we get if we apply that matrix to b1. The
1219                 # second row of the right multiplication matrix by b1
1220                 # is what we get when we apply that matrix to b2...
1221                 #
1222                 # IMPORTANT: this assumes that all vectors are COLUMN
1223                 # vectors, unlike our superclass (which uses row vectors).
1224                 for b_left in V.basis():
1225                     eja_b_left = self.parent()(b_left)
1226                     # Multiply in the original EJA, but then get the
1227                     # coordinates from the subalgebra in terms of its
1228                     # basis.
1229                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
1230                     b_right_rows.append(this_row)
1231                 b_right_matrix = matrix(F, b_right_rows)
1232                 mats.append(b_right_matrix)
1233
1234             # It's an algebra of polynomials in one element, and EJAs
1235             # are power-associative.
1236             #
1237             # TODO: choose generator names intelligently.
1238             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
1239
1240
1241         def subalgebra_idempotent(self):
1242             """
1243             Find an idempotent in the associative subalgebra I generate
1244             using Proposition 2.3.5 in Baes.
1245
1246             TESTS::
1247
1248                 sage: set_random_seed()
1249                 sage: J = random_eja()
1250                 sage: x = J.random_element()
1251                 sage: while x.is_nilpotent():
1252                 ....:     x = J.random_element()
1253                 sage: c = x.subalgebra_idempotent()
1254                 sage: c^2 == c
1255                 True
1256
1257             """
1258             if self.is_nilpotent():
1259                 raise ValueError("this only works with non-nilpotent elements!")
1260
1261             V = self.span_of_powers()
1262             J = self.subalgebra_generated_by()
1263             # Mis-design warning: the basis used for span_of_powers()
1264             # and subalgebra_generated_by() must be the same, and in
1265             # the same order!
1266             u = J(V.coordinates(self.vector()))
1267
1268             # The image of the matrix of left-u^m-multiplication
1269             # will be minimal for some natural number s...
1270             s = 0
1271             minimal_dim = V.dimension()
1272             for i in xrange(1, V.dimension()):
1273                 this_dim = (u**i).operator_matrix().image().dimension()
1274                 if this_dim < minimal_dim:
1275                     minimal_dim = this_dim
1276                     s = i
1277
1278             # Now minimal_matrix should correspond to the smallest
1279             # non-zero subspace in Baes's (or really, Koecher's)
1280             # proposition.
1281             #
1282             # However, we need to restrict the matrix to work on the
1283             # subspace... or do we? Can't we just solve, knowing that
1284             # A(c) = u^(s+1) should have a solution in the big space,
1285             # too?
1286             #
1287             # Beware, solve_right() means that we're using COLUMN vectors.
1288             # Our FiniteDimensionalAlgebraElement superclass uses rows.
1289             u_next = u**(s+1)
1290             A = u_next.operator_matrix()
1291             c_coordinates = A.solve_right(u_next.vector())
1292
1293             # Now c_coordinates is the idempotent we want, but it's in
1294             # the coordinate system of the subalgebra.
1295             #
1296             # We need the basis for J, but as elements of the parent algebra.
1297             #
1298             basis = [self.parent(v) for v in V.basis()]
1299             return self.parent().linear_combination(zip(c_coordinates, basis))
1300
1301
1302         def trace(self):
1303             """
1304             Return my trace, the sum of my eigenvalues.
1305
1306             EXAMPLES::
1307
1308                 sage: J = JordanSpinEJA(3)
1309                 sage: x = sum(J.gens())
1310                 sage: x.trace()
1311                 2
1312
1313             ::
1314
1315                 sage: J = RealCartesianProductEJA(5)
1316                 sage: J.one().trace()
1317                 5
1318
1319             TESTS:
1320
1321             The trace of an element is a real number::
1322
1323                 sage: set_random_seed()
1324                 sage: J = random_eja()
1325                 sage: J.random_element().trace() in J.base_ring()
1326                 True
1327
1328             """
1329             P = self.parent()
1330             r = P.rank()
1331             p = P._charpoly_coeff(r-1)
1332             # The _charpoly_coeff function already adds the factor of
1333             # -1 to ensure that _charpoly_coeff(r-1) is really what
1334             # appears in front of t^{r-1} in the charpoly. However,
1335             # we want the negative of THAT for the trace.
1336             return -p(*self.vector())
1337
1338
1339         def trace_inner_product(self, other):
1340             """
1341             Return the trace inner product of myself and ``other``.
1342
1343             TESTS:
1344
1345             The trace inner product is commutative::
1346
1347                 sage: set_random_seed()
1348                 sage: J = random_eja()
1349                 sage: x = J.random_element(); y = J.random_element()
1350                 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1351                 True
1352
1353             The trace inner product is bilinear::
1354
1355                 sage: set_random_seed()
1356                 sage: J = random_eja()
1357                 sage: x = J.random_element()
1358                 sage: y = J.random_element()
1359                 sage: z = J.random_element()
1360                 sage: a = QQ.random_element();
1361                 sage: actual = (a*(x+z)).trace_inner_product(y)
1362                 sage: expected = ( a*x.trace_inner_product(y) +
1363                 ....:              a*z.trace_inner_product(y) )
1364                 sage: actual == expected
1365                 True
1366                 sage: actual = x.trace_inner_product(a*(y+z))
1367                 sage: expected = ( a*x.trace_inner_product(y) +
1368                 ....:              a*x.trace_inner_product(z) )
1369                 sage: actual == expected
1370                 True
1371
1372             The trace inner product satisfies the compatibility
1373             condition in the definition of a Euclidean Jordan algebra::
1374
1375                 sage: set_random_seed()
1376                 sage: J = random_eja()
1377                 sage: x = J.random_element()
1378                 sage: y = J.random_element()
1379                 sage: z = J.random_element()
1380                 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1381                 True
1382
1383             """
1384             if not other in self.parent():
1385                 raise TypeError("'other' must live in the same algebra")
1386
1387             return (self*other).trace()
1388
1389
1390 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
1391     """
1392     Return the Euclidean Jordan Algebra corresponding to the set
1393     `R^n` under the Hadamard product.
1394
1395     Note: this is nothing more than the Cartesian product of ``n``
1396     copies of the spin algebra. Once Cartesian product algebras
1397     are implemented, this can go.
1398
1399     EXAMPLES:
1400
1401     This multiplication table can be verified by hand::
1402
1403         sage: J = RealCartesianProductEJA(3)
1404         sage: e0,e1,e2 = J.gens()
1405         sage: e0*e0
1406         e0
1407         sage: e0*e1
1408         0
1409         sage: e0*e2
1410         0
1411         sage: e1*e1
1412         e1
1413         sage: e1*e2
1414         0
1415         sage: e2*e2
1416         e2
1417
1418     """
1419     @staticmethod
1420     def __classcall_private__(cls, n, field=QQ):
1421         # The FiniteDimensionalAlgebra constructor takes a list of
1422         # matrices, the ith representing right multiplication by the ith
1423         # basis element in the vector space. So if e_1 = (1,0,0), then
1424         # right (Hadamard) multiplication of x by e_1 picks out the first
1425         # component of x; and likewise for the ith basis element e_i.
1426         Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i))
1427                for i in xrange(n) ]
1428
1429         fdeja = super(RealCartesianProductEJA, cls)
1430         return fdeja.__classcall_private__(cls, field, Qs, rank=n)
1431
1432     def inner_product(self, x, y):
1433         return _usual_ip(x,y)
1434
1435
1436 def random_eja():
1437     """
1438     Return a "random" finite-dimensional Euclidean Jordan Algebra.
1439
1440     ALGORITHM:
1441
1442     For now, we choose a random natural number ``n`` (greater than zero)
1443     and then give you back one of the following:
1444
1445       * The cartesian product of the rational numbers ``n`` times; this is
1446         ``QQ^n`` with the Hadamard product.
1447
1448       * The Jordan spin algebra on ``QQ^n``.
1449
1450       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1451         product.
1452
1453       * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1454         in the space of ``2n``-by-``2n`` real symmetric matrices.
1455
1456       * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1457         in the space of ``4n``-by-``4n`` real symmetric matrices.
1458
1459     Later this might be extended to return Cartesian products of the
1460     EJAs above.
1461
1462     TESTS::
1463
1464         sage: random_eja()
1465         Euclidean Jordan algebra of degree...
1466
1467     """
1468
1469     # The max_n component lets us choose different upper bounds on the
1470     # value "n" that gets passed to the constructor. This is needed
1471     # because e.g. R^{10} is reasonable to test, while the Hermitian
1472     # 10-by-10 quaternion matrices are not.
1473     (constructor, max_n) = choice([(RealCartesianProductEJA, 6),
1474                                    (JordanSpinEJA, 6),
1475                                    (RealSymmetricEJA, 5),
1476                                    (ComplexHermitianEJA, 4),
1477                                    (QuaternionHermitianEJA, 3)])
1478     n = ZZ.random_element(1, max_n)
1479     return constructor(n, field=QQ)
1480
1481
1482
1483 def _real_symmetric_basis(n, field=QQ):
1484     """
1485     Return a basis for the space of real symmetric n-by-n matrices.
1486     """
1487     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1488     # coordinates.
1489     S = []
1490     for i in xrange(n):
1491         for j in xrange(i+1):
1492             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1493             if i == j:
1494                 Sij = Eij
1495             else:
1496                 # Beware, orthogonal but not normalized!
1497                 Sij = Eij + Eij.transpose()
1498             S.append(Sij)
1499     return tuple(S)
1500
1501
1502 def _complex_hermitian_basis(n, field=QQ):
1503     """
1504     Returns a basis for the space of complex Hermitian n-by-n matrices.
1505
1506     TESTS::
1507
1508         sage: set_random_seed()
1509         sage: n = ZZ.random_element(1,5)
1510         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1511         True
1512
1513     """
1514     F = QuadraticField(-1, 'I')
1515     I = F.gen()
1516
1517     # This is like the symmetric case, but we need to be careful:
1518     #
1519     #   * We want conjugate-symmetry, not just symmetry.
1520     #   * The diagonal will (as a result) be real.
1521     #
1522     S = []
1523     for i in xrange(n):
1524         for j in xrange(i+1):
1525             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1526             if i == j:
1527                 Sij = _embed_complex_matrix(Eij)
1528                 S.append(Sij)
1529             else:
1530                 # Beware, orthogonal but not normalized! The second one
1531                 # has a minus because it's conjugated.
1532                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1533                 S.append(Sij_real)
1534                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1535                 S.append(Sij_imag)
1536     return tuple(S)
1537
1538
1539 def _quaternion_hermitian_basis(n, field=QQ):
1540     """
1541     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1542
1543     TESTS::
1544
1545         sage: set_random_seed()
1546         sage: n = ZZ.random_element(1,5)
1547         sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1548         True
1549
1550     """
1551     Q = QuaternionAlgebra(QQ,-1,-1)
1552     I,J,K = Q.gens()
1553
1554     # This is like the symmetric case, but we need to be careful:
1555     #
1556     #   * We want conjugate-symmetry, not just symmetry.
1557     #   * The diagonal will (as a result) be real.
1558     #
1559     S = []
1560     for i in xrange(n):
1561         for j in xrange(i+1):
1562             Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1563             if i == j:
1564                 Sij = _embed_quaternion_matrix(Eij)
1565                 S.append(Sij)
1566             else:
1567                 # Beware, orthogonal but not normalized! The second,
1568                 # third, and fourth ones have a minus because they're
1569                 # conjugated.
1570                 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1571                 S.append(Sij_real)
1572                 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1573                 S.append(Sij_I)
1574                 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1575                 S.append(Sij_J)
1576                 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1577                 S.append(Sij_K)
1578     return tuple(S)
1579
1580
1581 def _mat2vec(m):
1582         return vector(m.base_ring(), m.list())
1583
1584 def _vec2mat(v):
1585         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1586
1587 def _multiplication_table_from_matrix_basis(basis):
1588     """
1589     At least three of the five simple Euclidean Jordan algebras have the
1590     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1591     multiplication on the right is matrix multiplication. Given a basis
1592     for the underlying matrix space, this function returns a
1593     multiplication table (obtained by looping through the basis
1594     elements) for an algebra of those matrices. A reordered copy
1595     of the basis is also returned to work around the fact that
1596     the ``span()`` in this function will change the order of the basis
1597     from what we think it is, to... something else.
1598     """
1599     # In S^2, for example, we nominally have four coordinates even
1600     # though the space is of dimension three only. The vector space V
1601     # is supposed to hold the entire long vector, and the subspace W
1602     # of V will be spanned by the vectors that arise from symmetric
1603     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1604     field = basis[0].base_ring()
1605     dimension = basis[0].nrows()
1606
1607     V = VectorSpace(field, dimension**2)
1608     W = V.span( _mat2vec(s) for s in basis )
1609
1610     # Taking the span above reorders our basis (thanks, jerk!) so we
1611     # need to put our "matrix basis" in the same order as the
1612     # (reordered) vector basis.
1613     S = tuple( _vec2mat(b) for b in W.basis() )
1614
1615     Qs = []
1616     for s in S:
1617         # Brute force the multiplication-by-s matrix by looping
1618         # through all elements of the basis and doing the computation
1619         # to find out what the corresponding row should be. BEWARE:
1620         # these multiplication tables won't be symmetric! It therefore
1621         # becomes REALLY IMPORTANT that the underlying algebra
1622         # constructor uses ROW vectors and not COLUMN vectors. That's
1623         # why we're computing rows here and not columns.
1624         Q_rows = []
1625         for t in S:
1626             this_row = _mat2vec((s*t + t*s)/2)
1627             Q_rows.append(W.coordinates(this_row))
1628         Q = matrix(field, W.dimension(), Q_rows)
1629         Qs.append(Q)
1630
1631     return (Qs, S)
1632
1633
1634 def _embed_complex_matrix(M):
1635     """
1636     Embed the n-by-n complex matrix ``M`` into the space of real
1637     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1638     bi` to the block matrix ``[[a,b],[-b,a]]``.
1639
1640     EXAMPLES::
1641
1642         sage: F = QuadraticField(-1,'i')
1643         sage: x1 = F(4 - 2*i)
1644         sage: x2 = F(1 + 2*i)
1645         sage: x3 = F(-i)
1646         sage: x4 = F(6)
1647         sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1648         sage: _embed_complex_matrix(M)
1649         [ 4 -2| 1  2]
1650         [ 2  4|-2  1]
1651         [-----+-----]
1652         [ 0 -1| 6  0]
1653         [ 1  0| 0  6]
1654
1655     TESTS:
1656
1657     Embedding is a homomorphism (isomorphism, in fact)::
1658
1659         sage: set_random_seed()
1660         sage: n = ZZ.random_element(5)
1661         sage: F = QuadraticField(-1, 'i')
1662         sage: X = random_matrix(F, n)
1663         sage: Y = random_matrix(F, n)
1664         sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1665         sage: expected = _embed_complex_matrix(X*Y)
1666         sage: actual == expected
1667         True
1668
1669     """
1670     n = M.nrows()
1671     if M.ncols() != n:
1672         raise ValueError("the matrix 'M' must be square")
1673     field = M.base_ring()
1674     blocks = []
1675     for z in M.list():
1676         a = z.real()
1677         b = z.imag()
1678         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1679
1680     # We can drop the imaginaries here.
1681     return block_matrix(field.base_ring(), n, blocks)
1682
1683
1684 def _unembed_complex_matrix(M):
1685     """
1686     The inverse of _embed_complex_matrix().
1687
1688     EXAMPLES::
1689
1690         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1691         ....:                 [-2,  1,  -4,  3],
1692         ....:                 [ 9,  10, 11, 12],
1693         ....:                 [-10, 9, -12, 11] ])
1694         sage: _unembed_complex_matrix(A)
1695         [  2*i + 1   4*i + 3]
1696         [ 10*i + 9 12*i + 11]
1697
1698     TESTS:
1699
1700     Unembedding is the inverse of embedding::
1701
1702         sage: set_random_seed()
1703         sage: F = QuadraticField(-1, 'i')
1704         sage: M = random_matrix(F, 3)
1705         sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1706         True
1707
1708     """
1709     n = ZZ(M.nrows())
1710     if M.ncols() != n:
1711         raise ValueError("the matrix 'M' must be square")
1712     if not n.mod(2).is_zero():
1713         raise ValueError("the matrix 'M' must be a complex embedding")
1714
1715     F = QuadraticField(-1, 'i')
1716     i = F.gen()
1717
1718     # Go top-left to bottom-right (reading order), converting every
1719     # 2-by-2 block we see to a single complex element.
1720     elements = []
1721     for k in xrange(n/2):
1722         for j in xrange(n/2):
1723             submat = M[2*k:2*k+2,2*j:2*j+2]
1724             if submat[0,0] != submat[1,1]:
1725                 raise ValueError('bad on-diagonal submatrix')
1726             if submat[0,1] != -submat[1,0]:
1727                 raise ValueError('bad off-diagonal submatrix')
1728             z = submat[0,0] + submat[0,1]*i
1729             elements.append(z)
1730
1731     return matrix(F, n/2, elements)
1732
1733
1734 def _embed_quaternion_matrix(M):
1735     """
1736     Embed the n-by-n quaternion matrix ``M`` into the space of real
1737     matrices of size 4n-by-4n by first sending each quaternion entry
1738     `z = a + bi + cj + dk` to the block-complex matrix
1739     ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1740     a real matrix.
1741
1742     EXAMPLES::
1743
1744         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1745         sage: i,j,k = Q.gens()
1746         sage: x = 1 + 2*i + 3*j + 4*k
1747         sage: M = matrix(Q, 1, [[x]])
1748         sage: _embed_quaternion_matrix(M)
1749         [ 1  2  3  4]
1750         [-2  1 -4  3]
1751         [-3  4  1 -2]
1752         [-4 -3  2  1]
1753
1754     Embedding is a homomorphism (isomorphism, in fact)::
1755
1756         sage: set_random_seed()
1757         sage: n = ZZ.random_element(5)
1758         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1759         sage: X = random_matrix(Q, n)
1760         sage: Y = random_matrix(Q, n)
1761         sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1762         sage: expected = _embed_quaternion_matrix(X*Y)
1763         sage: actual == expected
1764         True
1765
1766     """
1767     quaternions = M.base_ring()
1768     n = M.nrows()
1769     if M.ncols() != n:
1770         raise ValueError("the matrix 'M' must be square")
1771
1772     F = QuadraticField(-1, 'i')
1773     i = F.gen()
1774
1775     blocks = []
1776     for z in M.list():
1777         t = z.coefficient_tuple()
1778         a = t[0]
1779         b = t[1]
1780         c = t[2]
1781         d = t[3]
1782         cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1783                                     [-c + d*i, a - b*i]])
1784         blocks.append(_embed_complex_matrix(cplx_matrix))
1785
1786     # We should have real entries by now, so use the realest field
1787     # we've got for the return value.
1788     return block_matrix(quaternions.base_ring(), n, blocks)
1789
1790
1791 def _unembed_quaternion_matrix(M):
1792     """
1793     The inverse of _embed_quaternion_matrix().
1794
1795     EXAMPLES::
1796
1797         sage: M = matrix(QQ, [[ 1,  2,  3,  4],
1798         ....:                 [-2,  1, -4,  3],
1799         ....:                 [-3,  4,  1, -2],
1800         ....:                 [-4, -3,  2,  1]])
1801         sage: _unembed_quaternion_matrix(M)
1802         [1 + 2*i + 3*j + 4*k]
1803
1804     TESTS:
1805
1806     Unembedding is the inverse of embedding::
1807
1808         sage: set_random_seed()
1809         sage: Q = QuaternionAlgebra(QQ, -1, -1)
1810         sage: M = random_matrix(Q, 3)
1811         sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1812         True
1813
1814     """
1815     n = ZZ(M.nrows())
1816     if M.ncols() != n:
1817         raise ValueError("the matrix 'M' must be square")
1818     if not n.mod(4).is_zero():
1819         raise ValueError("the matrix 'M' must be a complex embedding")
1820
1821     Q = QuaternionAlgebra(QQ,-1,-1)
1822     i,j,k = Q.gens()
1823
1824     # Go top-left to bottom-right (reading order), converting every
1825     # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1826     # quaternion block.
1827     elements = []
1828     for l in xrange(n/4):
1829         for m in xrange(n/4):
1830             submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1831             if submat[0,0] != submat[1,1].conjugate():
1832                 raise ValueError('bad on-diagonal submatrix')
1833             if submat[0,1] != -submat[1,0].conjugate():
1834                 raise ValueError('bad off-diagonal submatrix')
1835             z  = submat[0,0].real() + submat[0,0].imag()*i
1836             z += submat[0,1].real()*j + submat[0,1].imag()*k
1837             elements.append(z)
1838
1839     return matrix(Q, n/4, elements)
1840
1841
1842 # The usual inner product on R^n.
1843 def _usual_ip(x,y):
1844     return x.vector().inner_product(y.vector())
1845
1846 # The inner product used for the real symmetric simple EJA.
1847 # We keep it as a separate function because e.g. the complex
1848 # algebra uses the same inner product, except divided by 2.
1849 def _matrix_ip(X,Y):
1850     X_mat = X.natural_representation()
1851     Y_mat = Y.natural_representation()
1852     return (X_mat*Y_mat).trace()
1853
1854
1855 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
1856     """
1857     The rank-n simple EJA consisting of real symmetric n-by-n
1858     matrices, the usual symmetric Jordan product, and the trace inner
1859     product. It has dimension `(n^2 + n)/2` over the reals.
1860
1861     EXAMPLES::
1862
1863         sage: J = RealSymmetricEJA(2)
1864         sage: e0, e1, e2 = J.gens()
1865         sage: e0*e0
1866         e0
1867         sage: e1*e1
1868         e0 + e2
1869         sage: e2*e2
1870         e2
1871
1872     TESTS:
1873
1874     The degree of this algebra is `(n^2 + n) / 2`::
1875
1876         sage: set_random_seed()
1877         sage: n = ZZ.random_element(1,5)
1878         sage: J = RealSymmetricEJA(n)
1879         sage: J.degree() == (n^2 + n)/2
1880         True
1881
1882     The Jordan multiplication is what we think it is::
1883
1884         sage: set_random_seed()
1885         sage: n = ZZ.random_element(1,5)
1886         sage: J = RealSymmetricEJA(n)
1887         sage: x = J.random_element()
1888         sage: y = J.random_element()
1889         sage: actual = (x*y).natural_representation()
1890         sage: X = x.natural_representation()
1891         sage: Y = y.natural_representation()
1892         sage: expected = (X*Y + Y*X)/2
1893         sage: actual == expected
1894         True
1895         sage: J(expected) == x*y
1896         True
1897
1898     """
1899     @staticmethod
1900     def __classcall_private__(cls, n, field=QQ):
1901         S = _real_symmetric_basis(n, field=field)
1902         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1903
1904         fdeja = super(RealSymmetricEJA, cls)
1905         return fdeja.__classcall_private__(cls,
1906                                            field,
1907                                            Qs,
1908                                            rank=n,
1909                                            natural_basis=T)
1910
1911     def inner_product(self, x, y):
1912         return _matrix_ip(x,y)
1913
1914
1915 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1916     """
1917     The rank-n simple EJA consisting of complex Hermitian n-by-n
1918     matrices over the real numbers, the usual symmetric Jordan product,
1919     and the real-part-of-trace inner product. It has dimension `n^2` over
1920     the reals.
1921
1922     TESTS:
1923
1924     The degree of this algebra is `n^2`::
1925
1926         sage: set_random_seed()
1927         sage: n = ZZ.random_element(1,5)
1928         sage: J = ComplexHermitianEJA(n)
1929         sage: J.degree() == n^2
1930         True
1931
1932     The Jordan multiplication is what we think it is::
1933
1934         sage: set_random_seed()
1935         sage: n = ZZ.random_element(1,5)
1936         sage: J = ComplexHermitianEJA(n)
1937         sage: x = J.random_element()
1938         sage: y = J.random_element()
1939         sage: actual = (x*y).natural_representation()
1940         sage: X = x.natural_representation()
1941         sage: Y = y.natural_representation()
1942         sage: expected = (X*Y + Y*X)/2
1943         sage: actual == expected
1944         True
1945         sage: J(expected) == x*y
1946         True
1947
1948     """
1949     @staticmethod
1950     def __classcall_private__(cls, n, field=QQ):
1951         S = _complex_hermitian_basis(n)
1952         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1953
1954         fdeja = super(ComplexHermitianEJA, cls)
1955         return fdeja.__classcall_private__(cls,
1956                                            field,
1957                                            Qs,
1958                                            rank=n,
1959                                            natural_basis=T)
1960
1961     def inner_product(self, x, y):
1962         # Since a+bi on the diagonal is represented as
1963         #
1964         #   a + bi  = [  a  b  ]
1965         #             [ -b  a  ],
1966         #
1967         # we'll double-count the "a" entries if we take the trace of
1968         # the embedding.
1969         return _matrix_ip(x,y)/2
1970
1971
1972 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1973     """
1974     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1975     matrices, the usual symmetric Jordan product, and the
1976     real-part-of-trace inner product. It has dimension `2n^2 - n` over
1977     the reals.
1978
1979     TESTS:
1980
1981     The degree of this algebra is `n^2`::
1982
1983         sage: set_random_seed()
1984         sage: n = ZZ.random_element(1,5)
1985         sage: J = QuaternionHermitianEJA(n)
1986         sage: J.degree() == 2*(n^2) - n
1987         True
1988
1989     The Jordan multiplication is what we think it is::
1990
1991         sage: set_random_seed()
1992         sage: n = ZZ.random_element(1,5)
1993         sage: J = QuaternionHermitianEJA(n)
1994         sage: x = J.random_element()
1995         sage: y = J.random_element()
1996         sage: actual = (x*y).natural_representation()
1997         sage: X = x.natural_representation()
1998         sage: Y = y.natural_representation()
1999         sage: expected = (X*Y + Y*X)/2
2000         sage: actual == expected
2001         True
2002         sage: J(expected) == x*y
2003         True
2004
2005     """
2006     @staticmethod
2007     def __classcall_private__(cls, n, field=QQ):
2008         S = _quaternion_hermitian_basis(n)
2009         (Qs, T) = _multiplication_table_from_matrix_basis(S)
2010
2011         fdeja = super(QuaternionHermitianEJA, cls)
2012         return fdeja.__classcall_private__(cls,
2013                                            field,
2014                                            Qs,
2015                                            rank=n,
2016                                            natural_basis=T)
2017
2018     def inner_product(self, x, y):
2019         # Since a+bi+cj+dk on the diagonal is represented as
2020         #
2021         #   a + bi +cj + dk = [  a  b  c  d]
2022         #                     [ -b  a -d  c]
2023         #                     [ -c  d  a -b]
2024         #                     [ -d -c  b  a],
2025         #
2026         # we'll quadruple-count the "a" entries if we take the trace of
2027         # the embedding.
2028         return _matrix_ip(x,y)/4
2029
2030
2031 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
2032     """
2033     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2034     with the usual inner product and jordan product ``x*y =
2035     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2036     the reals.
2037
2038     EXAMPLES:
2039
2040     This multiplication table can be verified by hand::
2041
2042         sage: J = JordanSpinEJA(4)
2043         sage: e0,e1,e2,e3 = J.gens()
2044         sage: e0*e0
2045         e0
2046         sage: e0*e1
2047         e1
2048         sage: e0*e2
2049         e2
2050         sage: e0*e3
2051         e3
2052         sage: e1*e2
2053         0
2054         sage: e1*e3
2055         0
2056         sage: e2*e3
2057         0
2058
2059     """
2060     @staticmethod
2061     def __classcall_private__(cls, n, field=QQ):
2062         Qs = []
2063         id_matrix = identity_matrix(field, n)
2064         for i in xrange(n):
2065             ei = id_matrix.column(i)
2066             Qi = zero_matrix(field, n)
2067             Qi.set_row(0, ei)
2068             Qi.set_column(0, ei)
2069             Qi += diagonal_matrix(n, [ei[0]]*n)
2070             # The addition of the diagonal matrix adds an extra ei[0] in the
2071             # upper-left corner of the matrix.
2072             Qi[0,0] = Qi[0,0] * ~field(2)
2073             Qs.append(Qi)
2074
2075         # The rank of the spin algebra is two, unless we're in a
2076         # one-dimensional ambient space (because the rank is bounded by
2077         # the ambient dimension).
2078         fdeja = super(JordanSpinEJA, cls)
2079         return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2))
2080
2081     def inner_product(self, x, y):
2082         return _usual_ip(x,y)