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