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
  15 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
  16     @staticmethod
  17     def __classcall_private__(cls,
  18                               field,
  19                               mult_table,
  20                               names='e',
  21                               assume_associative=False,
  22                               category=None,
  23                               rank=None,
  24                               natural_basis=None):
  25         n = len(mult_table)
  26         mult_table = [b.base_extend(field) for b in mult_table]
  27         for b in mult_table:
  28             b.set_immutable()
  29             if not (is_Matrix(b) and b.dimensions() == (n, n)):
  30                 raise ValueError("input is not a multiplication table")
  31         mult_table = tuple(mult_table)
  32
  33         cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
  34         cat.or_subcategory(category)
  35         if assume_associative:
  36             cat = cat.Associative()
  37
  38         names = normalize_names(n, names)
  39
  40         fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls)
  41         return fda.__classcall__(cls,
  42                                  field,
  43                                  mult_table,
  44                                  assume_associative=assume_associative,
  45                                  names=names,
  46                                  category=cat,
  47                                  rank=rank,
  48                                  natural_basis=natural_basis)
  49
  50
  51     def __init__(self,
  52                  field,
  53                  mult_table,
  54                  names='e',
  55                  assume_associative=False,
  56                  category=None,
  57                  rank=None,
  58                  natural_basis=None):
  59         """
  60         EXAMPLES:
  61
  62         By definition, Jordan multiplication commutes::
  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*y == y*x
  69             True
  70
  71         """
  72         self._charpoly = None # for caching
  73         self._rank = rank
  74         self._natural_basis = natural_basis
  75         self._multiplication_table = mult_table
  76         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
  77         fda.__init__(field,
  78                      mult_table,
  79                      names=names,
  80                      category=category)
  81
  82
  83     def _repr_(self):
  84         """
  85         Return a string representation of ``self``.
  86         """
  87         fmt = "Euclidean Jordan algebra of degree {} over {}"
  88         return fmt.format(self.degree(), self.base_ring())
  89
  90
  91
  92     def characteristic_polynomial(self):
  93         """
  94         EXAMPLES:
  95
  96         The characteristic polynomial in the spin algebra is given in
  97         Alizadeh, Example 11.11::
  98
  99             sage: J = JordanSpinEJA(3)
 100             sage: p = J.characteristic_polynomial(); p
 101             X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
 102             sage: xvec = J.one().vector()
 103             sage: p(*xvec)
 104             t^2 - 2*t + 1
 105
 106         """
 107         if self._charpoly is not None:
 108             return self._charpoly
 109
 110         r = self.rank()
 111         n = self.dimension()
 112
 113         # Now switch to the polynomial rings.
 114         names = ['X' + str(i) for i in range(1,n+1)]
 115         R = PolynomialRing(self.base_ring(), names)
 116         J = FiniteDimensionalEuclideanJordanAlgebra(R,
 117                                                     self._multiplication_table,
 118                                                     rank=r)
 119
 120         def e(k):
 121             # The coordinates of e_k with respect to the basis B.
 122             # But, the e_k are elements of B...
 123             return identity_matrix(J.base_ring(), n).column(k-1).column()
 124
 125         # A matrix implementation 1
 126         x = J(vector(R, R.gens()))
 127         l1 = [column_matrix((x**k).vector()) for k in range(r)]
 128         l2 = [e(k) for k in range(r+1, n+1)]
 129         A_of_x = block_matrix(1, n, (l1 + l2))
 130         xr = (x**r).vector()
 131         a = []
 132         denominator = A_of_x.det() # This is constant
 133         for i in range(n):
 134             A_cols = A_of_x.columns()
 135             A_cols[i] = xr
 136             numerator = column_matrix(A_of_x.base_ring(), A_cols).det()
 137             ai = numerator/denominator
 138             a.append(ai)
 139
 140         # We go to a bit of trouble here to reorder the
 141         # indeterminates, so that it's easier to evaluate the
 142         # characteristic polynomial at x's coordinates and get back
 143         # something in terms of t, which is what we want.
 144         S = PolynomialRing(self.base_ring(),'t')
 145         t = S.gen(0)
 146         S = PolynomialRing(S, R.variable_names())
 147         t = S(t)
 148
 149         # We're relying on the theory here to ensure that each entry
 150         # a[i] is indeed back in R, and the added negative signs are
 151         # to make the whole expression sum to zero.
 152         a = [R(-ai) for ai in a] # corresponds to powerx x^0 through x^(r-1)
 153
 154         # Note: all entries past the rth should be zero. The
 155         # coefficient of the highest power (x^r) is 1, but it doesn't
 156         # appear in the solution vector which contains coefficients
 157         # for the other powers (to make them sum to x^r).
 158         if (r < n):
 159             a[r] = 1 # corresponds to x^r
 160         else:
 161             # When the rank is equal to the dimension, trying to
 162             # assign a[r] goes out-of-bounds.
 163             a.append(1) # corresponds to x^r
 164
 165         self._charpoly = sum( a[k]*(t**k) for k in range(len(a)) )
 166         return self._charpoly
 167
 168
 169     def inner_product(self, x, y):
 170         """
 171         The inner product associated with this Euclidean Jordan algebra.
 172
 173         Defaults to the trace inner product, but can be overridden by
 174         subclasses if they are sure that the necessary properties are
 175         satisfied.
 176
 177         EXAMPLES:
 178
 179         The inner product must satisfy its axiom for this algebra to truly
 180         be a Euclidean Jordan Algebra::
 181
 182             sage: set_random_seed()
 183             sage: J = random_eja()
 184             sage: x = J.random_element()
 185             sage: y = J.random_element()
 186             sage: z = J.random_element()
 187             sage: (x*y).inner_product(z) == y.inner_product(x*z)
 188             True
 189
 190         """
 191         if (not x in self) or (not y in self):
 192             raise TypeError("arguments must live in this algebra")
 193         return x.trace_inner_product(y)
 194
 195
 196     def natural_basis(self):
 197         """
 198         Return a more-natural representation of this algebra's basis.
 199
 200         Every finite-dimensional Euclidean Jordan Algebra is a direct
 201         sum of five simple algebras, four of which comprise Hermitian
 202         matrices. This method returns the original "natural" basis
 203         for our underlying vector space. (Typically, the natural basis
 204         is used to construct the multiplication table in the first place.)
 205
 206         Note that this will always return a matrix. The standard basis
 207         in `R^n` will be returned as `n`-by-`1` column matrices.
 208
 209         EXAMPLES::
 210
 211             sage: J = RealSymmetricEJA(2)
 212             sage: J.basis()
 213             Family (e0, e1, e2)
 214             sage: J.natural_basis()
 215             (
 216             [1 0]  [0 1]  [0 0]
 217             [0 0], [1 0], [0 1]
 218             )
 219
 220         ::
 221
 222             sage: J = JordanSpinEJA(2)
 223             sage: J.basis()
 224             Family (e0, e1)
 225             sage: J.natural_basis()
 226             (
 227             [1]  [0]
 228             [0], [1]
 229             )
 230
 231         """
 232         if self._natural_basis is None:
 233             return tuple( b.vector().column() for b in self.basis() )
 234         else:
 235             return self._natural_basis
 236
 237
 238     def rank(self):
 239         """
 240         Return the rank of this EJA.
 241         """
 242         if self._rank is None:
 243             raise ValueError("no rank specified at genesis")
 244         else:
 245             return self._rank
 246
 247
 248     class Element(FiniteDimensionalAlgebraElement):
 249         """
 250         An element of a Euclidean Jordan algebra.
 251         """
 252
 253         def __init__(self, A, elt=None):
 254             """
 255             EXAMPLES:
 256
 257             The identity in `S^n` is converted to the identity in the EJA::
 258
 259                 sage: J = RealSymmetricEJA(3)
 260                 sage: I = identity_matrix(QQ,3)
 261                 sage: J(I) == J.one()
 262                 True
 263
 264             This skew-symmetric matrix can't be represented in the EJA::
 265
 266                 sage: J = RealSymmetricEJA(3)
 267                 sage: A = matrix(QQ,3, lambda i,j: i-j)
 268                 sage: J(A)
 269                 Traceback (most recent call last):
 270                 ...
 271                 ArithmeticError: vector is not in free module
 272
 273             """
 274             # Goal: if we're given a matrix, and if it lives in our
 275             # parent algebra's "natural ambient space," convert it
 276             # into an algebra element.
 277             #
 278             # The catch is, we make a recursive call after converting
 279             # the given matrix into a vector that lives in the algebra.
 280             # This we need to try the parent class initializer first,
 281             # to avoid recursing forever if we're given something that
 282             # already fits into the algebra, but also happens to live
 283             # in the parent's "natural ambient space" (this happens with
 284             # vectors in R^n).
 285             try:
 286                 FiniteDimensionalAlgebraElement.__init__(self, A, elt)
 287             except ValueError:
 288                 natural_basis = A.natural_basis()
 289                 if elt in natural_basis[0].matrix_space():
 290                     # Thanks for nothing! Matrix spaces aren't vector
 291                     # spaces in Sage, so we have to figure out its
 292                     # natural-basis coordinates ourselves.
 293                     V = VectorSpace(elt.base_ring(), elt.nrows()**2)
 294                     W = V.span( _mat2vec(s) for s in natural_basis )
 295                     coords =  W.coordinates(_mat2vec(elt))
 296                     FiniteDimensionalAlgebraElement.__init__(self, A, coords)
 297
 298         def __pow__(self, n):
 299             """
 300             Return ``self`` raised to the power ``n``.
 301
 302             Jordan algebras are always power-associative; see for
 303             example Faraut and Koranyi, Proposition II.1.2 (ii).
 304
 305             .. WARNING:
 306
 307                 We have to override this because our superclass uses row vectors
 308                 instead of column vectors! We, on the other hand, assume column
 309                 vectors everywhere.
 310
 311             EXAMPLES::
 312
 313                 sage: set_random_seed()
 314                 sage: x = random_eja().random_element()
 315                 sage: x.operator_matrix()*x.vector() == (x^2).vector()
 316                 True
 317
 318             A few examples of power-associativity::
 319
 320                 sage: set_random_seed()
 321                 sage: x = random_eja().random_element()
 322                 sage: x*(x*x)*(x*x) == x^5
 323                 True
 324                 sage: (x*x)*(x*x*x) == x^5
 325                 True
 326
 327             We also know that powers operator-commute (Koecher, Chapter
 328             III, Corollary 1)::
 329
 330                 sage: set_random_seed()
 331                 sage: x = random_eja().random_element()
 332                 sage: m = ZZ.random_element(0,10)
 333                 sage: n = ZZ.random_element(0,10)
 334                 sage: Lxm = (x^m).operator_matrix()
 335                 sage: Lxn = (x^n).operator_matrix()
 336                 sage: Lxm*Lxn == Lxn*Lxm
 337                 True
 338
 339             """
 340             A = self.parent()
 341             if n == 0:
 342                 return A.one()
 343             elif n == 1:
 344                 return self
 345             else:
 346                 return A( (self.operator_matrix()**(n-1))*self.vector() )
 347
 348
 349         def apply_univariate_polynomial(self, p):
 350             """
 351             Apply the univariate polynomial ``p`` to this element.
 352
 353             A priori, SageMath won't allow us to apply a univariate
 354             polynomial to an element of an EJA, because we don't know
 355             that EJAs are rings (they are usually not associative). Of
 356             course, we know that EJAs are power-associative, so the
 357             operation is ultimately kosher. This function sidesteps
 358             the CAS to get the answer we want and expect.
 359
 360             EXAMPLES::
 361
 362                 sage: R = PolynomialRing(QQ, 't')
 363                 sage: t = R.gen(0)
 364                 sage: p = t^4 - t^3 + 5*t - 2
 365                 sage: J = RealCartesianProductEJA(5)
 366                 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
 367                 True
 368
 369             TESTS:
 370
 371             We should always get back an element of the algebra::
 372
 373                 sage: set_random_seed()
 374                 sage: p = PolynomialRing(QQ, 't').random_element()
 375                 sage: J = random_eja()
 376                 sage: x = J.random_element()
 377                 sage: x.apply_univariate_polynomial(p) in J
 378                 True
 379
 380             """
 381             if len(p.variables()) > 1:
 382                 raise ValueError("not a univariate polynomial")
 383             P = self.parent()
 384             R = P.base_ring()
 385             # Convert the coeficcients to the parent's base ring,
 386             # because a priori they might live in an (unnecessarily)
 387             # larger ring for which P.sum() would fail below.
 388             cs = [ R(c) for c in p.coefficients(sparse=False) ]
 389             return P.sum( cs[k]*(self**k) for k in range(len(cs)) )
 390
 391
 392         def characteristic_polynomial(self):
 393             """
 394             Return the characteristic polynomial of this element.
 395
 396             EXAMPLES:
 397
 398             The rank of `R^3` is three, and the minimal polynomial of
 399             the identity element is `(t-1)` from which it follows that
 400             the characteristic polynomial should be `(t-1)^3`::
 401
 402                 sage: J = RealCartesianProductEJA(3)
 403                 sage: J.one().characteristic_polynomial()
 404                 t^3 - 3*t^2 + 3*t - 1
 405
 406             Likewise, the characteristic of the zero element in the
 407             rank-three algebra `R^{n}` should be `t^{3}`::
 408
 409                 sage: J = RealCartesianProductEJA(3)
 410                 sage: J.zero().characteristic_polynomial()
 411                 t^3
 412
 413             The characteristic polynomial of an element should evaluate
 414             to zero on that element::
 415
 416                 sage: set_random_seed()
 417                 sage: x = RealCartesianProductEJA(3).random_element()
 418                 sage: p = x.characteristic_polynomial()
 419                 sage: x.apply_univariate_polynomial(p)
 420                 0
 421
 422             """
 423             p = self.parent().characteristic_polynomial()
 424             return p(*self.vector())
 425
 426
 427         def inner_product(self, other):
 428             """
 429             Return the parent algebra's inner product of myself and ``other``.
 430
 431             EXAMPLES:
 432
 433             The inner product in the Jordan spin algebra is the usual
 434             inner product on `R^n` (this example only works because the
 435             basis for the Jordan algebra is the standard basis in `R^n`)::
 436
 437                 sage: J = JordanSpinEJA(3)
 438                 sage: x = vector(QQ,[1,2,3])
 439                 sage: y = vector(QQ,[4,5,6])
 440                 sage: x.inner_product(y)
 441                 32
 442                 sage: J(x).inner_product(J(y))
 443                 32
 444
 445             The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
 446             multiplication is the usual matrix multiplication in `S^n`,
 447             so the inner product of the identity matrix with itself
 448             should be the `n`::
 449
 450                 sage: J = RealSymmetricEJA(3)
 451                 sage: J.one().inner_product(J.one())
 452                 3
 453
 454             Likewise, the inner product on `C^n` is `<X,Y> =
 455             Re(trace(X*Y))`, where we must necessarily take the real
 456             part because the product of Hermitian matrices may not be
 457             Hermitian::
 458
 459                 sage: J = ComplexHermitianEJA(3)
 460                 sage: J.one().inner_product(J.one())
 461                 3
 462
 463             Ditto for the quaternions::
 464
 465                 sage: J = QuaternionHermitianEJA(3)
 466                 sage: J.one().inner_product(J.one())
 467                 3
 468
 469             TESTS:
 470
 471             Ensure that we can always compute an inner product, and that
 472             it gives us back a real number::
 473
 474                 sage: set_random_seed()
 475                 sage: J = random_eja()
 476                 sage: x = J.random_element()
 477                 sage: y = J.random_element()
 478                 sage: x.inner_product(y) in RR
 479                 True
 480
 481             """
 482             P = self.parent()
 483             if not other in P:
 484                 raise TypeError("'other' must live in the same algebra")
 485
 486             return P.inner_product(self, other)
 487
 488
 489         def operator_commutes_with(self, other):
 490             """
 491             Return whether or not this element operator-commutes
 492             with ``other``.
 493
 494             EXAMPLES:
 495
 496             The definition of a Jordan algebra says that any element
 497             operator-commutes with its square::
 498
 499                 sage: set_random_seed()
 500                 sage: x = random_eja().random_element()
 501                 sage: x.operator_commutes_with(x^2)
 502                 True
 503
 504             TESTS:
 505
 506             Test Lemma 1 from Chapter III of Koecher::
 507
 508                 sage: set_random_seed()
 509                 sage: J = random_eja()
 510                 sage: u = J.random_element()
 511                 sage: v = J.random_element()
 512                 sage: lhs = u.operator_commutes_with(u*v)
 513                 sage: rhs = v.operator_commutes_with(u^2)
 514                 sage: lhs == rhs
 515                 True
 516
 517             """
 518             if not other in self.parent():
 519                 raise TypeError("'other' must live in the same algebra")
 520
 521             A = self.operator_matrix()
 522             B = other.operator_matrix()
 523             return (A*B == B*A)
 524
 525
 526         def det(self):
 527             """
 528             Return my determinant, the product of my eigenvalues.
 529
 530             EXAMPLES::
 531
 532                 sage: J = JordanSpinEJA(2)
 533                 sage: e0,e1 = J.gens()
 534                 sage: x = e0 + e1
 535                 sage: x.det()
 536                 0
 537                 sage: J = JordanSpinEJA(3)
 538                 sage: e0,e1,e2 = J.gens()
 539                 sage: x = e0 + e1 + e2
 540                 sage: x.det()
 541                 -1
 542
 543             """
 544             cs = self.characteristic_polynomial().coefficients(sparse=False)
 545             r = len(cs) - 1
 546             if r >= 0:
 547                 return cs[0] * (-1)**r
 548             else:
 549                 raise ValueError('charpoly had no coefficients')
 550
 551
 552         def inverse(self):
 553             """
 554             Return the Jordan-multiplicative inverse of this element.
 555
 556             We can't use the superclass method because it relies on the
 557             algebra being associative.
 558
 559             EXAMPLES:
 560
 561             The inverse in the spin factor algebra is given in Alizadeh's
 562             Example 11.11::
 563
 564                 sage: set_random_seed()
 565                 sage: n = ZZ.random_element(1,10)
 566                 sage: J = JordanSpinEJA(n)
 567                 sage: x = J.random_element()
 568                 sage: while x.is_zero():
 569                 ....:     x = J.random_element()
 570                 sage: x_vec = x.vector()
 571                 sage: x0 = x_vec[0]
 572                 sage: x_bar = x_vec[1:]
 573                 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
 574                 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
 575                 sage: x_inverse = coeff*inv_vec
 576                 sage: x.inverse() == J(x_inverse)
 577                 True
 578
 579             TESTS:
 580
 581             The identity element is its own inverse::
 582
 583                 sage: set_random_seed()
 584                 sage: J = random_eja()
 585                 sage: J.one().inverse() == J.one()
 586                 True
 587
 588             If an element has an inverse, it acts like one. TODO: this
 589             can be a lot less ugly once ``is_invertible`` doesn't crash
 590             on irregular elements::
 591
 592                 sage: set_random_seed()
 593                 sage: J = random_eja()
 594                 sage: x = J.random_element()
 595                 sage: try:
 596                 ....:     x.inverse()*x == J.one()
 597                 ....: except:
 598                 ....:     True
 599                 True
 600
 601             """
 602             if self.parent().is_associative():
 603                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 604                 return elt.inverse()
 605
 606             # TODO: we can do better once the call to is_invertible()
 607             # doesn't crash on irregular elements.
 608             #if not self.is_invertible():
 609             #    raise ValueError('element is not invertible')
 610
 611             # We do this a little different than the usual recursive
 612             # call to a finite-dimensional algebra element, because we
 613             # wind up with an inverse that lives in the subalgebra and
 614             # we need information about the parent to convert it back.
 615             V = self.span_of_powers()
 616             assoc_subalg = self.subalgebra_generated_by()
 617             # Mis-design warning: the basis used for span_of_powers()
 618             # and subalgebra_generated_by() must be the same, and in
 619             # the same order!
 620             elt = assoc_subalg(V.coordinates(self.vector()))
 621
 622             # This will be in the subalgebra's coordinates...
 623             fda_elt = FiniteDimensionalAlgebraElement(assoc_subalg, elt)
 624             subalg_inverse = fda_elt.inverse()
 625
 626             # So we have to convert back...
 627             basis = [ self.parent(v) for v in V.basis() ]
 628             pairs = zip(subalg_inverse.vector(), basis)
 629             return self.parent().linear_combination(pairs)
 630
 631
 632         def is_invertible(self):
 633             """
 634             Return whether or not this element is invertible.
 635
 636             We can't use the superclass method because it relies on
 637             the algebra being associative.
 638
 639             ALGORITHM:
 640
 641             The usual way to do this is to check if the determinant is
 642             zero, but we need the characteristic polynomial for the
 643             determinant. The minimal polynomial is a lot easier to get,
 644             so we use Corollary 2 in Chapter V of Koecher to check
 645             whether or not the paren't algebra's zero element is a root
 646             of this element's minimal polynomial.
 647
 648             TESTS:
 649
 650             The identity element is always invertible::
 651
 652                 sage: set_random_seed()
 653                 sage: J = random_eja()
 654                 sage: J.one().is_invertible()
 655                 True
 656
 657             The zero element is never invertible::
 658
 659                 sage: set_random_seed()
 660                 sage: J = random_eja()
 661                 sage: J.zero().is_invertible()
 662                 False
 663
 664             """
 665             zero = self.parent().zero()
 666             p = self.minimal_polynomial()
 667             return not (p(zero) == zero)
 668
 669
 670         def is_nilpotent(self):
 671             """
 672             Return whether or not some power of this element is zero.
 673
 674             The superclass method won't work unless we're in an
 675             associative algebra, and we aren't. However, we generate
 676             an assocoative subalgebra and we're nilpotent there if and
 677             only if we're nilpotent here (probably).
 678
 679             TESTS:
 680
 681             The identity element is never nilpotent::
 682
 683                 sage: set_random_seed()
 684                 sage: random_eja().one().is_nilpotent()
 685                 False
 686
 687             The additive identity is always nilpotent::
 688
 689                 sage: set_random_seed()
 690                 sage: random_eja().zero().is_nilpotent()
 691                 True
 692
 693             """
 694             # The element we're going to call "is_nilpotent()" on.
 695             # Either myself, interpreted as an element of a finite-
 696             # dimensional algebra, or an element of an associative
 697             # subalgebra.
 698             elt = None
 699
 700             if self.parent().is_associative():
 701                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 702             else:
 703                 V = self.span_of_powers()
 704                 assoc_subalg = self.subalgebra_generated_by()
 705                 # Mis-design warning: the basis used for span_of_powers()
 706                 # and subalgebra_generated_by() must be the same, and in
 707                 # the same order!
 708                 elt = assoc_subalg(V.coordinates(self.vector()))
 709
 710             # Recursive call, but should work since elt lives in an
 711             # associative algebra.
 712             return elt.is_nilpotent()
 713
 714
 715         def is_regular(self):
 716             """
 717             Return whether or not this is a regular element.
 718
 719             EXAMPLES:
 720
 721             The identity element always has degree one, but any element
 722             linearly-independent from it is regular::
 723
 724                 sage: J = JordanSpinEJA(5)
 725                 sage: J.one().is_regular()
 726                 False
 727                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 728                 sage: for x in J.gens():
 729                 ....:     (J.one() + x).is_regular()
 730                 False
 731                 True
 732                 True
 733                 True
 734                 True
 735
 736             """
 737             return self.degree() == self.parent().rank()
 738
 739
 740         def degree(self):
 741             """
 742             Compute the degree of this element the straightforward way
 743             according to the definition; by appending powers to a list
 744             and figuring out its dimension (that is, whether or not
 745             they're linearly dependent).
 746
 747             EXAMPLES::
 748
 749                 sage: J = JordanSpinEJA(4)
 750                 sage: J.one().degree()
 751                 1
 752                 sage: e0,e1,e2,e3 = J.gens()
 753                 sage: (e0 - e1).degree()
 754                 2
 755
 756             In the spin factor algebra (of rank two), all elements that
 757             aren't multiples of the identity are regular::
 758
 759                 sage: set_random_seed()
 760                 sage: n = ZZ.random_element(1,10)
 761                 sage: J = JordanSpinEJA(n)
 762                 sage: x = J.random_element()
 763                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 764                 True
 765
 766             """
 767             return self.span_of_powers().dimension()
 768
 769
 770         def minimal_polynomial(self):
 771             """
 772             Return the minimal polynomial of this element,
 773             as a function of the variable `t`.
 774
 775             ALGORITHM:
 776
 777             We restrict ourselves to the associative subalgebra
 778             generated by this element, and then return the minimal
 779             polynomial of this element's operator matrix (in that
 780             subalgebra). This works by Baes Proposition 2.3.16.
 781
 782             TESTS:
 783
 784             The minimal polynomial of the identity and zero elements are
 785             always the same::
 786
 787                 sage: set_random_seed()
 788                 sage: J = random_eja()
 789                 sage: J.one().minimal_polynomial()
 790                 t - 1
 791                 sage: J.zero().minimal_polynomial()
 792                 t
 793
 794             The degree of an element is (by one definition) the degree
 795             of its minimal polynomial::
 796
 797                 sage: set_random_seed()
 798                 sage: x = random_eja().random_element()
 799                 sage: x.degree() == x.minimal_polynomial().degree()
 800                 True
 801
 802             The minimal polynomial and the characteristic polynomial coincide
 803             and are known (see Alizadeh, Example 11.11) for all elements of
 804             the spin factor algebra that aren't scalar multiples of the
 805             identity::
 806
 807                 sage: set_random_seed()
 808                 sage: n = ZZ.random_element(2,10)
 809                 sage: J = JordanSpinEJA(n)
 810                 sage: y = J.random_element()
 811                 sage: while y == y.coefficient(0)*J.one():
 812                 ....:     y = J.random_element()
 813                 sage: y0 = y.vector()[0]
 814                 sage: y_bar = y.vector()[1:]
 815                 sage: actual = y.minimal_polynomial()
 816                 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
 817                 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
 818                 sage: bool(actual == expected)
 819                 True
 820
 821             The minimal polynomial should always kill its element::
 822
 823                 sage: set_random_seed()
 824                 sage: x = random_eja().random_element()
 825                 sage: p = x.minimal_polynomial()
 826                 sage: x.apply_univariate_polynomial(p)
 827                 0
 828
 829             """
 830             V = self.span_of_powers()
 831             assoc_subalg = self.subalgebra_generated_by()
 832             # Mis-design warning: the basis used for span_of_powers()
 833             # and subalgebra_generated_by() must be the same, and in
 834             # the same order!
 835             elt = assoc_subalg(V.coordinates(self.vector()))
 836
 837             # We get back a symbolic polynomial in 'x' but want a real
 838             # polynomial in 't'.
 839             p_of_x = elt.operator_matrix().minimal_polynomial()
 840             return p_of_x.change_variable_name('t')
 841
 842
 843         def natural_representation(self):
 844             """
 845             Return a more-natural representation of this element.
 846
 847             Every finite-dimensional Euclidean Jordan Algebra is a
 848             direct sum of five simple algebras, four of which comprise
 849             Hermitian matrices. This method returns the original
 850             "natural" representation of this element as a Hermitian
 851             matrix, if it has one. If not, you get the usual representation.
 852
 853             EXAMPLES::
 854
 855                 sage: J = ComplexHermitianEJA(3)
 856                 sage: J.one()
 857                 e0 + e5 + e8
 858                 sage: J.one().natural_representation()
 859                 [1 0 0 0 0 0]
 860                 [0 1 0 0 0 0]
 861                 [0 0 1 0 0 0]
 862                 [0 0 0 1 0 0]
 863                 [0 0 0 0 1 0]
 864                 [0 0 0 0 0 1]
 865
 866             ::
 867
 868                 sage: J = QuaternionHermitianEJA(3)
 869                 sage: J.one()
 870                 e0 + e9 + e14
 871                 sage: J.one().natural_representation()
 872                 [1 0 0 0 0 0 0 0 0 0 0 0]
 873                 [0 1 0 0 0 0 0 0 0 0 0 0]
 874                 [0 0 1 0 0 0 0 0 0 0 0 0]
 875                 [0 0 0 1 0 0 0 0 0 0 0 0]
 876                 [0 0 0 0 1 0 0 0 0 0 0 0]
 877                 [0 0 0 0 0 1 0 0 0 0 0 0]
 878                 [0 0 0 0 0 0 1 0 0 0 0 0]
 879                 [0 0 0 0 0 0 0 1 0 0 0 0]
 880                 [0 0 0 0 0 0 0 0 1 0 0 0]
 881                 [0 0 0 0 0 0 0 0 0 1 0 0]
 882                 [0 0 0 0 0 0 0 0 0 0 1 0]
 883                 [0 0 0 0 0 0 0 0 0 0 0 1]
 884
 885             """
 886             B = self.parent().natural_basis()
 887             W = B[0].matrix_space()
 888             return W.linear_combination(zip(self.vector(), B))
 889
 890
 891         def operator_matrix(self):
 892             """
 893             Return the matrix that represents left- (or right-)
 894             multiplication by this element in the parent algebra.
 895
 896             We have to override this because the superclass method
 897             returns a matrix that acts on row vectors (that is, on
 898             the right).
 899
 900             EXAMPLES:
 901
 902             Test the first polarization identity from my notes, Koecher Chapter
 903             III, or from Baes (2.3)::
 904
 905                 sage: set_random_seed()
 906                 sage: J = random_eja()
 907                 sage: x = J.random_element()
 908                 sage: y = J.random_element()
 909                 sage: Lx = x.operator_matrix()
 910                 sage: Ly = y.operator_matrix()
 911                 sage: Lxx = (x*x).operator_matrix()
 912                 sage: Lxy = (x*y).operator_matrix()
 913                 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
 914                 True
 915
 916             Test the second polarization identity from my notes or from
 917             Baes (2.4)::
 918
 919                 sage: set_random_seed()
 920                 sage: J = random_eja()
 921                 sage: x = J.random_element()
 922                 sage: y = J.random_element()
 923                 sage: z = J.random_element()
 924                 sage: Lx = x.operator_matrix()
 925                 sage: Ly = y.operator_matrix()
 926                 sage: Lz = z.operator_matrix()
 927                 sage: Lzy = (z*y).operator_matrix()
 928                 sage: Lxy = (x*y).operator_matrix()
 929                 sage: Lxz = (x*z).operator_matrix()
 930                 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
 931                 True
 932
 933             Test the third polarization identity from my notes or from
 934             Baes (2.5)::
 935
 936                 sage: set_random_seed()
 937                 sage: J = random_eja()
 938                 sage: u = J.random_element()
 939                 sage: y = J.random_element()
 940                 sage: z = J.random_element()
 941                 sage: Lu = u.operator_matrix()
 942                 sage: Ly = y.operator_matrix()
 943                 sage: Lz = z.operator_matrix()
 944                 sage: Lzy = (z*y).operator_matrix()
 945                 sage: Luy = (u*y).operator_matrix()
 946                 sage: Luz = (u*z).operator_matrix()
 947                 sage: Luyz = (u*(y*z)).operator_matrix()
 948                 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
 949                 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
 950                 sage: bool(lhs == rhs)
 951                 True
 952
 953             """
 954             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 955             return fda_elt.matrix().transpose()
 956
 957
 958         def quadratic_representation(self, other=None):
 959             """
 960             Return the quadratic representation of this element.
 961
 962             EXAMPLES:
 963
 964             The explicit form in the spin factor algebra is given by
 965             Alizadeh's Example 11.12::
 966
 967                 sage: set_random_seed()
 968                 sage: n = ZZ.random_element(1,10)
 969                 sage: J = JordanSpinEJA(n)
 970                 sage: x = J.random_element()
 971                 sage: x_vec = x.vector()
 972                 sage: x0 = x_vec[0]
 973                 sage: x_bar = x_vec[1:]
 974                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 975                 sage: B = 2*x0*x_bar.row()
 976                 sage: C = 2*x0*x_bar.column()
 977                 sage: D = identity_matrix(QQ, n-1)
 978                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 979                 sage: D = D + 2*x_bar.tensor_product(x_bar)
 980                 sage: Q = block_matrix(2,2,[A,B,C,D])
 981                 sage: Q == x.quadratic_representation()
 982                 True
 983
 984             Test all of the properties from Theorem 11.2 in Alizadeh::
 985
 986                 sage: set_random_seed()
 987                 sage: J = random_eja()
 988                 sage: x = J.random_element()
 989                 sage: y = J.random_element()
 990
 991             Property 1:
 992
 993                 sage: actual = x.quadratic_representation(y)
 994                 sage: expected = ( (x+y).quadratic_representation()
 995                 ....:              -x.quadratic_representation()
 996                 ....:              -y.quadratic_representation() ) / 2
 997                 sage: actual == expected
 998                 True
 999
1000             Property 2:
1001
1002                 sage: alpha = QQ.random_element()
1003                 sage: actual = (alpha*x).quadratic_representation()
1004                 sage: expected = (alpha^2)*x.quadratic_representation()
1005                 sage: actual == expected
1006                 True
1007
1008             Property 5:
1009
1010                 sage: Qy = y.quadratic_representation()
1011                 sage: actual = J(Qy*x.vector()).quadratic_representation()
1012                 sage: expected = Qy*x.quadratic_representation()*Qy
1013                 sage: actual == expected
1014                 True
1015
1016             Property 6:
1017
1018                 sage: k = ZZ.random_element(1,10)
1019                 sage: actual = (x^k).quadratic_representation()
1020                 sage: expected = (x.quadratic_representation())^k
1021                 sage: actual == expected
1022                 True
1023
1024             """
1025             if other is None:
1026                 other=self
1027             elif not other in self.parent():
1028                 raise TypeError("'other' must live in the same algebra")
1029
1030             L = self.operator_matrix()
1031             M = other.operator_matrix()
1032             return ( L*M + M*L - (self*other).operator_matrix() )
1033
1034
1035         def span_of_powers(self):
1036             """
1037             Return the vector space spanned by successive powers of
1038             this element.
1039             """
1040             # The dimension of the subalgebra can't be greater than
1041             # the big algebra, so just put everything into a list
1042             # and let span() get rid of the excess.
1043             #
1044             # We do the extra ambient_vector_space() in case we're messing
1045             # with polynomials and the direct parent is a module.
1046             V = self.vector().parent().ambient_vector_space()
1047             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
1048
1049
1050         def subalgebra_generated_by(self):
1051             """
1052             Return the associative subalgebra of the parent EJA generated
1053             by this element.
1054
1055             TESTS::
1056
1057                 sage: set_random_seed()
1058                 sage: x = random_eja().random_element()
1059                 sage: x.subalgebra_generated_by().is_associative()
1060                 True
1061
1062             Squaring in the subalgebra should be the same thing as
1063             squaring in the superalgebra::
1064
1065                 sage: set_random_seed()
1066                 sage: x = random_eja().random_element()
1067                 sage: u = x.subalgebra_generated_by().random_element()
1068                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
1069                 True
1070
1071             """
1072             # First get the subspace spanned by the powers of myself...
1073             V = self.span_of_powers()
1074             F = self.base_ring()
1075
1076             # Now figure out the entries of the right-multiplication
1077             # matrix for the successive basis elements b0, b1,... of
1078             # that subspace.
1079             mats = []
1080             for b_right in V.basis():
1081                 eja_b_right = self.parent()(b_right)
1082                 b_right_rows = []
1083                 # The first row of the right-multiplication matrix by
1084                 # b1 is what we get if we apply that matrix to b1. The
1085                 # second row of the right multiplication matrix by b1
1086                 # is what we get when we apply that matrix to b2...
1087                 #
1088                 # IMPORTANT: this assumes that all vectors are COLUMN
1089                 # vectors, unlike our superclass (which uses row vectors).
1090                 for b_left in V.basis():
1091                     eja_b_left = self.parent()(b_left)
1092                     # Multiply in the original EJA, but then get the
1093                     # coordinates from the subalgebra in terms of its
1094                     # basis.
1095                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
1096                     b_right_rows.append(this_row)
1097                 b_right_matrix = matrix(F, b_right_rows)
1098                 mats.append(b_right_matrix)
1099
1100             # It's an algebra of polynomials in one element, and EJAs
1101             # are power-associative.
1102             #
1103             # TODO: choose generator names intelligently.
1104             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
1105
1106
1107         def subalgebra_idempotent(self):
1108             """
1109             Find an idempotent in the associative subalgebra I generate
1110             using Proposition 2.3.5 in Baes.
1111
1112             TESTS::
1113
1114                 sage: set_random_seed()
1115                 sage: J = RealCartesianProductEJA(5)
1116                 sage: c = J.random_element().subalgebra_idempotent()
1117                 sage: c^2 == c
1118                 True
1119                 sage: J = JordanSpinEJA(5)
1120                 sage: c = J.random_element().subalgebra_idempotent()
1121                 sage: c^2 == c
1122                 True
1123
1124             """
1125             if self.is_nilpotent():
1126                 raise ValueError("this only works with non-nilpotent elements!")
1127
1128             V = self.span_of_powers()
1129             J = self.subalgebra_generated_by()
1130             # Mis-design warning: the basis used for span_of_powers()
1131             # and subalgebra_generated_by() must be the same, and in
1132             # the same order!
1133             u = J(V.coordinates(self.vector()))
1134
1135             # The image of the matrix of left-u^m-multiplication
1136             # will be minimal for some natural number s...
1137             s = 0
1138             minimal_dim = V.dimension()
1139             for i in xrange(1, V.dimension()):
1140                 this_dim = (u**i).operator_matrix().image().dimension()
1141                 if this_dim < minimal_dim:
1142                     minimal_dim = this_dim
1143                     s = i
1144
1145             # Now minimal_matrix should correspond to the smallest
1146             # non-zero subspace in Baes's (or really, Koecher's)
1147             # proposition.
1148             #
1149             # However, we need to restrict the matrix to work on the
1150             # subspace... or do we? Can't we just solve, knowing that
1151             # A(c) = u^(s+1) should have a solution in the big space,
1152             # too?
1153             #
1154             # Beware, solve_right() means that we're using COLUMN vectors.
1155             # Our FiniteDimensionalAlgebraElement superclass uses rows.
1156             u_next = u**(s+1)
1157             A = u_next.operator_matrix()
1158             c_coordinates = A.solve_right(u_next.vector())
1159
1160             # Now c_coordinates is the idempotent we want, but it's in
1161             # the coordinate system of the subalgebra.
1162             #
1163             # We need the basis for J, but as elements of the parent algebra.
1164             #
1165             basis = [self.parent(v) for v in V.basis()]
1166             return self.parent().linear_combination(zip(c_coordinates, basis))
1167
1168
1169         def trace(self):
1170             """
1171             Return my trace, the sum of my eigenvalues.
1172
1173             EXAMPLES::
1174
1175                 sage: J = JordanSpinEJA(3)
1176                 sage: e0,e1,e2 = J.gens()
1177                 sage: x = e0 + e1 + e2
1178                 sage: x.trace()
1179                 2
1180
1181             """
1182             cs = self.characteristic_polynomial().coefficients(sparse=False)
1183             if len(cs) >= 2:
1184                 return -1*cs[-2]
1185             else:
1186                 raise ValueError('charpoly had fewer than 2 coefficients')
1187
1188
1189         def trace_inner_product(self, other):
1190             """
1191             Return the trace inner product of myself and ``other``.
1192             """
1193             if not other in self.parent():
1194                 raise TypeError("'other' must live in the same algebra")
1195
1196             return (self*other).trace()
1197
1198
1199 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
1200     """
1201     Return the Euclidean Jordan Algebra corresponding to the set
1202     `R^n` under the Hadamard product.
1203
1204     Note: this is nothing more than the Cartesian product of ``n``
1205     copies of the spin algebra. Once Cartesian product algebras
1206     are implemented, this can go.
1207
1208     EXAMPLES:
1209
1210     This multiplication table can be verified by hand::
1211
1212         sage: J = RealCartesianProductEJA(3)
1213         sage: e0,e1,e2 = J.gens()
1214         sage: e0*e0
1215         e0
1216         sage: e0*e1
1217         0
1218         sage: e0*e2
1219         0
1220         sage: e1*e1
1221         e1
1222         sage: e1*e2
1223         0
1224         sage: e2*e2
1225         e2
1226
1227     """
1228     @staticmethod
1229     def __classcall_private__(cls, n, field=QQ):
1230         # The FiniteDimensionalAlgebra constructor takes a list of
1231         # matrices, the ith representing right multiplication by the ith
1232         # basis element in the vector space. So if e_1 = (1,0,0), then
1233         # right (Hadamard) multiplication of x by e_1 picks out the first
1234         # component of x; and likewise for the ith basis element e_i.
1235         Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i))
1236                for i in xrange(n) ]
1237
1238         fdeja = super(RealCartesianProductEJA, cls)
1239         return fdeja.__classcall_private__(cls, field, Qs, rank=n)
1240
1241     def inner_product(self, x, y):
1242         return _usual_ip(x,y)
1243
1244
1245 def random_eja():
1246     """
1247     Return a "random" finite-dimensional Euclidean Jordan Algebra.
1248
1249     ALGORITHM:
1250
1251     For now, we choose a random natural number ``n`` (greater than zero)
1252     and then give you back one of the following:
1253
1254       * The cartesian product of the rational numbers ``n`` times; this is
1255         ``QQ^n`` with the Hadamard product.
1256
1257       * The Jordan spin algebra on ``QQ^n``.
1258
1259       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1260         product.
1261
1262       * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1263         in the space of ``2n``-by-``2n`` real symmetric matrices.
1264
1265       * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1266         in the space of ``4n``-by-``4n`` real symmetric matrices.
1267
1268     Later this might be extended to return Cartesian products of the
1269     EJAs above.
1270
1271     TESTS::
1272
1273         sage: random_eja()
1274         Euclidean Jordan algebra of degree...
1275
1276     """
1277
1278     # The max_n component lets us choose different upper bounds on the
1279     # value "n" that gets passed to the constructor. This is needed
1280     # because e.g. R^{10} is reasonable to test, while the Hermitian
1281     # 10-by-10 quaternion matrices are not.
1282     (constructor, max_n) = choice([(RealCartesianProductEJA, 6),
1283                                    (JordanSpinEJA, 6),
1284                                    (RealSymmetricEJA, 5),
1285                                    (ComplexHermitianEJA, 4),
1286                                    (QuaternionHermitianEJA, 3)])
1287     n = ZZ.random_element(1, max_n)
1288     return constructor(n, field=QQ)
1289
1290
1291
1292 def _real_symmetric_basis(n, field=QQ):
1293     """
1294     Return a basis for the space of real symmetric n-by-n matrices.
1295     """
1296     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1297     # coordinates.
1298     S = []
1299     for i in xrange(n):
1300         for j in xrange(i+1):
1301             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1302             if i == j:
1303                 Sij = Eij
1304             else:
1305                 # Beware, orthogonal but not normalized!
1306                 Sij = Eij + Eij.transpose()
1307             S.append(Sij)
1308     return tuple(S)
1309
1310
1311 def _complex_hermitian_basis(n, field=QQ):
1312     """
1313     Returns a basis for the space of complex Hermitian n-by-n matrices.
1314
1315     TESTS::
1316
1317         sage: set_random_seed()
1318         sage: n = ZZ.random_element(1,5)
1319         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1320         True
1321
1322     """
1323     F = QuadraticField(-1, 'I')
1324     I = F.gen()
1325
1326     # This is like the symmetric case, but we need to be careful:
1327     #
1328     #   * We want conjugate-symmetry, not just symmetry.
1329     #   * The diagonal will (as a result) be real.
1330     #
1331     S = []
1332     for i in xrange(n):
1333         for j in xrange(i+1):
1334             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1335             if i == j:
1336                 Sij = _embed_complex_matrix(Eij)
1337                 S.append(Sij)
1338             else:
1339                 # Beware, orthogonal but not normalized! The second one
1340                 # has a minus because it's conjugated.
1341                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1342                 S.append(Sij_real)
1343                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1344                 S.append(Sij_imag)
1345     return tuple(S)
1346
1347
1348 def _quaternion_hermitian_basis(n, field=QQ):
1349     """
1350     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1351
1352     TESTS::
1353
1354         sage: set_random_seed()
1355         sage: n = ZZ.random_element(1,5)
1356         sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1357         True
1358
1359     """
1360     Q = QuaternionAlgebra(QQ,-1,-1)
1361     I,J,K = Q.gens()
1362
1363     # This is like the symmetric case, but we need to be careful:
1364     #
1365     #   * We want conjugate-symmetry, not just symmetry.
1366     #   * The diagonal will (as a result) be real.
1367     #
1368     S = []
1369     for i in xrange(n):
1370         for j in xrange(i+1):
1371             Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1372             if i == j:
1373                 Sij = _embed_quaternion_matrix(Eij)
1374                 S.append(Sij)
1375             else:
1376                 # Beware, orthogonal but not normalized! The second,
1377                 # third, and fourth ones have a minus because they're
1378                 # conjugated.
1379                 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1380                 S.append(Sij_real)
1381                 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1382                 S.append(Sij_I)
1383                 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1384                 S.append(Sij_J)
1385                 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1386                 S.append(Sij_K)
1387     return tuple(S)
1388
1389
1390 def _mat2vec(m):
1391         return vector(m.base_ring(), m.list())
1392
1393 def _vec2mat(v):
1394         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1395
1396 def _multiplication_table_from_matrix_basis(basis):
1397     """
1398     At least three of the five simple Euclidean Jordan algebras have the
1399     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1400     multiplication on the right is matrix multiplication. Given a basis
1401     for the underlying matrix space, this function returns a
1402     multiplication table (obtained by looping through the basis
1403     elements) for an algebra of those matrices. A reordered copy
1404     of the basis is also returned to work around the fact that
1405     the ``span()`` in this function will change the order of the basis
1406     from what we think it is, to... something else.
1407     """
1408     # In S^2, for example, we nominally have four coordinates even
1409     # though the space is of dimension three only. The vector space V
1410     # is supposed to hold the entire long vector, and the subspace W
1411     # of V will be spanned by the vectors that arise from symmetric
1412     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1413     field = basis[0].base_ring()
1414     dimension = basis[0].nrows()
1415
1416     V = VectorSpace(field, dimension**2)
1417     W = V.span( _mat2vec(s) for s in basis )
1418
1419     # Taking the span above reorders our basis (thanks, jerk!) so we
1420     # need to put our "matrix basis" in the same order as the
1421     # (reordered) vector basis.
1422     S = tuple( _vec2mat(b) for b in W.basis() )
1423
1424     Qs = []
1425     for s in S:
1426         # Brute force the multiplication-by-s matrix by looping
1427         # through all elements of the basis and doing the computation
1428         # to find out what the corresponding row should be. BEWARE:
1429         # these multiplication tables won't be symmetric! It therefore
1430         # becomes REALLY IMPORTANT that the underlying algebra
1431         # constructor uses ROW vectors and not COLUMN vectors. That's
1432         # why we're computing rows here and not columns.
1433         Q_rows = []
1434         for t in S:
1435             this_row = _mat2vec((s*t + t*s)/2)
1436             Q_rows.append(W.coordinates(this_row))
1437         Q = matrix(field, W.dimension(), Q_rows)
1438         Qs.append(Q)
1439
1440     return (Qs, S)
1441
1442
1443 def _embed_complex_matrix(M):
1444     """
1445     Embed the n-by-n complex matrix ``M`` into the space of real
1446     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1447     bi` to the block matrix ``[[a,b],[-b,a]]``.
1448
1449     EXAMPLES::
1450
1451         sage: F = QuadraticField(-1,'i')
1452         sage: x1 = F(4 - 2*i)
1453         sage: x2 = F(1 + 2*i)
1454         sage: x3 = F(-i)
1455         sage: x4 = F(6)
1456         sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1457         sage: _embed_complex_matrix(M)
1458         [ 4 -2| 1  2]
1459         [ 2  4|-2  1]
1460         [-----+-----]
1461         [ 0 -1| 6  0]
1462         [ 1  0| 0  6]
1463
1464     TESTS:
1465
1466     Embedding is a homomorphism (isomorphism, in fact)::
1467
1468         sage: set_random_seed()
1469         sage: n = ZZ.random_element(5)
1470         sage: F = QuadraticField(-1, 'i')
1471         sage: X = random_matrix(F, n)
1472         sage: Y = random_matrix(F, n)
1473         sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1474         sage: expected = _embed_complex_matrix(X*Y)
1475         sage: actual == expected
1476         True
1477
1478     """
1479     n = M.nrows()
1480     if M.ncols() != n:
1481         raise ValueError("the matrix 'M' must be square")
1482     field = M.base_ring()
1483     blocks = []
1484     for z in M.list():
1485         a = z.real()
1486         b = z.imag()
1487         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1488
1489     # We can drop the imaginaries here.
1490     return block_matrix(field.base_ring(), n, blocks)
1491
1492
1493 def _unembed_complex_matrix(M):
1494     """
1495     The inverse of _embed_complex_matrix().
1496
1497     EXAMPLES::
1498
1499         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1500         ....:                 [-2,  1,  -4,  3],
1501         ....:                 [ 9,  10, 11, 12],
1502         ....:                 [-10, 9, -12, 11] ])
1503         sage: _unembed_complex_matrix(A)
1504         [  2*i + 1   4*i + 3]
1505         [ 10*i + 9 12*i + 11]
1506
1507     TESTS:
1508
1509     Unembedding is the inverse of embedding::
1510
1511         sage: set_random_seed()
1512         sage: F = QuadraticField(-1, 'i')
1513         sage: M = random_matrix(F, 3)
1514         sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1515         True
1516
1517     """
1518     n = ZZ(M.nrows())
1519     if M.ncols() != n:
1520         raise ValueError("the matrix 'M' must be square")
1521     if not n.mod(2).is_zero():
1522         raise ValueError("the matrix 'M' must be a complex embedding")
1523
1524     F = QuadraticField(-1, 'i')
1525     i = F.gen()
1526
1527     # Go top-left to bottom-right (reading order), converting every
1528     # 2-by-2 block we see to a single complex element.
1529     elements = []
1530     for k in xrange(n/2):
1531         for j in xrange(n/2):
1532             submat = M[2*k:2*k+2,2*j:2*j+2]
1533             if submat[0,0] != submat[1,1]:
1534                 raise ValueError('bad on-diagonal submatrix')
1535             if submat[0,1] != -submat[1,0]:
1536                 raise ValueError('bad off-diagonal submatrix')
1537             z = submat[0,0] + submat[0,1]*i
1538             elements.append(z)
1539
1540     return matrix(F, n/2, elements)
1541
1542
1543 def _embed_quaternion_matrix(M):
1544     """
1545     Embed the n-by-n quaternion matrix ``M`` into the space of real
1546     matrices of size 4n-by-4n by first sending each quaternion entry
1547     `z = a + bi + cj + dk` to the block-complex matrix
1548     ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1549     a real matrix.
1550
1551     EXAMPLES::
1552
1553         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1554         sage: i,j,k = Q.gens()
1555         sage: x = 1 + 2*i + 3*j + 4*k
1556         sage: M = matrix(Q, 1, [[x]])
1557         sage: _embed_quaternion_matrix(M)
1558         [ 1  2  3  4]
1559         [-2  1 -4  3]
1560         [-3  4  1 -2]
1561         [-4 -3  2  1]
1562
1563     Embedding is a homomorphism (isomorphism, in fact)::
1564
1565         sage: set_random_seed()
1566         sage: n = ZZ.random_element(5)
1567         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1568         sage: X = random_matrix(Q, n)
1569         sage: Y = random_matrix(Q, n)
1570         sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1571         sage: expected = _embed_quaternion_matrix(X*Y)
1572         sage: actual == expected
1573         True
1574
1575     """
1576     quaternions = M.base_ring()
1577     n = M.nrows()
1578     if M.ncols() != n:
1579         raise ValueError("the matrix 'M' must be square")
1580
1581     F = QuadraticField(-1, 'i')
1582     i = F.gen()
1583
1584     blocks = []
1585     for z in M.list():
1586         t = z.coefficient_tuple()
1587         a = t[0]
1588         b = t[1]
1589         c = t[2]
1590         d = t[3]
1591         cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1592                                     [-c + d*i, a - b*i]])
1593         blocks.append(_embed_complex_matrix(cplx_matrix))
1594
1595     # We should have real entries by now, so use the realest field
1596     # we've got for the return value.
1597     return block_matrix(quaternions.base_ring(), n, blocks)
1598
1599
1600 def _unembed_quaternion_matrix(M):
1601     """
1602     The inverse of _embed_quaternion_matrix().
1603
1604     EXAMPLES::
1605
1606         sage: M = matrix(QQ, [[ 1,  2,  3,  4],
1607         ....:                 [-2,  1, -4,  3],
1608         ....:                 [-3,  4,  1, -2],
1609         ....:                 [-4, -3,  2,  1]])
1610         sage: _unembed_quaternion_matrix(M)
1611         [1 + 2*i + 3*j + 4*k]
1612
1613     TESTS:
1614
1615     Unembedding is the inverse of embedding::
1616
1617         sage: set_random_seed()
1618         sage: Q = QuaternionAlgebra(QQ, -1, -1)
1619         sage: M = random_matrix(Q, 3)
1620         sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1621         True
1622
1623     """
1624     n = ZZ(M.nrows())
1625     if M.ncols() != n:
1626         raise ValueError("the matrix 'M' must be square")
1627     if not n.mod(4).is_zero():
1628         raise ValueError("the matrix 'M' must be a complex embedding")
1629
1630     Q = QuaternionAlgebra(QQ,-1,-1)
1631     i,j,k = Q.gens()
1632
1633     # Go top-left to bottom-right (reading order), converting every
1634     # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1635     # quaternion block.
1636     elements = []
1637     for l in xrange(n/4):
1638         for m in xrange(n/4):
1639             submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1640             if submat[0,0] != submat[1,1].conjugate():
1641                 raise ValueError('bad on-diagonal submatrix')
1642             if submat[0,1] != -submat[1,0].conjugate():
1643                 raise ValueError('bad off-diagonal submatrix')
1644             z  = submat[0,0].real() + submat[0,0].imag()*i
1645             z += submat[0,1].real()*j + submat[0,1].imag()*k
1646             elements.append(z)
1647
1648     return matrix(Q, n/4, elements)
1649
1650
1651 # The usual inner product on R^n.
1652 def _usual_ip(x,y):
1653     return x.vector().inner_product(y.vector())
1654
1655 # The inner product used for the real symmetric simple EJA.
1656 # We keep it as a separate function because e.g. the complex
1657 # algebra uses the same inner product, except divided by 2.
1658 def _matrix_ip(X,Y):
1659     X_mat = X.natural_representation()
1660     Y_mat = Y.natural_representation()
1661     return (X_mat*Y_mat).trace()
1662
1663
1664 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
1665     """
1666     The rank-n simple EJA consisting of real symmetric n-by-n
1667     matrices, the usual symmetric Jordan product, and the trace inner
1668     product. It has dimension `(n^2 + n)/2` over the reals.
1669
1670     EXAMPLES::
1671
1672         sage: J = RealSymmetricEJA(2)
1673         sage: e0, e1, e2 = J.gens()
1674         sage: e0*e0
1675         e0
1676         sage: e1*e1
1677         e0 + e2
1678         sage: e2*e2
1679         e2
1680
1681     TESTS:
1682
1683     The degree of this algebra is `(n^2 + n) / 2`::
1684
1685         sage: set_random_seed()
1686         sage: n = ZZ.random_element(1,5)
1687         sage: J = RealSymmetricEJA(n)
1688         sage: J.degree() == (n^2 + n)/2
1689         True
1690
1691     The Jordan multiplication is what we think it is::
1692
1693         sage: set_random_seed()
1694         sage: n = ZZ.random_element(1,5)
1695         sage: J = RealSymmetricEJA(n)
1696         sage: x = J.random_element()
1697         sage: y = J.random_element()
1698         sage: actual = (x*y).natural_representation()
1699         sage: X = x.natural_representation()
1700         sage: Y = y.natural_representation()
1701         sage: expected = (X*Y + Y*X)/2
1702         sage: actual == expected
1703         True
1704         sage: J(expected) == x*y
1705         True
1706
1707     """
1708     @staticmethod
1709     def __classcall_private__(cls, n, field=QQ):
1710         S = _real_symmetric_basis(n, field=field)
1711         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1712
1713         fdeja = super(RealSymmetricEJA, cls)
1714         return fdeja.__classcall_private__(cls,
1715                                            field,
1716                                            Qs,
1717                                            rank=n,
1718                                            natural_basis=T)
1719
1720     def inner_product(self, x, y):
1721         return _matrix_ip(x,y)
1722
1723
1724 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1725     """
1726     The rank-n simple EJA consisting of complex Hermitian n-by-n
1727     matrices over the real numbers, the usual symmetric Jordan product,
1728     and the real-part-of-trace inner product. It has dimension `n^2` over
1729     the reals.
1730
1731     TESTS:
1732
1733     The degree of this algebra is `n^2`::
1734
1735         sage: set_random_seed()
1736         sage: n = ZZ.random_element(1,5)
1737         sage: J = ComplexHermitianEJA(n)
1738         sage: J.degree() == n^2
1739         True
1740
1741     The Jordan multiplication is what we think it is::
1742
1743         sage: set_random_seed()
1744         sage: n = ZZ.random_element(1,5)
1745         sage: J = ComplexHermitianEJA(n)
1746         sage: x = J.random_element()
1747         sage: y = J.random_element()
1748         sage: actual = (x*y).natural_representation()
1749         sage: X = x.natural_representation()
1750         sage: Y = y.natural_representation()
1751         sage: expected = (X*Y + Y*X)/2
1752         sage: actual == expected
1753         True
1754         sage: J(expected) == x*y
1755         True
1756
1757     """
1758     @staticmethod
1759     def __classcall_private__(cls, n, field=QQ):
1760         S = _complex_hermitian_basis(n)
1761         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1762
1763         fdeja = super(ComplexHermitianEJA, cls)
1764         return fdeja.__classcall_private__(cls,
1765                                            field,
1766                                            Qs,
1767                                            rank=n,
1768                                            natural_basis=T)
1769
1770     def inner_product(self, x, y):
1771         # Since a+bi on the diagonal is represented as
1772         #
1773         #   a + bi  = [  a  b  ]
1774         #             [ -b  a  ],
1775         #
1776         # we'll double-count the "a" entries if we take the trace of
1777         # the embedding.
1778         return _matrix_ip(x,y)/2
1779
1780
1781 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1782     """
1783     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1784     matrices, the usual symmetric Jordan product, and the
1785     real-part-of-trace inner product. It has dimension `2n^2 - n` over
1786     the reals.
1787
1788     TESTS:
1789
1790     The degree of this algebra is `n^2`::
1791
1792         sage: set_random_seed()
1793         sage: n = ZZ.random_element(1,5)
1794         sage: J = QuaternionHermitianEJA(n)
1795         sage: J.degree() == 2*(n^2) - n
1796         True
1797
1798     The Jordan multiplication is what we think it is::
1799
1800         sage: set_random_seed()
1801         sage: n = ZZ.random_element(1,5)
1802         sage: J = QuaternionHermitianEJA(n)
1803         sage: x = J.random_element()
1804         sage: y = J.random_element()
1805         sage: actual = (x*y).natural_representation()
1806         sage: X = x.natural_representation()
1807         sage: Y = y.natural_representation()
1808         sage: expected = (X*Y + Y*X)/2
1809         sage: actual == expected
1810         True
1811         sage: J(expected) == x*y
1812         True
1813
1814     """
1815     @staticmethod
1816     def __classcall_private__(cls, n, field=QQ):
1817         S = _quaternion_hermitian_basis(n)
1818         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1819
1820         fdeja = super(QuaternionHermitianEJA, cls)
1821         return fdeja.__classcall_private__(cls,
1822                                            field,
1823                                            Qs,
1824                                            rank=n,
1825                                            natural_basis=T)
1826
1827     def inner_product(self, x, y):
1828         # Since a+bi+cj+dk on the diagonal is represented as
1829         #
1830         #   a + bi +cj + dk = [  a  b  c  d]
1831         #                     [ -b  a -d  c]
1832         #                     [ -c  d  a -b]
1833         #                     [ -d -c  b  a],
1834         #
1835         # we'll quadruple-count the "a" entries if we take the trace of
1836         # the embedding.
1837         return _matrix_ip(x,y)/4
1838
1839
1840 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
1841     """
1842     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1843     with the usual inner product and jordan product ``x*y =
1844     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1845     the reals.
1846
1847     EXAMPLES:
1848
1849     This multiplication table can be verified by hand::
1850
1851         sage: J = JordanSpinEJA(4)
1852         sage: e0,e1,e2,e3 = J.gens()
1853         sage: e0*e0
1854         e0
1855         sage: e0*e1
1856         e1
1857         sage: e0*e2
1858         e2
1859         sage: e0*e3
1860         e3
1861         sage: e1*e2
1862         0
1863         sage: e1*e3
1864         0
1865         sage: e2*e3
1866         0
1867
1868     """
1869     @staticmethod
1870     def __classcall_private__(cls, n, field=QQ):
1871         Qs = []
1872         id_matrix = identity_matrix(field, n)
1873         for i in xrange(n):
1874             ei = id_matrix.column(i)
1875             Qi = zero_matrix(field, n)
1876             Qi.set_row(0, ei)
1877             Qi.set_column(0, ei)
1878             Qi += diagonal_matrix(n, [ei[0]]*n)
1879             # The addition of the diagonal matrix adds an extra ei[0] in the
1880             # upper-left corner of the matrix.
1881             Qi[0,0] = Qi[0,0] * ~field(2)
1882             Qs.append(Qi)
1883
1884         # The rank of the spin algebra is two, unless we're in a
1885         # one-dimensional ambient space (because the rank is bounded by
1886         # the ambient dimension).
1887         fdeja = super(JordanSpinEJA, cls)
1888         return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2))
1889
1890     def inner_product(self, x, y):
1891         return _usual_ip(x,y)