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