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