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
1176     # The max_n component lets us choose different upper bounds on the
1177     # value "n" that gets passed to the constructor. This is needed
1178     # because e.g. R^{10} is reasonable to test, while the Hermitian
1179     # 10-by-10 quaternion matrices are not.
1180     (constructor, max_n) = choice([(RealCartesianProductEJA, 6),
1181                                    (JordanSpinEJA, 6),
1182                                    (RealSymmetricEJA, 5),
1183                                    (ComplexHermitianEJA, 4),
1184                                    (QuaternionHermitianEJA, 3)])
1185     n = ZZ.random_element(1, max_n)
1186     return constructor(n, field=QQ)
1187
1188
1189
1190 def _real_symmetric_basis(n, field=QQ):
1191     """
1192     Return a basis for the space of real symmetric n-by-n matrices.
1193     """
1194     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1195     # coordinates.
1196     S = []
1197     for i in xrange(n):
1198         for j in xrange(i+1):
1199             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1200             if i == j:
1201                 Sij = Eij
1202             else:
1203                 # Beware, orthogonal but not normalized!
1204                 Sij = Eij + Eij.transpose()
1205             S.append(Sij)
1206     return tuple(S)
1207
1208
1209 def _complex_hermitian_basis(n, field=QQ):
1210     """
1211     Returns a basis for the space of complex Hermitian n-by-n matrices.
1212
1213     TESTS::
1214
1215         sage: set_random_seed()
1216         sage: n = ZZ.random_element(1,5)
1217         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1218         True
1219
1220     """
1221     F = QuadraticField(-1, 'I')
1222     I = F.gen()
1223
1224     # This is like the symmetric case, but we need to be careful:
1225     #
1226     #   * We want conjugate-symmetry, not just symmetry.
1227     #   * The diagonal will (as a result) be real.
1228     #
1229     S = []
1230     for i in xrange(n):
1231         for j in xrange(i+1):
1232             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1233             if i == j:
1234                 Sij = _embed_complex_matrix(Eij)
1235                 S.append(Sij)
1236             else:
1237                 # Beware, orthogonal but not normalized! The second one
1238                 # has a minus because it's conjugated.
1239                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1240                 S.append(Sij_real)
1241                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1242                 S.append(Sij_imag)
1243     return tuple(S)
1244
1245
1246 def _quaternion_hermitian_basis(n, field=QQ):
1247     """
1248     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1249
1250     TESTS::
1251
1252         sage: set_random_seed()
1253         sage: n = ZZ.random_element(1,5)
1254         sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1255         True
1256
1257     """
1258     Q = QuaternionAlgebra(QQ,-1,-1)
1259     I,J,K = Q.gens()
1260
1261     # This is like the symmetric case, but we need to be careful:
1262     #
1263     #   * We want conjugate-symmetry, not just symmetry.
1264     #   * The diagonal will (as a result) be real.
1265     #
1266     S = []
1267     for i in xrange(n):
1268         for j in xrange(i+1):
1269             Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1270             if i == j:
1271                 Sij = _embed_quaternion_matrix(Eij)
1272                 S.append(Sij)
1273             else:
1274                 # Beware, orthogonal but not normalized! The second,
1275                 # third, and fourth ones have a minus because they're
1276                 # conjugated.
1277                 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1278                 S.append(Sij_real)
1279                 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1280                 S.append(Sij_I)
1281                 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1282                 S.append(Sij_J)
1283                 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1284                 S.append(Sij_K)
1285     return tuple(S)
1286
1287
1288 def _mat2vec(m):
1289         return vector(m.base_ring(), m.list())
1290
1291 def _vec2mat(v):
1292         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1293
1294 def _multiplication_table_from_matrix_basis(basis):
1295     """
1296     At least three of the five simple Euclidean Jordan algebras have the
1297     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1298     multiplication on the right is matrix multiplication. Given a basis
1299     for the underlying matrix space, this function returns a
1300     multiplication table (obtained by looping through the basis
1301     elements) for an algebra of those matrices. A reordered copy
1302     of the basis is also returned to work around the fact that
1303     the ``span()`` in this function will change the order of the basis
1304     from what we think it is, to... something else.
1305     """
1306     # In S^2, for example, we nominally have four coordinates even
1307     # though the space is of dimension three only. The vector space V
1308     # is supposed to hold the entire long vector, and the subspace W
1309     # of V will be spanned by the vectors that arise from symmetric
1310     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1311     field = basis[0].base_ring()
1312     dimension = basis[0].nrows()
1313
1314     V = VectorSpace(field, dimension**2)
1315     W = V.span( _mat2vec(s) for s in basis )
1316
1317     # Taking the span above reorders our basis (thanks, jerk!) so we
1318     # need to put our "matrix basis" in the same order as the
1319     # (reordered) vector basis.
1320     S = tuple( _vec2mat(b) for b in W.basis() )
1321
1322     Qs = []
1323     for s in S:
1324         # Brute force the multiplication-by-s matrix by looping
1325         # through all elements of the basis and doing the computation
1326         # to find out what the corresponding row should be. BEWARE:
1327         # these multiplication tables won't be symmetric! It therefore
1328         # becomes REALLY IMPORTANT that the underlying algebra
1329         # constructor uses ROW vectors and not COLUMN vectors. That's
1330         # why we're computing rows here and not columns.
1331         Q_rows = []
1332         for t in S:
1333             this_row = _mat2vec((s*t + t*s)/2)
1334             Q_rows.append(W.coordinates(this_row))
1335         Q = matrix(field, W.dimension(), Q_rows)
1336         Qs.append(Q)
1337
1338     return (Qs, S)
1339
1340
1341 def _embed_complex_matrix(M):
1342     """
1343     Embed the n-by-n complex matrix ``M`` into the space of real
1344     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1345     bi` to the block matrix ``[[a,b],[-b,a]]``.
1346
1347     EXAMPLES::
1348
1349         sage: F = QuadraticField(-1,'i')
1350         sage: x1 = F(4 - 2*i)
1351         sage: x2 = F(1 + 2*i)
1352         sage: x3 = F(-i)
1353         sage: x4 = F(6)
1354         sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1355         sage: _embed_complex_matrix(M)
1356         [ 4 -2| 1  2]
1357         [ 2  4|-2  1]
1358         [-----+-----]
1359         [ 0 -1| 6  0]
1360         [ 1  0| 0  6]
1361
1362     TESTS:
1363
1364     Embedding is a homomorphism (isomorphism, in fact)::
1365
1366         sage: set_random_seed()
1367         sage: n = ZZ.random_element(5)
1368         sage: F = QuadraticField(-1, 'i')
1369         sage: X = random_matrix(F, n)
1370         sage: Y = random_matrix(F, n)
1371         sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1372         sage: expected = _embed_complex_matrix(X*Y)
1373         sage: actual == expected
1374         True
1375
1376     """
1377     n = M.nrows()
1378     if M.ncols() != n:
1379         raise ValueError("the matrix 'M' must be square")
1380     field = M.base_ring()
1381     blocks = []
1382     for z in M.list():
1383         a = z.real()
1384         b = z.imag()
1385         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1386
1387     # We can drop the imaginaries here.
1388     return block_matrix(field.base_ring(), n, blocks)
1389
1390
1391 def _unembed_complex_matrix(M):
1392     """
1393     The inverse of _embed_complex_matrix().
1394
1395     EXAMPLES::
1396
1397         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1398         ....:                 [-2,  1,  -4,  3],
1399         ....:                 [ 9,  10, 11, 12],
1400         ....:                 [-10, 9, -12, 11] ])
1401         sage: _unembed_complex_matrix(A)
1402         [  2*i + 1   4*i + 3]
1403         [ 10*i + 9 12*i + 11]
1404
1405     TESTS:
1406
1407     Unembedding is the inverse of embedding::
1408
1409         sage: set_random_seed()
1410         sage: F = QuadraticField(-1, 'i')
1411         sage: M = random_matrix(F, 3)
1412         sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1413         True
1414
1415     """
1416     n = ZZ(M.nrows())
1417     if M.ncols() != n:
1418         raise ValueError("the matrix 'M' must be square")
1419     if not n.mod(2).is_zero():
1420         raise ValueError("the matrix 'M' must be a complex embedding")
1421
1422     F = QuadraticField(-1, 'i')
1423     i = F.gen()
1424
1425     # Go top-left to bottom-right (reading order), converting every
1426     # 2-by-2 block we see to a single complex element.
1427     elements = []
1428     for k in xrange(n/2):
1429         for j in xrange(n/2):
1430             submat = M[2*k:2*k+2,2*j:2*j+2]
1431             if submat[0,0] != submat[1,1]:
1432                 raise ValueError('bad on-diagonal submatrix')
1433             if submat[0,1] != -submat[1,0]:
1434                 raise ValueError('bad off-diagonal submatrix')
1435             z = submat[0,0] + submat[0,1]*i
1436             elements.append(z)
1437
1438     return matrix(F, n/2, elements)
1439
1440
1441 def _embed_quaternion_matrix(M):
1442     """
1443     Embed the n-by-n quaternion matrix ``M`` into the space of real
1444     matrices of size 4n-by-4n by first sending each quaternion entry
1445     `z = a + bi + cj + dk` to the block-complex matrix
1446     ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1447     a real matrix.
1448
1449     EXAMPLES::
1450
1451         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1452         sage: i,j,k = Q.gens()
1453         sage: x = 1 + 2*i + 3*j + 4*k
1454         sage: M = matrix(Q, 1, [[x]])
1455         sage: _embed_quaternion_matrix(M)
1456         [ 1  2  3  4]
1457         [-2  1 -4  3]
1458         [-3  4  1 -2]
1459         [-4 -3  2  1]
1460
1461     Embedding is a homomorphism (isomorphism, in fact)::
1462
1463         sage: set_random_seed()
1464         sage: n = ZZ.random_element(5)
1465         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1466         sage: X = random_matrix(Q, n)
1467         sage: Y = random_matrix(Q, n)
1468         sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1469         sage: expected = _embed_quaternion_matrix(X*Y)
1470         sage: actual == expected
1471         True
1472
1473     """
1474     quaternions = M.base_ring()
1475     n = M.nrows()
1476     if M.ncols() != n:
1477         raise ValueError("the matrix 'M' must be square")
1478
1479     F = QuadraticField(-1, 'i')
1480     i = F.gen()
1481
1482     blocks = []
1483     for z in M.list():
1484         t = z.coefficient_tuple()
1485         a = t[0]
1486         b = t[1]
1487         c = t[2]
1488         d = t[3]
1489         cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1490                                     [-c + d*i, a - b*i]])
1491         blocks.append(_embed_complex_matrix(cplx_matrix))
1492
1493     # We should have real entries by now, so use the realest field
1494     # we've got for the return value.
1495     return block_matrix(quaternions.base_ring(), n, blocks)
1496
1497
1498 def _unembed_quaternion_matrix(M):
1499     """
1500     The inverse of _embed_quaternion_matrix().
1501
1502     EXAMPLES::
1503
1504         sage: M = matrix(QQ, [[ 1,  2,  3,  4],
1505         ....:                 [-2,  1, -4,  3],
1506         ....:                 [-3,  4,  1, -2],
1507         ....:                 [-4, -3,  2,  1]])
1508         sage: _unembed_quaternion_matrix(M)
1509         [1 + 2*i + 3*j + 4*k]
1510
1511     TESTS:
1512
1513     Unembedding is the inverse of embedding::
1514
1515         sage: set_random_seed()
1516         sage: Q = QuaternionAlgebra(QQ, -1, -1)
1517         sage: M = random_matrix(Q, 3)
1518         sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1519         True
1520
1521     """
1522     n = ZZ(M.nrows())
1523     if M.ncols() != n:
1524         raise ValueError("the matrix 'M' must be square")
1525     if not n.mod(4).is_zero():
1526         raise ValueError("the matrix 'M' must be a complex embedding")
1527
1528     Q = QuaternionAlgebra(QQ,-1,-1)
1529     i,j,k = Q.gens()
1530
1531     # Go top-left to bottom-right (reading order), converting every
1532     # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1533     # quaternion block.
1534     elements = []
1535     for l in xrange(n/4):
1536         for m in xrange(n/4):
1537             submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1538             if submat[0,0] != submat[1,1].conjugate():
1539                 raise ValueError('bad on-diagonal submatrix')
1540             if submat[0,1] != -submat[1,0].conjugate():
1541                 raise ValueError('bad off-diagonal submatrix')
1542             z  = submat[0,0].real() + submat[0,0].imag()*i
1543             z += submat[0,1].real()*j + submat[0,1].imag()*k
1544             elements.append(z)
1545
1546     return matrix(Q, n/4, elements)
1547
1548
1549 # The usual inner product on R^n.
1550 def _usual_ip(x,y):
1551     return x.vector().inner_product(y.vector())
1552
1553 # The inner product used for the real symmetric simple EJA.
1554 # We keep it as a separate function because e.g. the complex
1555 # algebra uses the same inner product, except divided by 2.
1556 def _matrix_ip(X,Y):
1557     X_mat = X.natural_representation()
1558     Y_mat = Y.natural_representation()
1559     return (X_mat*Y_mat).trace()
1560
1561
1562 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
1563     """
1564     The rank-n simple EJA consisting of real symmetric n-by-n
1565     matrices, the usual symmetric Jordan product, and the trace inner
1566     product. It has dimension `(n^2 + n)/2` over the reals.
1567
1568     EXAMPLES::
1569
1570         sage: J = RealSymmetricEJA(2)
1571         sage: e0, e1, e2 = J.gens()
1572         sage: e0*e0
1573         e0
1574         sage: e1*e1
1575         e0 + e2
1576         sage: e2*e2
1577         e2
1578
1579     TESTS:
1580
1581     The degree of this algebra is `(n^2 + n) / 2`::
1582
1583         sage: set_random_seed()
1584         sage: n = ZZ.random_element(1,5)
1585         sage: J = RealSymmetricEJA(n)
1586         sage: J.degree() == (n^2 + n)/2
1587         True
1588
1589     The Jordan multiplication is what we think it is::
1590
1591         sage: set_random_seed()
1592         sage: n = ZZ.random_element(1,5)
1593         sage: J = RealSymmetricEJA(n)
1594         sage: x = J.random_element()
1595         sage: y = J.random_element()
1596         sage: actual = (x*y).natural_representation()
1597         sage: X = x.natural_representation()
1598         sage: Y = y.natural_representation()
1599         sage: expected = (X*Y + Y*X)/2
1600         sage: actual == expected
1601         True
1602         sage: J(expected) == x*y
1603         True
1604
1605     """
1606     @staticmethod
1607     def __classcall_private__(cls, n, field=QQ):
1608         S = _real_symmetric_basis(n, field=field)
1609         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1610
1611         fdeja = super(RealSymmetricEJA, cls)
1612         return fdeja.__classcall_private__(cls,
1613                                            field,
1614                                            Qs,
1615                                            rank=n,
1616                                            natural_basis=T)
1617
1618     def inner_product(self, x, y):
1619         return _matrix_ip(x,y)
1620
1621
1622 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1623     """
1624     The rank-n simple EJA consisting of complex Hermitian n-by-n
1625     matrices over the real numbers, the usual symmetric Jordan product,
1626     and the real-part-of-trace inner product. It has dimension `n^2` over
1627     the reals.
1628
1629     TESTS:
1630
1631     The degree of this algebra is `n^2`::
1632
1633         sage: set_random_seed()
1634         sage: n = ZZ.random_element(1,5)
1635         sage: J = ComplexHermitianEJA(n)
1636         sage: J.degree() == n^2
1637         True
1638
1639     The Jordan multiplication is what we think it is::
1640
1641         sage: set_random_seed()
1642         sage: n = ZZ.random_element(1,5)
1643         sage: J = ComplexHermitianEJA(n)
1644         sage: x = J.random_element()
1645         sage: y = J.random_element()
1646         sage: actual = (x*y).natural_representation()
1647         sage: X = x.natural_representation()
1648         sage: Y = y.natural_representation()
1649         sage: expected = (X*Y + Y*X)/2
1650         sage: actual == expected
1651         True
1652         sage: J(expected) == x*y
1653         True
1654
1655     """
1656     @staticmethod
1657     def __classcall_private__(cls, n, field=QQ):
1658         S = _complex_hermitian_basis(n)
1659         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1660
1661         fdeja = super(ComplexHermitianEJA, cls)
1662         return fdeja.__classcall_private__(cls,
1663                                            field,
1664                                            Qs,
1665                                            rank=n,
1666                                            natural_basis=T)
1667
1668     def inner_product(self, x, y):
1669         # Since a+bi on the diagonal is represented as
1670         #
1671         #   a + bi  = [  a  b  ]
1672         #             [ -b  a  ],
1673         #
1674         # we'll double-count the "a" entries if we take the trace of
1675         # the embedding.
1676         return _matrix_ip(x,y)/2
1677
1678
1679 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1680     """
1681     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1682     matrices, the usual symmetric Jordan product, and the
1683     real-part-of-trace inner product. It has dimension `2n^2 - n` over
1684     the reals.
1685
1686     TESTS:
1687
1688     The degree of this algebra is `n^2`::
1689
1690         sage: set_random_seed()
1691         sage: n = ZZ.random_element(1,5)
1692         sage: J = QuaternionHermitianEJA(n)
1693         sage: J.degree() == 2*(n^2) - n
1694         True
1695
1696     The Jordan multiplication is what we think it is::
1697
1698         sage: set_random_seed()
1699         sage: n = ZZ.random_element(1,5)
1700         sage: J = QuaternionHermitianEJA(n)
1701         sage: x = J.random_element()
1702         sage: y = J.random_element()
1703         sage: actual = (x*y).natural_representation()
1704         sage: X = x.natural_representation()
1705         sage: Y = y.natural_representation()
1706         sage: expected = (X*Y + Y*X)/2
1707         sage: actual == expected
1708         True
1709         sage: J(expected) == x*y
1710         True
1711
1712     """
1713     @staticmethod
1714     def __classcall_private__(cls, n, field=QQ):
1715         S = _quaternion_hermitian_basis(n)
1716         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1717
1718         fdeja = super(QuaternionHermitianEJA, cls)
1719         return fdeja.__classcall_private__(cls,
1720                                            field,
1721                                            Qs,
1722                                            rank=n,
1723                                            natural_basis=T)
1724
1725     def inner_product(self, x, y):
1726         # Since a+bi+cj+dk on the diagonal is represented as
1727         #
1728         #   a + bi +cj + dk = [  a  b  c  d]
1729         #                     [ -b  a -d  c]
1730         #                     [ -c  d  a -b]
1731         #                     [ -d -c  b  a],
1732         #
1733         # we'll quadruple-count the "a" entries if we take the trace of
1734         # the embedding.
1735         return _matrix_ip(x,y)/4
1736
1737
1738 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
1739     """
1740     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1741     with the usual inner product and jordan product ``x*y =
1742     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1743     the reals.
1744
1745     EXAMPLES:
1746
1747     This multiplication table can be verified by hand::
1748
1749         sage: J = JordanSpinEJA(4)
1750         sage: e0,e1,e2,e3 = J.gens()
1751         sage: e0*e0
1752         e0
1753         sage: e0*e1
1754         e1
1755         sage: e0*e2
1756         e2
1757         sage: e0*e3
1758         e3
1759         sage: e1*e2
1760         0
1761         sage: e1*e3
1762         0
1763         sage: e2*e3
1764         0
1765
1766     """
1767     @staticmethod
1768     def __classcall_private__(cls, n, field=QQ):
1769         Qs = []
1770         id_matrix = identity_matrix(field, n)
1771         for i in xrange(n):
1772             ei = id_matrix.column(i)
1773             Qi = zero_matrix(field, n)
1774             Qi.set_row(0, ei)
1775             Qi.set_column(0, ei)
1776             Qi += diagonal_matrix(n, [ei[0]]*n)
1777             # The addition of the diagonal matrix adds an extra ei[0] in the
1778             # upper-left corner of the matrix.
1779             Qi[0,0] = Qi[0,0] * ~field(2)
1780             Qs.append(Qi)
1781
1782         # The rank of the spin algebra is two, unless we're in a
1783         # one-dimensional ambient space (because the rank is bounded by
1784         # the ambient dimension).
1785         fdeja = super(JordanSpinEJA, cls)
1786         return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2))
1787
1788     def inner_product(self, x, y):
1789         return _usual_ip(x,y)