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