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