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