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                               inner_product=None):
  26         n = len(mult_table)
  27         mult_table = [b.base_extend(field) for b in mult_table]
  28         for b in mult_table:
  29             b.set_immutable()
  30             if not (is_Matrix(b) and b.dimensions() == (n, n)):
  31                 raise ValueError("input is not a multiplication table")
  32         mult_table = tuple(mult_table)
  33
  34         cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
  35         cat.or_subcategory(category)
  36         if assume_associative:
  37             cat = cat.Associative()
  38
  39         names = normalize_names(n, names)
  40
  41         fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls)
  42         return fda.__classcall__(cls,
  43                                  field,
  44                                  mult_table,
  45                                  assume_associative=assume_associative,
  46                                  names=names,
  47                                  category=cat,
  48                                  rank=rank,
  49                                  natural_basis=natural_basis,
  50                                  inner_product=inner_product)
  51
  52
  53     def __init__(self, field,
  54                  mult_table,
  55                  names='e',
  56                  assume_associative=False,
  57                  category=None,
  58                  rank=None,
  59                  natural_basis=None,
  60                  inner_product=None):
  61         """
  62         EXAMPLES:
  63
  64         By definition, Jordan multiplication commutes::
  65
  66             sage: set_random_seed()
  67             sage: J = random_eja()
  68             sage: x = J.random_element()
  69             sage: y = J.random_element()
  70             sage: x*y == y*x
  71             True
  72
  73         """
  74         self._rank = rank
  75         self._natural_basis = natural_basis
  76         self._inner_product = inner_product
  77         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
  78         fda.__init__(field,
  79                      mult_table,
  80                      names=names,
  81                      category=category)
  82
  83
  84     def _repr_(self):
  85         """
  86         Return a string representation of ``self``.
  87         """
  88         fmt = "Euclidean Jordan algebra of degree {} over {}"
  89         return fmt.format(self.degree(), self.base_ring())
  90
  91
  92     def inner_product(self, x, y):
  93         """
  94         The inner product associated with this Euclidean Jordan algebra.
  95
  96         Will default to the trace inner product if nothing else.
  97         """
  98         if (not x in self) or (not y in self):
  99             raise TypeError("arguments must live in this algebra")
 100         if self._inner_product is None:
 101             return x.trace_inner_product(y)
 102         else:
 103             return self._inner_product(x,y)
 104
 105
 106     def natural_basis(self):
 107         """
 108         Return a more-natural representation of this algebra's basis.
 109
 110         Every finite-dimensional Euclidean Jordan Algebra is a direct
 111         sum of five simple algebras, four of which comprise Hermitian
 112         matrices. This method returns the original "natural" basis
 113         for our underlying vector space. (Typically, the natural basis
 114         is used to construct the multiplication table in the first place.)
 115
 116         Note that this will always return a matrix. The standard basis
 117         in `R^n` will be returned as `n`-by-`1` column matrices.
 118
 119         EXAMPLES::
 120
 121             sage: J = RealSymmetricSimpleEJA(2)
 122             sage: J.basis()
 123             Family (e0, e1, e2)
 124             sage: J.natural_basis()
 125             (
 126             [1 0]  [0 1]  [0 0]
 127             [0 0], [1 0], [0 1]
 128             )
 129
 130         ::
 131
 132             sage: J = JordanSpinSimpleEJA(2)
 133             sage: J.basis()
 134             Family (e0, e1)
 135             sage: J.natural_basis()
 136             (
 137             [1]  [0]
 138             [0], [1]
 139             )
 140
 141         """
 142         if self._natural_basis is None:
 143             return tuple( b.vector().column() for b in self.basis() )
 144         else:
 145             return self._natural_basis
 146
 147
 148     def rank(self):
 149         """
 150         Return the rank of this EJA.
 151         """
 152         if self._rank is None:
 153             raise ValueError("no rank specified at genesis")
 154         else:
 155             return self._rank
 156
 157
 158     class Element(FiniteDimensionalAlgebraElement):
 159         """
 160         An element of a Euclidean Jordan algebra.
 161         """
 162
 163         def __init__(self, A, elt=None):
 164             """
 165             EXAMPLES:
 166
 167             The identity in `S^n` is converted to the identity in the EJA::
 168
 169                 sage: J = RealSymmetricSimpleEJA(3)
 170                 sage: I = identity_matrix(QQ,3)
 171                 sage: J(I) == J.one()
 172                 True
 173
 174             This skew-symmetric matrix can't be represented in the EJA::
 175
 176                 sage: J = RealSymmetricSimpleEJA(3)
 177                 sage: A = matrix(QQ,3, lambda i,j: i-j)
 178                 sage: J(A)
 179                 Traceback (most recent call last):
 180                 ...
 181                 ArithmeticError: vector is not in free module
 182
 183             """
 184             # Goal: if we're given a matrix, and if it lives in our
 185             # parent algebra's "natural ambient space," convert it
 186             # into an algebra element.
 187             #
 188             # The catch is, we make a recursive call after converting
 189             # the given matrix into a vector that lives in the algebra.
 190             # This we need to try the parent class initializer first,
 191             # to avoid recursing forever if we're given something that
 192             # already fits into the algebra, but also happens to live
 193             # in the parent's "natural ambient space" (this happens with
 194             # vectors in R^n).
 195             try:
 196                 FiniteDimensionalAlgebraElement.__init__(self, A, elt)
 197             except ValueError:
 198                 natural_basis = A.natural_basis()
 199                 if elt in natural_basis[0].matrix_space():
 200                     # Thanks for nothing! Matrix spaces aren't vector
 201                     # spaces in Sage, so we have to figure out its
 202                     # natural-basis coordinates ourselves.
 203                     V = VectorSpace(elt.base_ring(), elt.nrows()**2)
 204                     W = V.span( _mat2vec(s) for s in natural_basis )
 205                     coords =  W.coordinates(_mat2vec(elt))
 206                     FiniteDimensionalAlgebraElement.__init__(self, A, coords)
 207
 208         def __pow__(self, n):
 209             """
 210             Return ``self`` raised to the power ``n``.
 211
 212             Jordan algebras are always power-associative; see for
 213             example Faraut and Koranyi, Proposition II.1.2 (ii).
 214
 215             .. WARNING:
 216
 217                 We have to override this because our superclass uses row vectors
 218                 instead of column vectors! We, on the other hand, assume column
 219                 vectors everywhere.
 220
 221             EXAMPLES::
 222
 223                 sage: set_random_seed()
 224                 sage: x = random_eja().random_element()
 225                 sage: x.operator_matrix()*x.vector() == (x^2).vector()
 226                 True
 227
 228             A few examples of power-associativity::
 229
 230                 sage: set_random_seed()
 231                 sage: x = random_eja().random_element()
 232                 sage: x*(x*x)*(x*x) == x^5
 233                 True
 234                 sage: (x*x)*(x*x*x) == x^5
 235                 True
 236
 237             We also know that powers operator-commute (Koecher, Chapter
 238             III, Corollary 1)::
 239
 240                 sage: set_random_seed()
 241                 sage: x = random_eja().random_element()
 242                 sage: m = ZZ.random_element(0,10)
 243                 sage: n = ZZ.random_element(0,10)
 244                 sage: Lxm = (x^m).operator_matrix()
 245                 sage: Lxn = (x^n).operator_matrix()
 246                 sage: Lxm*Lxn == Lxn*Lxm
 247                 True
 248
 249             """
 250             A = self.parent()
 251             if n == 0:
 252                 return A.one()
 253             elif n == 1:
 254                 return self
 255             else:
 256                 return A( (self.operator_matrix()**(n-1))*self.vector() )
 257
 258
 259         def characteristic_polynomial(self):
 260             """
 261             Return my characteristic polynomial (if I'm a regular
 262             element).
 263
 264             Eventually this should be implemented in terms of the parent
 265             algebra's characteristic polynomial that works for ALL
 266             elements.
 267             """
 268             if self.is_regular():
 269                 return self.minimal_polynomial()
 270             else:
 271                 raise NotImplementedError('irregular element')
 272
 273
 274         def inner_product(self, other):
 275             """
 276             Return the parent algebra's inner product of myself and ``other``.
 277
 278             EXAMPLES:
 279
 280             The inner product in the Jordan spin algebra is the usual
 281             inner product on `R^n` (this example only works because the
 282             basis for the Jordan algebra is the standard basis in `R^n`)::
 283
 284                 sage: J = JordanSpinSimpleEJA(3)
 285                 sage: x = vector(QQ,[1,2,3])
 286                 sage: y = vector(QQ,[4,5,6])
 287                 sage: x.inner_product(y)
 288                 32
 289                 sage: J(x).inner_product(J(y))
 290                 32
 291
 292             The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
 293             multiplication is the usual matrix multiplication in `S^n`,
 294             so the inner product of the identity matrix with itself
 295             should be the `n`::
 296
 297                 sage: J = RealSymmetricSimpleEJA(3)
 298                 sage: J.one().inner_product(J.one())
 299                 3
 300
 301             Likewise, the inner product on `C^n` is `<X,Y> =
 302             Re(trace(X*Y))`, where we must necessarily take the real
 303             part because the product of Hermitian matrices may not be
 304             Hermitian::
 305
 306                 sage: J = ComplexHermitianSimpleEJA(3)
 307                 sage: J.one().inner_product(J.one())
 308                 3
 309
 310             TESTS:
 311
 312             Ensure that we can always compute an inner product, and that
 313             it gives us back a real number::
 314
 315                 sage: set_random_seed()
 316                 sage: J = random_eja()
 317                 sage: x = J.random_element()
 318                 sage: y = J.random_element()
 319                 sage: x.inner_product(y) in RR
 320                 True
 321
 322             """
 323             P = self.parent()
 324             if not other in P:
 325                 raise TypeError("'other' must live in the same algebra")
 326
 327             return P.inner_product(self, other)
 328
 329
 330         def operator_commutes_with(self, other):
 331             """
 332             Return whether or not this element operator-commutes
 333             with ``other``.
 334
 335             EXAMPLES:
 336
 337             The definition of a Jordan algebra says that any element
 338             operator-commutes with its square::
 339
 340                 sage: set_random_seed()
 341                 sage: x = random_eja().random_element()
 342                 sage: x.operator_commutes_with(x^2)
 343                 True
 344
 345             TESTS:
 346
 347             Test Lemma 1 from Chapter III of Koecher::
 348
 349                 sage: set_random_seed()
 350                 sage: J = random_eja()
 351                 sage: u = J.random_element()
 352                 sage: v = J.random_element()
 353                 sage: lhs = u.operator_commutes_with(u*v)
 354                 sage: rhs = v.operator_commutes_with(u^2)
 355                 sage: lhs == rhs
 356                 True
 357
 358             """
 359             if not other in self.parent():
 360                 raise TypeError("'other' must live in the same algebra")
 361
 362             A = self.operator_matrix()
 363             B = other.operator_matrix()
 364             return (A*B == B*A)
 365
 366
 367         def det(self):
 368             """
 369             Return my determinant, the product of my eigenvalues.
 370
 371             EXAMPLES::
 372
 373                 sage: J = JordanSpinSimpleEJA(2)
 374                 sage: e0,e1 = J.gens()
 375                 sage: x = e0 + e1
 376                 sage: x.det()
 377                 0
 378                 sage: J = JordanSpinSimpleEJA(3)
 379                 sage: e0,e1,e2 = J.gens()
 380                 sage: x = e0 + e1 + e2
 381                 sage: x.det()
 382                 -1
 383
 384             """
 385             cs = self.characteristic_polynomial().coefficients(sparse=False)
 386             r = len(cs) - 1
 387             if r >= 0:
 388                 return cs[0] * (-1)**r
 389             else:
 390                 raise ValueError('charpoly had no coefficients')
 391
 392
 393         def inverse(self):
 394             """
 395             Return the Jordan-multiplicative inverse of this element.
 396
 397             We can't use the superclass method because it relies on the
 398             algebra being associative.
 399
 400             EXAMPLES:
 401
 402             The inverse in the spin factor algebra is given in Alizadeh's
 403             Example 11.11::
 404
 405                 sage: set_random_seed()
 406                 sage: n = ZZ.random_element(1,10)
 407                 sage: J = JordanSpinSimpleEJA(n)
 408                 sage: x = J.random_element()
 409                 sage: while x.is_zero():
 410                 ....:     x = J.random_element()
 411                 sage: x_vec = x.vector()
 412                 sage: x0 = x_vec[0]
 413                 sage: x_bar = x_vec[1:]
 414                 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
 415                 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
 416                 sage: x_inverse = coeff*inv_vec
 417                 sage: x.inverse() == J(x_inverse)
 418                 True
 419
 420             TESTS:
 421
 422             The identity element is its own inverse::
 423
 424                 sage: set_random_seed()
 425                 sage: J = random_eja()
 426                 sage: J.one().inverse() == J.one()
 427                 True
 428
 429             If an element has an inverse, it acts like one. TODO: this
 430             can be a lot less ugly once ``is_invertible`` doesn't crash
 431             on irregular elements::
 432
 433                 sage: set_random_seed()
 434                 sage: J = random_eja()
 435                 sage: x = J.random_element()
 436                 sage: try:
 437                 ....:     x.inverse()*x == J.one()
 438                 ....: except:
 439                 ....:     True
 440                 True
 441
 442             """
 443             if self.parent().is_associative():
 444                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 445                 return elt.inverse()
 446
 447             # TODO: we can do better once the call to is_invertible()
 448             # doesn't crash on irregular elements.
 449             #if not self.is_invertible():
 450             #    raise ValueError('element is not invertible')
 451
 452             # We do this a little different than the usual recursive
 453             # call to a finite-dimensional algebra element, because we
 454             # wind up with an inverse that lives in the subalgebra and
 455             # we need information about the parent to convert it back.
 456             V = self.span_of_powers()
 457             assoc_subalg = self.subalgebra_generated_by()
 458             # Mis-design warning: the basis used for span_of_powers()
 459             # and subalgebra_generated_by() must be the same, and in
 460             # the same order!
 461             elt = assoc_subalg(V.coordinates(self.vector()))
 462
 463             # This will be in the subalgebra's coordinates...
 464             fda_elt = FiniteDimensionalAlgebraElement(assoc_subalg, elt)
 465             subalg_inverse = fda_elt.inverse()
 466
 467             # So we have to convert back...
 468             basis = [ self.parent(v) for v in V.basis() ]
 469             pairs = zip(subalg_inverse.vector(), basis)
 470             return self.parent().linear_combination(pairs)
 471
 472
 473         def is_invertible(self):
 474             """
 475             Return whether or not this element is invertible.
 476
 477             We can't use the superclass method because it relies on
 478             the algebra being associative.
 479             """
 480             return not self.det().is_zero()
 481
 482
 483         def is_nilpotent(self):
 484             """
 485             Return whether or not some power of this element is zero.
 486
 487             The superclass method won't work unless we're in an
 488             associative algebra, and we aren't. However, we generate
 489             an assocoative subalgebra and we're nilpotent there if and
 490             only if we're nilpotent here (probably).
 491
 492             TESTS:
 493
 494             The identity element is never nilpotent::
 495
 496                 sage: set_random_seed()
 497                 sage: random_eja().one().is_nilpotent()
 498                 False
 499
 500             The additive identity is always nilpotent::
 501
 502                 sage: set_random_seed()
 503                 sage: random_eja().zero().is_nilpotent()
 504                 True
 505
 506             """
 507             # The element we're going to call "is_nilpotent()" on.
 508             # Either myself, interpreted as an element of a finite-
 509             # dimensional algebra, or an element of an associative
 510             # subalgebra.
 511             elt = None
 512
 513             if self.parent().is_associative():
 514                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 515             else:
 516                 V = self.span_of_powers()
 517                 assoc_subalg = self.subalgebra_generated_by()
 518                 # Mis-design warning: the basis used for span_of_powers()
 519                 # and subalgebra_generated_by() must be the same, and in
 520                 # the same order!
 521                 elt = assoc_subalg(V.coordinates(self.vector()))
 522
 523             # Recursive call, but should work since elt lives in an
 524             # associative algebra.
 525             return elt.is_nilpotent()
 526
 527
 528         def is_regular(self):
 529             """
 530             Return whether or not this is a regular element.
 531
 532             EXAMPLES:
 533
 534             The identity element always has degree one, but any element
 535             linearly-independent from it is regular::
 536
 537                 sage: J = JordanSpinSimpleEJA(5)
 538                 sage: J.one().is_regular()
 539                 False
 540                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 541                 sage: for x in J.gens():
 542                 ....:     (J.one() + x).is_regular()
 543                 False
 544                 True
 545                 True
 546                 True
 547                 True
 548
 549             """
 550             return self.degree() == self.parent().rank()
 551
 552
 553         def degree(self):
 554             """
 555             Compute the degree of this element the straightforward way
 556             according to the definition; by appending powers to a list
 557             and figuring out its dimension (that is, whether or not
 558             they're linearly dependent).
 559
 560             EXAMPLES::
 561
 562                 sage: J = JordanSpinSimpleEJA(4)
 563                 sage: J.one().degree()
 564                 1
 565                 sage: e0,e1,e2,e3 = J.gens()
 566                 sage: (e0 - e1).degree()
 567                 2
 568
 569             In the spin factor algebra (of rank two), all elements that
 570             aren't multiples of the identity are regular::
 571
 572                 sage: set_random_seed()
 573                 sage: n = ZZ.random_element(1,10)
 574                 sage: J = JordanSpinSimpleEJA(n)
 575                 sage: x = J.random_element()
 576                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 577                 True
 578
 579             """
 580             return self.span_of_powers().dimension()
 581
 582
 583         def minimal_polynomial(self):
 584             """
 585             EXAMPLES::
 586
 587                 sage: set_random_seed()
 588                 sage: x = random_eja().random_element()
 589                 sage: x.degree() == x.minimal_polynomial().degree()
 590                 True
 591
 592             ::
 593
 594                 sage: set_random_seed()
 595                 sage: x = random_eja().random_element()
 596                 sage: x.degree() == x.minimal_polynomial().degree()
 597                 True
 598
 599             The minimal polynomial and the characteristic polynomial coincide
 600             and are known (see Alizadeh, Example 11.11) for all elements of
 601             the spin factor algebra that aren't scalar multiples of the
 602             identity::
 603
 604                 sage: set_random_seed()
 605                 sage: n = ZZ.random_element(2,10)
 606                 sage: J = JordanSpinSimpleEJA(n)
 607                 sage: y = J.random_element()
 608                 sage: while y == y.coefficient(0)*J.one():
 609                 ....:     y = J.random_element()
 610                 sage: y0 = y.vector()[0]
 611                 sage: y_bar = y.vector()[1:]
 612                 sage: actual = y.minimal_polynomial()
 613                 sage: x = SR.symbol('x', domain='real')
 614                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 615                 sage: bool(actual == expected)
 616                 True
 617
 618             """
 619             # The element we're going to call "minimal_polynomial()" on.
 620             # Either myself, interpreted as an element of a finite-
 621             # dimensional algebra, or an element of an associative
 622             # subalgebra.
 623             elt = None
 624
 625             if self.parent().is_associative():
 626                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 627             else:
 628                 V = self.span_of_powers()
 629                 assoc_subalg = self.subalgebra_generated_by()
 630                 # Mis-design warning: the basis used for span_of_powers()
 631                 # and subalgebra_generated_by() must be the same, and in
 632                 # the same order!
 633                 elt = assoc_subalg(V.coordinates(self.vector()))
 634
 635             # Recursive call, but should work since elt lives in an
 636             # associative algebra.
 637             return elt.minimal_polynomial()
 638
 639
 640         def natural_representation(self):
 641             """
 642             Return a more-natural representation of this element.
 643
 644             Every finite-dimensional Euclidean Jordan Algebra is a
 645             direct sum of five simple algebras, four of which comprise
 646             Hermitian matrices. This method returns the original
 647             "natural" representation of this element as a Hermitian
 648             matrix, if it has one. If not, you get the usual representation.
 649
 650             EXAMPLES::
 651
 652                 sage: J = ComplexHermitianSimpleEJA(3)
 653                 sage: J.one()
 654                 e0 + e5 + e8
 655                 sage: J.one().natural_representation()
 656                 [1 0 0 0 0 0]
 657                 [0 1 0 0 0 0]
 658                 [0 0 1 0 0 0]
 659                 [0 0 0 1 0 0]
 660                 [0 0 0 0 1 0]
 661                 [0 0 0 0 0 1]
 662
 663             """
 664             B = self.parent().natural_basis()
 665             W = B[0].matrix_space()
 666             return W.linear_combination(zip(self.vector(), B))
 667
 668
 669         def operator_matrix(self):
 670             """
 671             Return the matrix that represents left- (or right-)
 672             multiplication by this element in the parent algebra.
 673
 674             We have to override this because the superclass method
 675             returns a matrix that acts on row vectors (that is, on
 676             the right).
 677
 678             EXAMPLES:
 679
 680             Test the first polarization identity from my notes, Koecher Chapter
 681             III, or from Baes (2.3)::
 682
 683                 sage: set_random_seed()
 684                 sage: J = random_eja()
 685                 sage: x = J.random_element()
 686                 sage: y = J.random_element()
 687                 sage: Lx = x.operator_matrix()
 688                 sage: Ly = y.operator_matrix()
 689                 sage: Lxx = (x*x).operator_matrix()
 690                 sage: Lxy = (x*y).operator_matrix()
 691                 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
 692                 True
 693
 694             Test the second polarization identity from my notes or from
 695             Baes (2.4)::
 696
 697                 sage: set_random_seed()
 698                 sage: J = random_eja()
 699                 sage: x = J.random_element()
 700                 sage: y = J.random_element()
 701                 sage: z = J.random_element()
 702                 sage: Lx = x.operator_matrix()
 703                 sage: Ly = y.operator_matrix()
 704                 sage: Lz = z.operator_matrix()
 705                 sage: Lzy = (z*y).operator_matrix()
 706                 sage: Lxy = (x*y).operator_matrix()
 707                 sage: Lxz = (x*z).operator_matrix()
 708                 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
 709                 True
 710
 711             Test the third polarization identity from my notes or from
 712             Baes (2.5)::
 713
 714                 sage: set_random_seed()
 715                 sage: J = random_eja()
 716                 sage: u = J.random_element()
 717                 sage: y = J.random_element()
 718                 sage: z = J.random_element()
 719                 sage: Lu = u.operator_matrix()
 720                 sage: Ly = y.operator_matrix()
 721                 sage: Lz = z.operator_matrix()
 722                 sage: Lzy = (z*y).operator_matrix()
 723                 sage: Luy = (u*y).operator_matrix()
 724                 sage: Luz = (u*z).operator_matrix()
 725                 sage: Luyz = (u*(y*z)).operator_matrix()
 726                 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
 727                 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
 728                 sage: bool(lhs == rhs)
 729                 True
 730
 731             """
 732             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 733             return fda_elt.matrix().transpose()
 734
 735
 736         def quadratic_representation(self, other=None):
 737             """
 738             Return the quadratic representation of this element.
 739
 740             EXAMPLES:
 741
 742             The explicit form in the spin factor algebra is given by
 743             Alizadeh's Example 11.12::
 744
 745                 sage: set_random_seed()
 746                 sage: n = ZZ.random_element(1,10)
 747                 sage: J = JordanSpinSimpleEJA(n)
 748                 sage: x = J.random_element()
 749                 sage: x_vec = x.vector()
 750                 sage: x0 = x_vec[0]
 751                 sage: x_bar = x_vec[1:]
 752                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 753                 sage: B = 2*x0*x_bar.row()
 754                 sage: C = 2*x0*x_bar.column()
 755                 sage: D = identity_matrix(QQ, n-1)
 756                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 757                 sage: D = D + 2*x_bar.tensor_product(x_bar)
 758                 sage: Q = block_matrix(2,2,[A,B,C,D])
 759                 sage: Q == x.quadratic_representation()
 760                 True
 761
 762             Test all of the properties from Theorem 11.2 in Alizadeh::
 763
 764                 sage: set_random_seed()
 765                 sage: J = random_eja()
 766                 sage: x = J.random_element()
 767                 sage: y = J.random_element()
 768
 769             Property 1:
 770
 771                 sage: actual = x.quadratic_representation(y)
 772                 sage: expected = ( (x+y).quadratic_representation()
 773                 ....:              -x.quadratic_representation()
 774                 ....:              -y.quadratic_representation() ) / 2
 775                 sage: actual == expected
 776                 True
 777
 778             Property 2:
 779
 780                 sage: alpha = QQ.random_element()
 781                 sage: actual = (alpha*x).quadratic_representation()
 782                 sage: expected = (alpha^2)*x.quadratic_representation()
 783                 sage: actual == expected
 784                 True
 785
 786             Property 5:
 787
 788                 sage: Qy = y.quadratic_representation()
 789                 sage: actual = J(Qy*x.vector()).quadratic_representation()
 790                 sage: expected = Qy*x.quadratic_representation()*Qy
 791                 sage: actual == expected
 792                 True
 793
 794             Property 6:
 795
 796                 sage: k = ZZ.random_element(1,10)
 797                 sage: actual = (x^k).quadratic_representation()
 798                 sage: expected = (x.quadratic_representation())^k
 799                 sage: actual == expected
 800                 True
 801
 802             """
 803             if other is None:
 804                 other=self
 805             elif not other in self.parent():
 806                 raise TypeError("'other' must live in the same algebra")
 807
 808             L = self.operator_matrix()
 809             M = other.operator_matrix()
 810             return ( L*M + M*L - (self*other).operator_matrix() )
 811
 812
 813         def span_of_powers(self):
 814             """
 815             Return the vector space spanned by successive powers of
 816             this element.
 817             """
 818             # The dimension of the subalgebra can't be greater than
 819             # the big algebra, so just put everything into a list
 820             # and let span() get rid of the excess.
 821             V = self.vector().parent()
 822             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 823
 824
 825         def subalgebra_generated_by(self):
 826             """
 827             Return the associative subalgebra of the parent EJA generated
 828             by this element.
 829
 830             TESTS::
 831
 832                 sage: set_random_seed()
 833                 sage: x = random_eja().random_element()
 834                 sage: x.subalgebra_generated_by().is_associative()
 835                 True
 836
 837             Squaring in the subalgebra should be the same thing as
 838             squaring in the superalgebra::
 839
 840                 sage: set_random_seed()
 841                 sage: x = random_eja().random_element()
 842                 sage: u = x.subalgebra_generated_by().random_element()
 843                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
 844                 True
 845
 846             """
 847             # First get the subspace spanned by the powers of myself...
 848             V = self.span_of_powers()
 849             F = self.base_ring()
 850
 851             # Now figure out the entries of the right-multiplication
 852             # matrix for the successive basis elements b0, b1,... of
 853             # that subspace.
 854             mats = []
 855             for b_right in V.basis():
 856                 eja_b_right = self.parent()(b_right)
 857                 b_right_rows = []
 858                 # The first row of the right-multiplication matrix by
 859                 # b1 is what we get if we apply that matrix to b1. The
 860                 # second row of the right multiplication matrix by b1
 861                 # is what we get when we apply that matrix to b2...
 862                 #
 863                 # IMPORTANT: this assumes that all vectors are COLUMN
 864                 # vectors, unlike our superclass (which uses row vectors).
 865                 for b_left in V.basis():
 866                     eja_b_left = self.parent()(b_left)
 867                     # Multiply in the original EJA, but then get the
 868                     # coordinates from the subalgebra in terms of its
 869                     # basis.
 870                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 871                     b_right_rows.append(this_row)
 872                 b_right_matrix = matrix(F, b_right_rows)
 873                 mats.append(b_right_matrix)
 874
 875             # It's an algebra of polynomials in one element, and EJAs
 876             # are power-associative.
 877             #
 878             # TODO: choose generator names intelligently.
 879             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 880
 881
 882         def subalgebra_idempotent(self):
 883             """
 884             Find an idempotent in the associative subalgebra I generate
 885             using Proposition 2.3.5 in Baes.
 886
 887             TESTS::
 888
 889                 sage: set_random_seed()
 890                 sage: J = eja_rn(5)
 891                 sage: c = J.random_element().subalgebra_idempotent()
 892                 sage: c^2 == c
 893                 True
 894                 sage: J = JordanSpinSimpleEJA(5)
 895                 sage: c = J.random_element().subalgebra_idempotent()
 896                 sage: c^2 == c
 897                 True
 898
 899             """
 900             if self.is_nilpotent():
 901                 raise ValueError("this only works with non-nilpotent elements!")
 902
 903             V = self.span_of_powers()
 904             J = self.subalgebra_generated_by()
 905             # Mis-design warning: the basis used for span_of_powers()
 906             # and subalgebra_generated_by() must be the same, and in
 907             # the same order!
 908             u = J(V.coordinates(self.vector()))
 909
 910             # The image of the matrix of left-u^m-multiplication
 911             # will be minimal for some natural number s...
 912             s = 0
 913             minimal_dim = V.dimension()
 914             for i in xrange(1, V.dimension()):
 915                 this_dim = (u**i).operator_matrix().image().dimension()
 916                 if this_dim < minimal_dim:
 917                     minimal_dim = this_dim
 918                     s = i
 919
 920             # Now minimal_matrix should correspond to the smallest
 921             # non-zero subspace in Baes's (or really, Koecher's)
 922             # proposition.
 923             #
 924             # However, we need to restrict the matrix to work on the
 925             # subspace... or do we? Can't we just solve, knowing that
 926             # A(c) = u^(s+1) should have a solution in the big space,
 927             # too?
 928             #
 929             # Beware, solve_right() means that we're using COLUMN vectors.
 930             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 931             u_next = u**(s+1)
 932             A = u_next.operator_matrix()
 933             c_coordinates = A.solve_right(u_next.vector())
 934
 935             # Now c_coordinates is the idempotent we want, but it's in
 936             # the coordinate system of the subalgebra.
 937             #
 938             # We need the basis for J, but as elements of the parent algebra.
 939             #
 940             basis = [self.parent(v) for v in V.basis()]
 941             return self.parent().linear_combination(zip(c_coordinates, basis))
 942
 943
 944         def trace(self):
 945             """
 946             Return my trace, the sum of my eigenvalues.
 947
 948             EXAMPLES::
 949
 950                 sage: J = JordanSpinSimpleEJA(3)
 951                 sage: e0,e1,e2 = J.gens()
 952                 sage: x = e0 + e1 + e2
 953                 sage: x.trace()
 954                 2
 955
 956             """
 957             cs = self.characteristic_polynomial().coefficients(sparse=False)
 958             if len(cs) >= 2:
 959                 return -1*cs[-2]
 960             else:
 961                 raise ValueError('charpoly had fewer than 2 coefficients')
 962
 963
 964         def trace_inner_product(self, other):
 965             """
 966             Return the trace inner product of myself and ``other``.
 967             """
 968             if not other in self.parent():
 969                 raise TypeError("'other' must live in the same algebra")
 970
 971             return (self*other).trace()
 972
 973
 974 def eja_rn(dimension, field=QQ):
 975     """
 976     Return the Euclidean Jordan Algebra corresponding to the set
 977     `R^n` under the Hadamard product.
 978
 979     EXAMPLES:
 980
 981     This multiplication table can be verified by hand::
 982
 983         sage: J = eja_rn(3)
 984         sage: e0,e1,e2 = J.gens()
 985         sage: e0*e0
 986         e0
 987         sage: e0*e1
 988         0
 989         sage: e0*e2
 990         0
 991         sage: e1*e1
 992         e1
 993         sage: e1*e2
 994         0
 995         sage: e2*e2
 996         e2
 997
 998     """
 999     # The FiniteDimensionalAlgebra constructor takes a list of
1000     # matrices, the ith representing right multiplication by the ith
1001     # basis element in the vector space. So if e_1 = (1,0,0), then
1002     # right (Hadamard) multiplication of x by e_1 picks out the first
1003     # component of x; and likewise for the ith basis element e_i.
1004     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
1005            for i in xrange(dimension) ]
1006
1007     return FiniteDimensionalEuclideanJordanAlgebra(field,
1008                                                    Qs,
1009                                                    rank=dimension,
1010                                                    inner_product=_usual_ip)
1011
1012
1013
1014 def random_eja():
1015     """
1016     Return a "random" finite-dimensional Euclidean Jordan Algebra.
1017
1018     ALGORITHM:
1019
1020     For now, we choose a random natural number ``n`` (greater than zero)
1021     and then give you back one of the following:
1022
1023       * The cartesian product of the rational numbers ``n`` times; this is
1024         ``QQ^n`` with the Hadamard product.
1025
1026       * The Jordan spin algebra on ``QQ^n``.
1027
1028       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1029         product.
1030
1031     Later this might be extended to return Cartesian products of the
1032     EJAs above.
1033
1034     TESTS::
1035
1036         sage: random_eja()
1037         Euclidean Jordan algebra of degree...
1038
1039     """
1040     n = ZZ.random_element(1,5)
1041     constructor = choice([eja_rn,
1042                           JordanSpinSimpleEJA,
1043                           RealSymmetricSimpleEJA,
1044                           ComplexHermitianSimpleEJA])
1045     return constructor(n, field=QQ)
1046
1047
1048
1049 def _real_symmetric_basis(n, field=QQ):
1050     """
1051     Return a basis for the space of real symmetric n-by-n matrices.
1052     """
1053     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1054     # coordinates.
1055     S = []
1056     for i in xrange(n):
1057         for j in xrange(i+1):
1058             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1059             if i == j:
1060                 Sij = Eij
1061             else:
1062                 # Beware, orthogonal but not normalized!
1063                 Sij = Eij + Eij.transpose()
1064             S.append(Sij)
1065     return tuple(S)
1066
1067
1068 def _complex_hermitian_basis(n, field=QQ):
1069     """
1070     Returns a basis for the space of complex Hermitian n-by-n matrices.
1071
1072     TESTS::
1073
1074         sage: set_random_seed()
1075         sage: n = ZZ.random_element(1,5)
1076         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1077         True
1078
1079     """
1080     F = QuadraticField(-1, 'I')
1081     I = F.gen()
1082
1083     # This is like the symmetric case, but we need to be careful:
1084     #
1085     #   * We want conjugate-symmetry, not just symmetry.
1086     #   * The diagonal will (as a result) be real.
1087     #
1088     S = []
1089     for i in xrange(n):
1090         for j in xrange(i+1):
1091             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1092             if i == j:
1093                 Sij = _embed_complex_matrix(Eij)
1094                 S.append(Sij)
1095             else:
1096                 # Beware, orthogonal but not normalized! The second one
1097                 # has a minus because it's conjugated.
1098                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1099                 S.append(Sij_real)
1100                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1101                 S.append(Sij_imag)
1102     return tuple(S)
1103
1104
1105 def _mat2vec(m):
1106         return vector(m.base_ring(), m.list())
1107
1108 def _vec2mat(v):
1109         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1110
1111 def _multiplication_table_from_matrix_basis(basis):
1112     """
1113     At least three of the five simple Euclidean Jordan algebras have the
1114     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1115     multiplication on the right is matrix multiplication. Given a basis
1116     for the underlying matrix space, this function returns a
1117     multiplication table (obtained by looping through the basis
1118     elements) for an algebra of those matrices. A reordered copy
1119     of the basis is also returned to work around the fact that
1120     the ``span()`` in this function will change the order of the basis
1121     from what we think it is, to... something else.
1122     """
1123     # In S^2, for example, we nominally have four coordinates even
1124     # though the space is of dimension three only. The vector space V
1125     # is supposed to hold the entire long vector, and the subspace W
1126     # of V will be spanned by the vectors that arise from symmetric
1127     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1128     field = basis[0].base_ring()
1129     dimension = basis[0].nrows()
1130
1131     V = VectorSpace(field, dimension**2)
1132     W = V.span( _mat2vec(s) for s in basis )
1133
1134     # Taking the span above reorders our basis (thanks, jerk!) so we
1135     # need to put our "matrix basis" in the same order as the
1136     # (reordered) vector basis.
1137     S = tuple( _vec2mat(b) for b in W.basis() )
1138
1139     Qs = []
1140     for s in S:
1141         # Brute force the multiplication-by-s matrix by looping
1142         # through all elements of the basis and doing the computation
1143         # to find out what the corresponding row should be. BEWARE:
1144         # these multiplication tables won't be symmetric! It therefore
1145         # becomes REALLY IMPORTANT that the underlying algebra
1146         # constructor uses ROW vectors and not COLUMN vectors. That's
1147         # why we're computing rows here and not columns.
1148         Q_rows = []
1149         for t in S:
1150             this_row = _mat2vec((s*t + t*s)/2)
1151             Q_rows.append(W.coordinates(this_row))
1152         Q = matrix(field, W.dimension(), Q_rows)
1153         Qs.append(Q)
1154
1155     return (Qs, S)
1156
1157
1158 def _embed_complex_matrix(M):
1159     """
1160     Embed the n-by-n complex matrix ``M`` into the space of real
1161     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1162     bi` to the block matrix ``[[a,b],[-b,a]]``.
1163
1164     EXAMPLES::
1165
1166         sage: F = QuadraticField(-1,'i')
1167         sage: x1 = F(4 - 2*i)
1168         sage: x2 = F(1 + 2*i)
1169         sage: x3 = F(-i)
1170         sage: x4 = F(6)
1171         sage: M = matrix(F,2,[x1,x2,x3,x4])
1172         sage: _embed_complex_matrix(M)
1173         [ 4  2| 1 -2]
1174         [-2  4| 2  1]
1175         [-----+-----]
1176         [ 0  1| 6  0]
1177         [-1  0| 0  6]
1178
1179     """
1180     n = M.nrows()
1181     if M.ncols() != n:
1182         raise ValueError("the matrix 'M' must be square")
1183     field = M.base_ring()
1184     blocks = []
1185     for z in M.list():
1186         a = z.real()
1187         b = z.imag()
1188         blocks.append(matrix(field, 2, [[a,-b],[b,a]]))
1189
1190     # We can drop the imaginaries here.
1191     return block_matrix(field.base_ring(), n, blocks)
1192
1193
1194 def _unembed_complex_matrix(M):
1195     """
1196     The inverse of _embed_complex_matrix().
1197
1198     EXAMPLES::
1199
1200         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1201         ....:                 [-2,  1,  -4,  3],
1202         ....:                 [ 9,  10, 11, 12],
1203         ....:                 [-10, 9, -12, 11] ])
1204         sage: _unembed_complex_matrix(A)
1205         [  -2*i + 1   -4*i + 3]
1206         [ -10*i + 9 -12*i + 11]
1207     """
1208     n = ZZ(M.nrows())
1209     if M.ncols() != n:
1210         raise ValueError("the matrix 'M' must be square")
1211     if not n.mod(2).is_zero():
1212         raise ValueError("the matrix 'M' must be a complex embedding")
1213
1214     F = QuadraticField(-1, 'i')
1215     i = F.gen()
1216
1217     # Go top-left to bottom-right (reading order), converting every
1218     # 2-by-2 block we see to a single complex element.
1219     elements = []
1220     for k in xrange(n/2):
1221         for j in xrange(n/2):
1222             submat = M[2*k:2*k+2,2*j:2*j+2]
1223             if submat[0,0] != submat[1,1]:
1224                 raise ValueError('bad real submatrix')
1225             if submat[0,1] != -submat[1,0]:
1226                 raise ValueError('bad imag submatrix')
1227             z = submat[0,0] + submat[1,0]*i
1228             elements.append(z)
1229
1230     return matrix(F, n/2, elements)
1231
1232 # The usual inner product on R^n.
1233 def _usual_ip(x,y):
1234     return x.vector().inner_product(y.vector())
1235
1236 # The inner product used for the real symmetric simple EJA.
1237 # We keep it as a separate function because e.g. the complex
1238 # algebra uses the same inner product, except divided by 2.
1239 def _matrix_ip(X,Y):
1240     X_mat = X.natural_representation()
1241     Y_mat = Y.natural_representation()
1242     return (X_mat*Y_mat).trace()
1243
1244
1245 def RealSymmetricSimpleEJA(n, field=QQ):
1246     """
1247     The rank-n simple EJA consisting of real symmetric n-by-n
1248     matrices, the usual symmetric Jordan product, and the trace inner
1249     product. It has dimension `(n^2 + n)/2` over the reals.
1250
1251     EXAMPLES::
1252
1253         sage: J = RealSymmetricSimpleEJA(2)
1254         sage: e0, e1, e2 = J.gens()
1255         sage: e0*e0
1256         e0
1257         sage: e1*e1
1258         e0 + e2
1259         sage: e2*e2
1260         e2
1261
1262     TESTS:
1263
1264     The degree of this algebra is `(n^2 + n) / 2`::
1265
1266         sage: set_random_seed()
1267         sage: n = ZZ.random_element(1,5)
1268         sage: J = RealSymmetricSimpleEJA(n)
1269         sage: J.degree() == (n^2 + n)/2
1270         True
1271
1272     """
1273     S = _real_symmetric_basis(n, field=field)
1274     (Qs, T) = _multiplication_table_from_matrix_basis(S)
1275
1276     return FiniteDimensionalEuclideanJordanAlgebra(field,
1277                                                    Qs,
1278                                                    rank=n,
1279                                                    natural_basis=T,
1280                                                    inner_product=_matrix_ip)
1281
1282
1283 def ComplexHermitianSimpleEJA(n, field=QQ):
1284     """
1285     The rank-n simple EJA consisting of complex Hermitian n-by-n
1286     matrices over the real numbers, the usual symmetric Jordan product,
1287     and the real-part-of-trace inner product. It has dimension `n^2` over
1288     the reals.
1289
1290     TESTS:
1291
1292     The degree of this algebra is `n^2`::
1293
1294         sage: set_random_seed()
1295         sage: n = ZZ.random_element(1,5)
1296         sage: J = ComplexHermitianSimpleEJA(n)
1297         sage: J.degree() == n^2
1298         True
1299
1300     """
1301     S = _complex_hermitian_basis(n)
1302     (Qs, T) = _multiplication_table_from_matrix_basis(S)
1303
1304     # Since a+bi on the diagonal is represented as
1305     #
1306     #   a + bi  = [  a  b  ]
1307     #             [ -b  a  ],
1308     #
1309     # we'll double-count the "a" entries if we take the trace of
1310     # the embedding.
1311     ip = lambda X,Y: _matrix_ip(X,Y)/2
1312
1313     return FiniteDimensionalEuclideanJordanAlgebra(field,
1314                                                    Qs,
1315                                                    rank=n,
1316                                                    natural_basis=T,
1317                                                    inner_product=ip)
1318
1319
1320 def QuaternionHermitianSimpleEJA(n):
1321     """
1322     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1323     matrices, the usual symmetric Jordan product, and the
1324     real-part-of-trace inner product. It has dimension `2n^2 - n` over
1325     the reals.
1326     """
1327     pass
1328
1329 def OctonionHermitianSimpleEJA(n):
1330     """
1331     This shit be crazy. It has dimension 27 over the reals.
1332     """
1333     n = 3
1334     pass
1335
1336 def JordanSpinSimpleEJA(n, field=QQ):
1337     """
1338     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1339     with the usual inner product and jordan product ``x*y =
1340     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1341     the reals.
1342
1343     EXAMPLES:
1344
1345     This multiplication table can be verified by hand::
1346
1347         sage: J = JordanSpinSimpleEJA(4)
1348         sage: e0,e1,e2,e3 = J.gens()
1349         sage: e0*e0
1350         e0
1351         sage: e0*e1
1352         e1
1353         sage: e0*e2
1354         e2
1355         sage: e0*e3
1356         e3
1357         sage: e1*e2
1358         0
1359         sage: e1*e3
1360         0
1361         sage: e2*e3
1362         0
1363
1364     In one dimension, this is the reals under multiplication::
1365
1366         sage: J1 = JordanSpinSimpleEJA(1)
1367         sage: J2 = eja_rn(1)
1368         sage: J1 == J2
1369         True
1370
1371     """
1372     Qs = []
1373     id_matrix = identity_matrix(field, n)
1374     for i in xrange(n):
1375         ei = id_matrix.column(i)
1376         Qi = zero_matrix(field, n)
1377         Qi.set_row(0, ei)
1378         Qi.set_column(0, ei)
1379         Qi += diagonal_matrix(n, [ei[0]]*n)
1380         # The addition of the diagonal matrix adds an extra ei[0] in the
1381         # upper-left corner of the matrix.
1382         Qi[0,0] = Qi[0,0] * ~field(2)
1383         Qs.append(Qi)
1384
1385     # The rank of the spin factor algebra is two, UNLESS we're in a
1386     # one-dimensional ambient space (the rank is bounded by the
1387     # ambient dimension).
1388     return FiniteDimensionalEuclideanJordanAlgebra(field,
1389                                                    Qs,
1390                                                    rank=min(n,2),
1391                                                    inner_product=_usual_ip)