mjo/eja/euclidean_jordan_algebra.py

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