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