mjo/eja/euclidean_jordan_algebra.py

   1 """
   2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
   3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
   4 are used in optimization, and have some additional nice methods beyond
   5 what can be supported in a general Jordan Algebra.
   6 """
   7
   8 from sage.categories.magmatic_algebras import MagmaticAlgebras
   9 from sage.structure.element import is_Matrix
  10 from sage.structure.category_object import normalize_names
  11
  12 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra
  13 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element import FiniteDimensionalAlgebraElement
  14 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 = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
 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()
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_matrix()
1381                 sage: Lxx = (x*x).operator_matrix()
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 (
1403                 ....:     Qx*x.inverse().vector() == x.vector() )
1404                 True
1405
1406                 sage: not x.is_invertible() or (
1407                 ....:   Qx.inverse()
1408                 ....:   ==
1409                 ....:   x.inverse().quadratic_representation() )
1410                 True
1411
1412                 sage: Qxy*(J.one().vector()) == (x*y).vector()
1413                 True
1414
1415             Property 4:
1416
1417                 sage: not x.is_invertible() or (
1418                 ....:   x.quadratic_representation(x.inverse())*Qx
1419                 ....:   == Qx*x.quadratic_representation(x.inverse()) )
1420                 True
1421
1422                 sage: not x.is_invertible() or (
1423                 ....:   x.quadratic_representation(x.inverse())*Qx
1424                 ....:   ==
1425                 ....:   2*x.operator_matrix()*Qex - Qx )
1426                 True
1427
1428                 sage: 2*x.operator_matrix()*Qex - Qx == Lxx
1429                 True
1430
1431             Property 5:
1432
1433                 sage: J(Qy*x.vector()).quadratic_representation() == Qy*Qx*Qy
1434                 True
1435
1436             Property 6:
1437
1438                 sage: Qxn == (Qx)^n
1439                 True
1440
1441             Property 7:
1442
1443                 sage: not x.is_invertible() or (
1444                 ....:   Qx*x.inverse().operator_matrix() == Lx )
1445                 True
1446
1447             Property 8:
1448
1449                 sage: not x.operator_commutes_with(y) or (
1450                 ....:   J(Qx*y.vector())^n == J(Qxn*(y^n).vector()) )
1451                 True
1452
1453             """
1454             if other is None:
1455                 other=self
1456             elif not other in self.parent():
1457                 raise TypeError("'other' must live in the same algebra")
1458
1459             L = self.operator_matrix()
1460             M = other.operator_matrix()
1461             return ( L*M + M*L - (self*other).operator_matrix() )
1462
1463
1464         def span_of_powers(self):
1465             """
1466             Return the vector space spanned by successive powers of
1467             this element.
1468             """
1469             # The dimension of the subalgebra can't be greater than
1470             # the big algebra, so just put everything into a list
1471             # and let span() get rid of the excess.
1472             #
1473             # We do the extra ambient_vector_space() in case we're messing
1474             # with polynomials and the direct parent is a module.
1475             V = self.vector().parent().ambient_vector_space()
1476             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
1477
1478
1479         def subalgebra_generated_by(self):
1480             """
1481             Return the associative subalgebra of the parent EJA generated
1482             by this element.
1483
1484             TESTS::
1485
1486                 sage: set_random_seed()
1487                 sage: x = random_eja().random_element()
1488                 sage: x.subalgebra_generated_by().is_associative()
1489                 True
1490
1491             Squaring in the subalgebra should be the same thing as
1492             squaring in the superalgebra::
1493
1494                 sage: set_random_seed()
1495                 sage: x = random_eja().random_element()
1496                 sage: u = x.subalgebra_generated_by().random_element()
1497                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
1498                 True
1499
1500             """
1501             # First get the subspace spanned by the powers of myself...
1502             V = self.span_of_powers()
1503             F = self.base_ring()
1504
1505             # Now figure out the entries of the right-multiplication
1506             # matrix for the successive basis elements b0, b1,... of
1507             # that subspace.
1508             mats = []
1509             for b_right in V.basis():
1510                 eja_b_right = self.parent()(b_right)
1511                 b_right_rows = []
1512                 # The first row of the right-multiplication matrix by
1513                 # b1 is what we get if we apply that matrix to b1. The
1514                 # second row of the right multiplication matrix by b1
1515                 # is what we get when we apply that matrix to b2...
1516                 #
1517                 # IMPORTANT: this assumes that all vectors are COLUMN
1518                 # vectors, unlike our superclass (which uses row vectors).
1519                 for b_left in V.basis():
1520                     eja_b_left = self.parent()(b_left)
1521                     # Multiply in the original EJA, but then get the
1522                     # coordinates from the subalgebra in terms of its
1523                     # basis.
1524                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
1525                     b_right_rows.append(this_row)
1526                 b_right_matrix = matrix(F, b_right_rows)
1527                 mats.append(b_right_matrix)
1528
1529             # It's an algebra of polynomials in one element, and EJAs
1530             # are power-associative.
1531             #
1532             # TODO: choose generator names intelligently.
1533             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
1534
1535
1536         def subalgebra_idempotent(self):
1537             """
1538             Find an idempotent in the associative subalgebra I generate
1539             using Proposition 2.3.5 in Baes.
1540
1541             TESTS::
1542
1543                 sage: set_random_seed()
1544                 sage: J = random_eja()
1545                 sage: x = J.random_element()
1546                 sage: while x.is_nilpotent():
1547                 ....:     x = J.random_element()
1548                 sage: c = x.subalgebra_idempotent()
1549                 sage: c^2 == c
1550                 True
1551
1552             """
1553             if self.is_nilpotent():
1554                 raise ValueError("this only works with non-nilpotent elements!")
1555
1556             V = self.span_of_powers()
1557             J = self.subalgebra_generated_by()
1558             # Mis-design warning: the basis used for span_of_powers()
1559             # and subalgebra_generated_by() must be the same, and in
1560             # the same order!
1561             u = J(V.coordinates(self.vector()))
1562
1563             # The image of the matrix of left-u^m-multiplication
1564             # will be minimal for some natural number s...
1565             s = 0
1566             minimal_dim = V.dimension()
1567             for i in xrange(1, V.dimension()):
1568                 this_dim = (u**i).operator_matrix().image().dimension()
1569                 if this_dim < minimal_dim:
1570                     minimal_dim = this_dim
1571                     s = i
1572
1573             # Now minimal_matrix should correspond to the smallest
1574             # non-zero subspace in Baes's (or really, Koecher's)
1575             # proposition.
1576             #
1577             # However, we need to restrict the matrix to work on the
1578             # subspace... or do we? Can't we just solve, knowing that
1579             # A(c) = u^(s+1) should have a solution in the big space,
1580             # too?
1581             #
1582             # Beware, solve_right() means that we're using COLUMN vectors.
1583             # Our FiniteDimensionalAlgebraElement superclass uses rows.
1584             u_next = u**(s+1)
1585             A = u_next.operator_matrix()
1586             c_coordinates = A.solve_right(u_next.vector())
1587
1588             # Now c_coordinates is the idempotent we want, but it's in
1589             # the coordinate system of the subalgebra.
1590             #
1591             # We need the basis for J, but as elements of the parent algebra.
1592             #
1593             basis = [self.parent(v) for v in V.basis()]
1594             return self.parent().linear_combination(zip(c_coordinates, basis))
1595
1596
1597         def trace(self):
1598             """
1599             Return my trace, the sum of my eigenvalues.
1600
1601             EXAMPLES::
1602
1603                 sage: J = JordanSpinEJA(3)
1604                 sage: x = sum(J.gens())
1605                 sage: x.trace()
1606                 2
1607
1608             ::
1609
1610                 sage: J = RealCartesianProductEJA(5)
1611                 sage: J.one().trace()
1612                 5
1613
1614             TESTS:
1615
1616             The trace of an element is a real number::
1617
1618                 sage: set_random_seed()
1619                 sage: J = random_eja()
1620                 sage: J.random_element().trace() in J.base_ring()
1621                 True
1622
1623             """
1624             P = self.parent()
1625             r = P.rank()
1626             p = P._charpoly_coeff(r-1)
1627             # The _charpoly_coeff function already adds the factor of
1628             # -1 to ensure that _charpoly_coeff(r-1) is really what
1629             # appears in front of t^{r-1} in the charpoly. However,
1630             # we want the negative of THAT for the trace.
1631             return -p(*self.vector())
1632
1633
1634         def trace_inner_product(self, other):
1635             """
1636             Return the trace inner product of myself and ``other``.
1637
1638             TESTS:
1639
1640             The trace inner product is commutative::
1641
1642                 sage: set_random_seed()
1643                 sage: J = random_eja()
1644                 sage: x = J.random_element(); y = J.random_element()
1645                 sage: x.trace_inner_product(y) == y.trace_inner_product(x)
1646                 True
1647
1648             The trace inner product is bilinear::
1649
1650                 sage: set_random_seed()
1651                 sage: J = random_eja()
1652                 sage: x = J.random_element()
1653                 sage: y = J.random_element()
1654                 sage: z = J.random_element()
1655                 sage: a = QQ.random_element();
1656                 sage: actual = (a*(x+z)).trace_inner_product(y)
1657                 sage: expected = ( a*x.trace_inner_product(y) +
1658                 ....:              a*z.trace_inner_product(y) )
1659                 sage: actual == expected
1660                 True
1661                 sage: actual = x.trace_inner_product(a*(y+z))
1662                 sage: expected = ( a*x.trace_inner_product(y) +
1663                 ....:              a*x.trace_inner_product(z) )
1664                 sage: actual == expected
1665                 True
1666
1667             The trace inner product satisfies the compatibility
1668             condition in the definition of a Euclidean Jordan algebra::
1669
1670                 sage: set_random_seed()
1671                 sage: J = random_eja()
1672                 sage: x = J.random_element()
1673                 sage: y = J.random_element()
1674                 sage: z = J.random_element()
1675                 sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
1676                 True
1677
1678             """
1679             if not other in self.parent():
1680                 raise TypeError("'other' must live in the same algebra")
1681
1682             return (self*other).trace()
1683
1684
1685 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
1686     """
1687     Return the Euclidean Jordan Algebra corresponding to the set
1688     `R^n` under the Hadamard product.
1689
1690     Note: this is nothing more than the Cartesian product of ``n``
1691     copies of the spin algebra. Once Cartesian product algebras
1692     are implemented, this can go.
1693
1694     EXAMPLES:
1695
1696     This multiplication table can be verified by hand::
1697
1698         sage: J = RealCartesianProductEJA(3)
1699         sage: e0,e1,e2 = J.gens()
1700         sage: e0*e0
1701         e0
1702         sage: e0*e1
1703         0
1704         sage: e0*e2
1705         0
1706         sage: e1*e1
1707         e1
1708         sage: e1*e2
1709         0
1710         sage: e2*e2
1711         e2
1712
1713     """
1714     @staticmethod
1715     def __classcall_private__(cls, n, field=QQ):
1716         # The FiniteDimensionalAlgebra constructor takes a list of
1717         # matrices, the ith representing right multiplication by the ith
1718         # basis element in the vector space. So if e_1 = (1,0,0), then
1719         # right (Hadamard) multiplication of x by e_1 picks out the first
1720         # component of x; and likewise for the ith basis element e_i.
1721         Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i))
1722                for i in xrange(n) ]
1723
1724         fdeja = super(RealCartesianProductEJA, cls)
1725         return fdeja.__classcall_private__(cls, field, Qs, rank=n)
1726
1727     def inner_product(self, x, y):
1728         return _usual_ip(x,y)
1729
1730
1731 def random_eja():
1732     """
1733     Return a "random" finite-dimensional Euclidean Jordan Algebra.
1734
1735     ALGORITHM:
1736
1737     For now, we choose a random natural number ``n`` (greater than zero)
1738     and then give you back one of the following:
1739
1740       * The cartesian product of the rational numbers ``n`` times; this is
1741         ``QQ^n`` with the Hadamard product.
1742
1743       * The Jordan spin algebra on ``QQ^n``.
1744
1745       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1746         product.
1747
1748       * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1749         in the space of ``2n``-by-``2n`` real symmetric matrices.
1750
1751       * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1752         in the space of ``4n``-by-``4n`` real symmetric matrices.
1753
1754     Later this might be extended to return Cartesian products of the
1755     EJAs above.
1756
1757     TESTS::
1758
1759         sage: random_eja()
1760         Euclidean Jordan algebra of degree...
1761
1762     """
1763
1764     # The max_n component lets us choose different upper bounds on the
1765     # value "n" that gets passed to the constructor. This is needed
1766     # because e.g. R^{10} is reasonable to test, while the Hermitian
1767     # 10-by-10 quaternion matrices are not.
1768     (constructor, max_n) = choice([(RealCartesianProductEJA, 6),
1769                                    (JordanSpinEJA, 6),
1770                                    (RealSymmetricEJA, 5),
1771                                    (ComplexHermitianEJA, 4),
1772                                    (QuaternionHermitianEJA, 3)])
1773     n = ZZ.random_element(1, max_n)
1774     return constructor(n, field=QQ)
1775
1776
1777
1778 def _real_symmetric_basis(n, field=QQ):
1779     """
1780     Return a basis for the space of real symmetric n-by-n matrices.
1781     """
1782     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1783     # coordinates.
1784     S = []
1785     for i in xrange(n):
1786         for j in xrange(i+1):
1787             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1788             if i == j:
1789                 Sij = Eij
1790             else:
1791                 # Beware, orthogonal but not normalized!
1792                 Sij = Eij + Eij.transpose()
1793             S.append(Sij)
1794     return tuple(S)
1795
1796
1797 def _complex_hermitian_basis(n, field=QQ):
1798     """
1799     Returns a basis for the space of complex Hermitian n-by-n matrices.
1800
1801     TESTS::
1802
1803         sage: set_random_seed()
1804         sage: n = ZZ.random_element(1,5)
1805         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1806         True
1807
1808     """
1809     F = QuadraticField(-1, 'I')
1810     I = F.gen()
1811
1812     # This is like the symmetric case, but we need to be careful:
1813     #
1814     #   * We want conjugate-symmetry, not just symmetry.
1815     #   * The diagonal will (as a result) be real.
1816     #
1817     S = []
1818     for i in xrange(n):
1819         for j in xrange(i+1):
1820             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1821             if i == j:
1822                 Sij = _embed_complex_matrix(Eij)
1823                 S.append(Sij)
1824             else:
1825                 # Beware, orthogonal but not normalized! The second one
1826                 # has a minus because it's conjugated.
1827                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1828                 S.append(Sij_real)
1829                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1830                 S.append(Sij_imag)
1831     return tuple(S)
1832
1833
1834 def _quaternion_hermitian_basis(n, field=QQ):
1835     """
1836     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1837
1838     TESTS::
1839
1840         sage: set_random_seed()
1841         sage: n = ZZ.random_element(1,5)
1842         sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1843         True
1844
1845     """
1846     Q = QuaternionAlgebra(QQ,-1,-1)
1847     I,J,K = Q.gens()
1848
1849     # This is like the symmetric case, but we need to be careful:
1850     #
1851     #   * We want conjugate-symmetry, not just symmetry.
1852     #   * The diagonal will (as a result) be real.
1853     #
1854     S = []
1855     for i in xrange(n):
1856         for j in xrange(i+1):
1857             Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1858             if i == j:
1859                 Sij = _embed_quaternion_matrix(Eij)
1860                 S.append(Sij)
1861             else:
1862                 # Beware, orthogonal but not normalized! The second,
1863                 # third, and fourth ones have a minus because they're
1864                 # conjugated.
1865                 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1866                 S.append(Sij_real)
1867                 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1868                 S.append(Sij_I)
1869                 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1870                 S.append(Sij_J)
1871                 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1872                 S.append(Sij_K)
1873     return tuple(S)
1874
1875
1876 def _mat2vec(m):
1877         return vector(m.base_ring(), m.list())
1878
1879 def _vec2mat(v):
1880         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1881
1882 def _multiplication_table_from_matrix_basis(basis):
1883     """
1884     At least three of the five simple Euclidean Jordan algebras have the
1885     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1886     multiplication on the right is matrix multiplication. Given a basis
1887     for the underlying matrix space, this function returns a
1888     multiplication table (obtained by looping through the basis
1889     elements) for an algebra of those matrices. A reordered copy
1890     of the basis is also returned to work around the fact that
1891     the ``span()`` in this function will change the order of the basis
1892     from what we think it is, to... something else.
1893     """
1894     # In S^2, for example, we nominally have four coordinates even
1895     # though the space is of dimension three only. The vector space V
1896     # is supposed to hold the entire long vector, and the subspace W
1897     # of V will be spanned by the vectors that arise from symmetric
1898     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1899     field = basis[0].base_ring()
1900     dimension = basis[0].nrows()
1901
1902     V = VectorSpace(field, dimension**2)
1903     W = V.span( _mat2vec(s) for s in basis )
1904
1905     # Taking the span above reorders our basis (thanks, jerk!) so we
1906     # need to put our "matrix basis" in the same order as the
1907     # (reordered) vector basis.
1908     S = tuple( _vec2mat(b) for b in W.basis() )
1909
1910     Qs = []
1911     for s in S:
1912         # Brute force the multiplication-by-s matrix by looping
1913         # through all elements of the basis and doing the computation
1914         # to find out what the corresponding row should be. BEWARE:
1915         # these multiplication tables won't be symmetric! It therefore
1916         # becomes REALLY IMPORTANT that the underlying algebra
1917         # constructor uses ROW vectors and not COLUMN vectors. That's
1918         # why we're computing rows here and not columns.
1919         Q_rows = []
1920         for t in S:
1921             this_row = _mat2vec((s*t + t*s)/2)
1922             Q_rows.append(W.coordinates(this_row))
1923         Q = matrix(field, W.dimension(), Q_rows)
1924         Qs.append(Q)
1925
1926     return (Qs, S)
1927
1928
1929 def _embed_complex_matrix(M):
1930     """
1931     Embed the n-by-n complex matrix ``M`` into the space of real
1932     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1933     bi` to the block matrix ``[[a,b],[-b,a]]``.
1934
1935     EXAMPLES::
1936
1937         sage: F = QuadraticField(-1,'i')
1938         sage: x1 = F(4 - 2*i)
1939         sage: x2 = F(1 + 2*i)
1940         sage: x3 = F(-i)
1941         sage: x4 = F(6)
1942         sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1943         sage: _embed_complex_matrix(M)
1944         [ 4 -2| 1  2]
1945         [ 2  4|-2  1]
1946         [-----+-----]
1947         [ 0 -1| 6  0]
1948         [ 1  0| 0  6]
1949
1950     TESTS:
1951
1952     Embedding is a homomorphism (isomorphism, in fact)::
1953
1954         sage: set_random_seed()
1955         sage: n = ZZ.random_element(5)
1956         sage: F = QuadraticField(-1, 'i')
1957         sage: X = random_matrix(F, n)
1958         sage: Y = random_matrix(F, n)
1959         sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1960         sage: expected = _embed_complex_matrix(X*Y)
1961         sage: actual == expected
1962         True
1963
1964     """
1965     n = M.nrows()
1966     if M.ncols() != n:
1967         raise ValueError("the matrix 'M' must be square")
1968     field = M.base_ring()
1969     blocks = []
1970     for z in M.list():
1971         a = z.real()
1972         b = z.imag()
1973         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1974
1975     # We can drop the imaginaries here.
1976     return block_matrix(field.base_ring(), n, blocks)
1977
1978
1979 def _unembed_complex_matrix(M):
1980     """
1981     The inverse of _embed_complex_matrix().
1982
1983     EXAMPLES::
1984
1985         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1986         ....:                 [-2,  1,  -4,  3],
1987         ....:                 [ 9,  10, 11, 12],
1988         ....:                 [-10, 9, -12, 11] ])
1989         sage: _unembed_complex_matrix(A)
1990         [  2*i + 1   4*i + 3]
1991         [ 10*i + 9 12*i + 11]
1992
1993     TESTS:
1994
1995     Unembedding is the inverse of embedding::
1996
1997         sage: set_random_seed()
1998         sage: F = QuadraticField(-1, 'i')
1999         sage: M = random_matrix(F, 3)
2000         sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
2001         True
2002
2003     """
2004     n = ZZ(M.nrows())
2005     if M.ncols() != n:
2006         raise ValueError("the matrix 'M' must be square")
2007     if not n.mod(2).is_zero():
2008         raise ValueError("the matrix 'M' must be a complex embedding")
2009
2010     F = QuadraticField(-1, 'i')
2011     i = F.gen()
2012
2013     # Go top-left to bottom-right (reading order), converting every
2014     # 2-by-2 block we see to a single complex element.
2015     elements = []
2016     for k in xrange(n/2):
2017         for j in xrange(n/2):
2018             submat = M[2*k:2*k+2,2*j:2*j+2]
2019             if submat[0,0] != submat[1,1]:
2020                 raise ValueError('bad on-diagonal submatrix')
2021             if submat[0,1] != -submat[1,0]:
2022                 raise ValueError('bad off-diagonal submatrix')
2023             z = submat[0,0] + submat[0,1]*i
2024             elements.append(z)
2025
2026     return matrix(F, n/2, elements)
2027
2028
2029 def _embed_quaternion_matrix(M):
2030     """
2031     Embed the n-by-n quaternion matrix ``M`` into the space of real
2032     matrices of size 4n-by-4n by first sending each quaternion entry
2033     `z = a + bi + cj + dk` to the block-complex matrix
2034     ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
2035     a real matrix.
2036
2037     EXAMPLES::
2038
2039         sage: Q = QuaternionAlgebra(QQ,-1,-1)
2040         sage: i,j,k = Q.gens()
2041         sage: x = 1 + 2*i + 3*j + 4*k
2042         sage: M = matrix(Q, 1, [[x]])
2043         sage: _embed_quaternion_matrix(M)
2044         [ 1  2  3  4]
2045         [-2  1 -4  3]
2046         [-3  4  1 -2]
2047         [-4 -3  2  1]
2048
2049     Embedding is a homomorphism (isomorphism, in fact)::
2050
2051         sage: set_random_seed()
2052         sage: n = ZZ.random_element(5)
2053         sage: Q = QuaternionAlgebra(QQ,-1,-1)
2054         sage: X = random_matrix(Q, n)
2055         sage: Y = random_matrix(Q, n)
2056         sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
2057         sage: expected = _embed_quaternion_matrix(X*Y)
2058         sage: actual == expected
2059         True
2060
2061     """
2062     quaternions = M.base_ring()
2063     n = M.nrows()
2064     if M.ncols() != n:
2065         raise ValueError("the matrix 'M' must be square")
2066
2067     F = QuadraticField(-1, 'i')
2068     i = F.gen()
2069
2070     blocks = []
2071     for z in M.list():
2072         t = z.coefficient_tuple()
2073         a = t[0]
2074         b = t[1]
2075         c = t[2]
2076         d = t[3]
2077         cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
2078                                     [-c + d*i, a - b*i]])
2079         blocks.append(_embed_complex_matrix(cplx_matrix))
2080
2081     # We should have real entries by now, so use the realest field
2082     # we've got for the return value.
2083     return block_matrix(quaternions.base_ring(), n, blocks)
2084
2085
2086 def _unembed_quaternion_matrix(M):
2087     """
2088     The inverse of _embed_quaternion_matrix().
2089
2090     EXAMPLES::
2091
2092         sage: M = matrix(QQ, [[ 1,  2,  3,  4],
2093         ....:                 [-2,  1, -4,  3],
2094         ....:                 [-3,  4,  1, -2],
2095         ....:                 [-4, -3,  2,  1]])
2096         sage: _unembed_quaternion_matrix(M)
2097         [1 + 2*i + 3*j + 4*k]
2098
2099     TESTS:
2100
2101     Unembedding is the inverse of embedding::
2102
2103         sage: set_random_seed()
2104         sage: Q = QuaternionAlgebra(QQ, -1, -1)
2105         sage: M = random_matrix(Q, 3)
2106         sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
2107         True
2108
2109     """
2110     n = ZZ(M.nrows())
2111     if M.ncols() != n:
2112         raise ValueError("the matrix 'M' must be square")
2113     if not n.mod(4).is_zero():
2114         raise ValueError("the matrix 'M' must be a complex embedding")
2115
2116     Q = QuaternionAlgebra(QQ,-1,-1)
2117     i,j,k = Q.gens()
2118
2119     # Go top-left to bottom-right (reading order), converting every
2120     # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2121     # quaternion block.
2122     elements = []
2123     for l in xrange(n/4):
2124         for m in xrange(n/4):
2125             submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
2126             if submat[0,0] != submat[1,1].conjugate():
2127                 raise ValueError('bad on-diagonal submatrix')
2128             if submat[0,1] != -submat[1,0].conjugate():
2129                 raise ValueError('bad off-diagonal submatrix')
2130             z  = submat[0,0].real() + submat[0,0].imag()*i
2131             z += submat[0,1].real()*j + submat[0,1].imag()*k
2132             elements.append(z)
2133
2134     return matrix(Q, n/4, elements)
2135
2136
2137 # The usual inner product on R^n.
2138 def _usual_ip(x,y):
2139     return x.vector().inner_product(y.vector())
2140
2141 # The inner product used for the real symmetric simple EJA.
2142 # We keep it as a separate function because e.g. the complex
2143 # algebra uses the same inner product, except divided by 2.
2144 def _matrix_ip(X,Y):
2145     X_mat = X.natural_representation()
2146     Y_mat = Y.natural_representation()
2147     return (X_mat*Y_mat).trace()
2148
2149
2150 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
2151     """
2152     The rank-n simple EJA consisting of real symmetric n-by-n
2153     matrices, the usual symmetric Jordan product, and the trace inner
2154     product. It has dimension `(n^2 + n)/2` over the reals.
2155
2156     EXAMPLES::
2157
2158         sage: J = RealSymmetricEJA(2)
2159         sage: e0, e1, e2 = J.gens()
2160         sage: e0*e0
2161         e0
2162         sage: e1*e1
2163         e0 + e2
2164         sage: e2*e2
2165         e2
2166
2167     TESTS:
2168
2169     The degree of this algebra is `(n^2 + n) / 2`::
2170
2171         sage: set_random_seed()
2172         sage: n = ZZ.random_element(1,5)
2173         sage: J = RealSymmetricEJA(n)
2174         sage: J.degree() == (n^2 + n)/2
2175         True
2176
2177     The Jordan multiplication is what we think it is::
2178
2179         sage: set_random_seed()
2180         sage: n = ZZ.random_element(1,5)
2181         sage: J = RealSymmetricEJA(n)
2182         sage: x = J.random_element()
2183         sage: y = J.random_element()
2184         sage: actual = (x*y).natural_representation()
2185         sage: X = x.natural_representation()
2186         sage: Y = y.natural_representation()
2187         sage: expected = (X*Y + Y*X)/2
2188         sage: actual == expected
2189         True
2190         sage: J(expected) == x*y
2191         True
2192
2193     """
2194     @staticmethod
2195     def __classcall_private__(cls, n, field=QQ):
2196         S = _real_symmetric_basis(n, field=field)
2197         (Qs, T) = _multiplication_table_from_matrix_basis(S)
2198
2199         fdeja = super(RealSymmetricEJA, cls)
2200         return fdeja.__classcall_private__(cls,
2201                                            field,
2202                                            Qs,
2203                                            rank=n,
2204                                            natural_basis=T)
2205
2206     def inner_product(self, x, y):
2207         return _matrix_ip(x,y)
2208
2209
2210 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
2211     """
2212     The rank-n simple EJA consisting of complex Hermitian n-by-n
2213     matrices over the real numbers, the usual symmetric Jordan product,
2214     and the real-part-of-trace inner product. It has dimension `n^2` over
2215     the reals.
2216
2217     TESTS:
2218
2219     The degree of this algebra is `n^2`::
2220
2221         sage: set_random_seed()
2222         sage: n = ZZ.random_element(1,5)
2223         sage: J = ComplexHermitianEJA(n)
2224         sage: J.degree() == n^2
2225         True
2226
2227     The Jordan multiplication is what we think it is::
2228
2229         sage: set_random_seed()
2230         sage: n = ZZ.random_element(1,5)
2231         sage: J = ComplexHermitianEJA(n)
2232         sage: x = J.random_element()
2233         sage: y = J.random_element()
2234         sage: actual = (x*y).natural_representation()
2235         sage: X = x.natural_representation()
2236         sage: Y = y.natural_representation()
2237         sage: expected = (X*Y + Y*X)/2
2238         sage: actual == expected
2239         True
2240         sage: J(expected) == x*y
2241         True
2242
2243     """
2244     @staticmethod
2245     def __classcall_private__(cls, n, field=QQ):
2246         S = _complex_hermitian_basis(n)
2247         (Qs, T) = _multiplication_table_from_matrix_basis(S)
2248
2249         fdeja = super(ComplexHermitianEJA, cls)
2250         return fdeja.__classcall_private__(cls,
2251                                            field,
2252                                            Qs,
2253                                            rank=n,
2254                                            natural_basis=T)
2255
2256     def inner_product(self, x, y):
2257         # Since a+bi on the diagonal is represented as
2258         #
2259         #   a + bi  = [  a  b  ]
2260         #             [ -b  a  ],
2261         #
2262         # we'll double-count the "a" entries if we take the trace of
2263         # the embedding.
2264         return _matrix_ip(x,y)/2
2265
2266
2267 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
2268     """
2269     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2270     matrices, the usual symmetric Jordan product, and the
2271     real-part-of-trace inner product. It has dimension `2n^2 - n` over
2272     the reals.
2273
2274     TESTS:
2275
2276     The degree of this algebra is `n^2`::
2277
2278         sage: set_random_seed()
2279         sage: n = ZZ.random_element(1,5)
2280         sage: J = QuaternionHermitianEJA(n)
2281         sage: J.degree() == 2*(n^2) - n
2282         True
2283
2284     The Jordan multiplication is what we think it is::
2285
2286         sage: set_random_seed()
2287         sage: n = ZZ.random_element(1,5)
2288         sage: J = QuaternionHermitianEJA(n)
2289         sage: x = J.random_element()
2290         sage: y = J.random_element()
2291         sage: actual = (x*y).natural_representation()
2292         sage: X = x.natural_representation()
2293         sage: Y = y.natural_representation()
2294         sage: expected = (X*Y + Y*X)/2
2295         sage: actual == expected
2296         True
2297         sage: J(expected) == x*y
2298         True
2299
2300     """
2301     @staticmethod
2302     def __classcall_private__(cls, n, field=QQ):
2303         S = _quaternion_hermitian_basis(n)
2304         (Qs, T) = _multiplication_table_from_matrix_basis(S)
2305
2306         fdeja = super(QuaternionHermitianEJA, cls)
2307         return fdeja.__classcall_private__(cls,
2308                                            field,
2309                                            Qs,
2310                                            rank=n,
2311                                            natural_basis=T)
2312
2313     def inner_product(self, x, y):
2314         # Since a+bi+cj+dk on the diagonal is represented as
2315         #
2316         #   a + bi +cj + dk = [  a  b  c  d]
2317         #                     [ -b  a -d  c]
2318         #                     [ -c  d  a -b]
2319         #                     [ -d -c  b  a],
2320         #
2321         # we'll quadruple-count the "a" entries if we take the trace of
2322         # the embedding.
2323         return _matrix_ip(x,y)/4
2324
2325
2326 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
2327     """
2328     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2329     with the usual inner product and jordan product ``x*y =
2330     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2331     the reals.
2332
2333     EXAMPLES:
2334
2335     This multiplication table can be verified by hand::
2336
2337         sage: J = JordanSpinEJA(4)
2338         sage: e0,e1,e2,e3 = J.gens()
2339         sage: e0*e0
2340         e0
2341         sage: e0*e1
2342         e1
2343         sage: e0*e2
2344         e2
2345         sage: e0*e3
2346         e3
2347         sage: e1*e2
2348         0
2349         sage: e1*e3
2350         0
2351         sage: e2*e3
2352         0
2353
2354     """
2355     @staticmethod
2356     def __classcall_private__(cls, n, field=QQ):
2357         Qs = []
2358         id_matrix = identity_matrix(field, n)
2359         for i in xrange(n):
2360             ei = id_matrix.column(i)
2361             Qi = zero_matrix(field, n)
2362             Qi.set_row(0, ei)
2363             Qi.set_column(0, ei)
2364             Qi += diagonal_matrix(n, [ei[0]]*n)
2365             # The addition of the diagonal matrix adds an extra ei[0] in the
2366             # upper-left corner of the matrix.
2367             Qi[0,0] = Qi[0,0] * ~field(2)
2368             Qs.append(Qi)
2369
2370         # The rank of the spin algebra is two, unless we're in a
2371         # one-dimensional ambient space (because the rank is bounded by
2372         # the ambient dimension).
2373         fdeja = super(JordanSpinEJA, cls)
2374         return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2))
2375
2376     def inner_product(self, x, y):
2377         return _usual_ip(x,y)