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