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