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