mjo/eja/euclidean_jordan_algebra.py

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