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