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