mjo/eja/euclidean_jordan_algebra.py

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