mjo/eja/euclidean_jordan_algebra.py

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