mjo/eja/euclidean_jordan_algebra.py

   1 """
   2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
   3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
   4 are used in optimization, and have some additional nice methods beyond
   5 what can be supported in a general Jordan Algebra.
   6 """
   7
   8 from sage.categories.magmatic_algebras import MagmaticAlgebras
   9 from sage.structure.element import is_Matrix
  10 from sage.structure.category_object import normalize_names
  11
  12 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra
  13 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element import FiniteDimensionalAlgebraElement
  14
  15 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
  16     @staticmethod
  17     def __classcall_private__(cls,
  18                               field,
  19                               mult_table,
  20                               names='e',
  21                               assume_associative=False,
  22                               category=None,
  23                               rank=None,
  24                               natural_basis=None):
  25         n = len(mult_table)
  26         mult_table = [b.base_extend(field) for b in mult_table]
  27         for b in mult_table:
  28             b.set_immutable()
  29             if not (is_Matrix(b) and b.dimensions() == (n, n)):
  30                 raise ValueError("input is not a multiplication table")
  31         mult_table = tuple(mult_table)
  32
  33         cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
  34         cat.or_subcategory(category)
  35         if assume_associative:
  36             cat = cat.Associative()
  37
  38         names = normalize_names(n, names)
  39
  40         fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls)
  41         return fda.__classcall__(cls,
  42                                  field,
  43                                  mult_table,
  44                                  assume_associative=assume_associative,
  45                                  names=names,
  46                                  category=cat,
  47                                  rank=rank,
  48                                  natural_basis=natural_basis)
  49
  50
  51     def __init__(self,
  52                  field,
  53                  mult_table,
  54                  names='e',
  55                  assume_associative=False,
  56                  category=None,
  57                  rank=None,
  58                  natural_basis=None):
  59         """
  60         EXAMPLES:
  61
  62         By definition, Jordan multiplication commutes::
  63
  64             sage: set_random_seed()
  65             sage: J = random_eja()
  66             sage: x = J.random_element()
  67             sage: y = J.random_element()
  68             sage: x*y == y*x
  69             True
  70
  71         """
  72         self._rank = rank
  73         self._natural_basis = natural_basis
  74         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
  75         fda.__init__(field,
  76                      mult_table,
  77                      names=names,
  78                      category=category)
  79
  80
  81     def _repr_(self):
  82         """
  83         Return a string representation of ``self``.
  84         """
  85         fmt = "Euclidean Jordan algebra of degree {} over {}"
  86         return fmt.format(self.degree(), self.base_ring())
  87
  88
  89     def inner_product(self, x, y):
  90         """
  91         The inner product associated with this Euclidean Jordan algebra.
  92
  93         Defaults to the trace inner product, but can be overridden by
  94         subclasses if they are sure that the necessary properties are
  95         satisfied.
  96
  97         EXAMPLES:
  98
  99         The inner product must satisfy its axiom for this algebra to truly
 100         be a Euclidean Jordan Algebra::
 101
 102             sage: set_random_seed()
 103             sage: J = random_eja()
 104             sage: x = J.random_element()
 105             sage: y = J.random_element()
 106             sage: z = J.random_element()
 107             sage: (x*y).inner_product(z) == y.inner_product(x*z)
 108             True
 109
 110         """
 111         if (not x in self) or (not y in self):
 112             raise TypeError("arguments must live in this algebra")
 113         return x.trace_inner_product(y)
 114
 115
 116     def natural_basis(self):
 117         """
 118         Return a more-natural representation of this algebra's basis.
 119
 120         Every finite-dimensional Euclidean Jordan Algebra is a direct
 121         sum of five simple algebras, four of which comprise Hermitian
 122         matrices. This method returns the original "natural" basis
 123         for our underlying vector space. (Typically, the natural basis
 124         is used to construct the multiplication table in the first place.)
 125
 126         Note that this will always return a matrix. The standard basis
 127         in `R^n` will be returned as `n`-by-`1` column matrices.
 128
 129         EXAMPLES::
 130
 131             sage: J = RealSymmetricEJA(2)
 132             sage: J.basis()
 133             Family (e0, e1, e2)
 134             sage: J.natural_basis()
 135             (
 136             [1 0]  [0 1]  [0 0]
 137             [0 0], [1 0], [0 1]
 138             )
 139
 140         ::
 141
 142             sage: J = JordanSpinEJA(2)
 143             sage: J.basis()
 144             Family (e0, e1)
 145             sage: J.natural_basis()
 146             (
 147             [1]  [0]
 148             [0], [1]
 149             )
 150
 151         """
 152         if self._natural_basis is None:
 153             return tuple( b.vector().column() for b in self.basis() )
 154         else:
 155             return self._natural_basis
 156
 157
 158     def rank(self):
 159         """
 160         Return the rank of this EJA.
 161         """
 162         if self._rank is None:
 163             raise ValueError("no rank specified at genesis")
 164         else:
 165             return self._rank
 166
 167
 168     class Element(FiniteDimensionalAlgebraElement):
 169         """
 170         An element of a Euclidean Jordan algebra.
 171         """
 172
 173         def __init__(self, A, elt=None):
 174             """
 175             EXAMPLES:
 176
 177             The identity in `S^n` is converted to the identity in the EJA::
 178
 179                 sage: J = RealSymmetricEJA(3)
 180                 sage: I = identity_matrix(QQ,3)
 181                 sage: J(I) == J.one()
 182                 True
 183
 184             This skew-symmetric matrix can't be represented in the EJA::
 185
 186                 sage: J = RealSymmetricEJA(3)
 187                 sage: A = matrix(QQ,3, lambda i,j: i-j)
 188                 sage: J(A)
 189                 Traceback (most recent call last):
 190                 ...
 191                 ArithmeticError: vector is not in free module
 192
 193             """
 194             # Goal: if we're given a matrix, and if it lives in our
 195             # parent algebra's "natural ambient space," convert it
 196             # into an algebra element.
 197             #
 198             # The catch is, we make a recursive call after converting
 199             # the given matrix into a vector that lives in the algebra.
 200             # This we need to try the parent class initializer first,
 201             # to avoid recursing forever if we're given something that
 202             # already fits into the algebra, but also happens to live
 203             # in the parent's "natural ambient space" (this happens with
 204             # vectors in R^n).
 205             try:
 206                 FiniteDimensionalAlgebraElement.__init__(self, A, elt)
 207             except ValueError:
 208                 natural_basis = A.natural_basis()
 209                 if elt in natural_basis[0].matrix_space():
 210                     # Thanks for nothing! Matrix spaces aren't vector
 211                     # spaces in Sage, so we have to figure out its
 212                     # natural-basis coordinates ourselves.
 213                     V = VectorSpace(elt.base_ring(), elt.nrows()**2)
 214                     W = V.span( _mat2vec(s) for s in natural_basis )
 215                     coords =  W.coordinates(_mat2vec(elt))
 216                     FiniteDimensionalAlgebraElement.__init__(self, A, coords)
 217
 218         def __pow__(self, n):
 219             """
 220             Return ``self`` raised to the power ``n``.
 221
 222             Jordan algebras are always power-associative; see for
 223             example Faraut and Koranyi, Proposition II.1.2 (ii).
 224
 225             .. WARNING:
 226
 227                 We have to override this because our superclass uses row vectors
 228                 instead of column vectors! We, on the other hand, assume column
 229                 vectors everywhere.
 230
 231             EXAMPLES::
 232
 233                 sage: set_random_seed()
 234                 sage: x = random_eja().random_element()
 235                 sage: x.operator_matrix()*x.vector() == (x^2).vector()
 236                 True
 237
 238             A few examples of power-associativity::
 239
 240                 sage: set_random_seed()
 241                 sage: x = random_eja().random_element()
 242                 sage: x*(x*x)*(x*x) == x^5
 243                 True
 244                 sage: (x*x)*(x*x*x) == x^5
 245                 True
 246
 247             We also know that powers operator-commute (Koecher, Chapter
 248             III, Corollary 1)::
 249
 250                 sage: set_random_seed()
 251                 sage: x = random_eja().random_element()
 252                 sage: m = ZZ.random_element(0,10)
 253                 sage: n = ZZ.random_element(0,10)
 254                 sage: Lxm = (x^m).operator_matrix()
 255                 sage: Lxn = (x^n).operator_matrix()
 256                 sage: Lxm*Lxn == Lxn*Lxm
 257                 True
 258
 259             """
 260             A = self.parent()
 261             if n == 0:
 262                 return A.one()
 263             elif n == 1:
 264                 return self
 265             else:
 266                 return A( (self.operator_matrix()**(n-1))*self.vector() )
 267
 268
 269         def characteristic_polynomial(self):
 270             """
 271             Return my characteristic polynomial (if I'm a regular
 272             element).
 273
 274             Eventually this should be implemented in terms of the parent
 275             algebra's characteristic polynomial that works for ALL
 276             elements.
 277             """
 278             if self.is_regular():
 279                 return self.minimal_polynomial()
 280             else:
 281                 raise NotImplementedError('irregular element')
 282
 283
 284         def inner_product(self, other):
 285             """
 286             Return the parent algebra's inner product of myself and ``other``.
 287
 288             EXAMPLES:
 289
 290             The inner product in the Jordan spin algebra is the usual
 291             inner product on `R^n` (this example only works because the
 292             basis for the Jordan algebra is the standard basis in `R^n`)::
 293
 294                 sage: J = JordanSpinEJA(3)
 295                 sage: x = vector(QQ,[1,2,3])
 296                 sage: y = vector(QQ,[4,5,6])
 297                 sage: x.inner_product(y)
 298                 32
 299                 sage: J(x).inner_product(J(y))
 300                 32
 301
 302             The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
 303             multiplication is the usual matrix multiplication in `S^n`,
 304             so the inner product of the identity matrix with itself
 305             should be the `n`::
 306
 307                 sage: J = RealSymmetricEJA(3)
 308                 sage: J.one().inner_product(J.one())
 309                 3
 310
 311             Likewise, the inner product on `C^n` is `<X,Y> =
 312             Re(trace(X*Y))`, where we must necessarily take the real
 313             part because the product of Hermitian matrices may not be
 314             Hermitian::
 315
 316                 sage: J = ComplexHermitianEJA(3)
 317                 sage: J.one().inner_product(J.one())
 318                 3
 319
 320             Ditto for the quaternions::
 321
 322                 sage: J = QuaternionHermitianEJA(3)
 323                 sage: J.one().inner_product(J.one())
 324                 3
 325
 326             TESTS:
 327
 328             Ensure that we can always compute an inner product, and that
 329             it gives us back a real number::
 330
 331                 sage: set_random_seed()
 332                 sage: J = random_eja()
 333                 sage: x = J.random_element()
 334                 sage: y = J.random_element()
 335                 sage: x.inner_product(y) in RR
 336                 True
 337
 338             """
 339             P = self.parent()
 340             if not other in P:
 341                 raise TypeError("'other' must live in the same algebra")
 342
 343             return P.inner_product(self, other)
 344
 345
 346         def operator_commutes_with(self, other):
 347             """
 348             Return whether or not this element operator-commutes
 349             with ``other``.
 350
 351             EXAMPLES:
 352
 353             The definition of a Jordan algebra says that any element
 354             operator-commutes with its square::
 355
 356                 sage: set_random_seed()
 357                 sage: x = random_eja().random_element()
 358                 sage: x.operator_commutes_with(x^2)
 359                 True
 360
 361             TESTS:
 362
 363             Test Lemma 1 from Chapter III of Koecher::
 364
 365                 sage: set_random_seed()
 366                 sage: J = random_eja()
 367                 sage: u = J.random_element()
 368                 sage: v = J.random_element()
 369                 sage: lhs = u.operator_commutes_with(u*v)
 370                 sage: rhs = v.operator_commutes_with(u^2)
 371                 sage: lhs == rhs
 372                 True
 373
 374             """
 375             if not other in self.parent():
 376                 raise TypeError("'other' must live in the same algebra")
 377
 378             A = self.operator_matrix()
 379             B = other.operator_matrix()
 380             return (A*B == B*A)
 381
 382
 383         def det(self):
 384             """
 385             Return my determinant, the product of my eigenvalues.
 386
 387             EXAMPLES::
 388
 389                 sage: J = JordanSpinEJA(2)
 390                 sage: e0,e1 = J.gens()
 391                 sage: x = e0 + e1
 392                 sage: x.det()
 393                 0
 394                 sage: J = JordanSpinEJA(3)
 395                 sage: e0,e1,e2 = J.gens()
 396                 sage: x = e0 + e1 + e2
 397                 sage: x.det()
 398                 -1
 399
 400             """
 401             cs = self.characteristic_polynomial().coefficients(sparse=False)
 402             r = len(cs) - 1
 403             if r >= 0:
 404                 return cs[0] * (-1)**r
 405             else:
 406                 raise ValueError('charpoly had no coefficients')
 407
 408
 409         def inverse(self):
 410             """
 411             Return the Jordan-multiplicative inverse of this element.
 412
 413             We can't use the superclass method because it relies on the
 414             algebra being associative.
 415
 416             EXAMPLES:
 417
 418             The inverse in the spin factor algebra is given in Alizadeh's
 419             Example 11.11::
 420
 421                 sage: set_random_seed()
 422                 sage: n = ZZ.random_element(1,10)
 423                 sage: J = JordanSpinEJA(n)
 424                 sage: x = J.random_element()
 425                 sage: while x.is_zero():
 426                 ....:     x = J.random_element()
 427                 sage: x_vec = x.vector()
 428                 sage: x0 = x_vec[0]
 429                 sage: x_bar = x_vec[1:]
 430                 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
 431                 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
 432                 sage: x_inverse = coeff*inv_vec
 433                 sage: x.inverse() == J(x_inverse)
 434                 True
 435
 436             TESTS:
 437
 438             The identity element is its own inverse::
 439
 440                 sage: set_random_seed()
 441                 sage: J = random_eja()
 442                 sage: J.one().inverse() == J.one()
 443                 True
 444
 445             If an element has an inverse, it acts like one. TODO: this
 446             can be a lot less ugly once ``is_invertible`` doesn't crash
 447             on irregular elements::
 448
 449                 sage: set_random_seed()
 450                 sage: J = random_eja()
 451                 sage: x = J.random_element()
 452                 sage: try:
 453                 ....:     x.inverse()*x == J.one()
 454                 ....: except:
 455                 ....:     True
 456                 True
 457
 458             """
 459             if self.parent().is_associative():
 460                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 461                 return elt.inverse()
 462
 463             # TODO: we can do better once the call to is_invertible()
 464             # doesn't crash on irregular elements.
 465             #if not self.is_invertible():
 466             #    raise ValueError('element is not invertible')
 467
 468             # We do this a little different than the usual recursive
 469             # call to a finite-dimensional algebra element, because we
 470             # wind up with an inverse that lives in the subalgebra and
 471             # we need information about the parent to convert it back.
 472             V = self.span_of_powers()
 473             assoc_subalg = self.subalgebra_generated_by()
 474             # Mis-design warning: the basis used for span_of_powers()
 475             # and subalgebra_generated_by() must be the same, and in
 476             # the same order!
 477             elt = assoc_subalg(V.coordinates(self.vector()))
 478
 479             # This will be in the subalgebra's coordinates...
 480             fda_elt = FiniteDimensionalAlgebraElement(assoc_subalg, elt)
 481             subalg_inverse = fda_elt.inverse()
 482
 483             # So we have to convert back...
 484             basis = [ self.parent(v) for v in V.basis() ]
 485             pairs = zip(subalg_inverse.vector(), basis)
 486             return self.parent().linear_combination(pairs)
 487
 488
 489         def is_invertible(self):
 490             """
 491             Return whether or not this element is invertible.
 492
 493             We can't use the superclass method because it relies on
 494             the algebra being associative.
 495             """
 496             return not self.det().is_zero()
 497
 498
 499         def is_nilpotent(self):
 500             """
 501             Return whether or not some power of this element is zero.
 502
 503             The superclass method won't work unless we're in an
 504             associative algebra, and we aren't. However, we generate
 505             an assocoative subalgebra and we're nilpotent there if and
 506             only if we're nilpotent here (probably).
 507
 508             TESTS:
 509
 510             The identity element is never nilpotent::
 511
 512                 sage: set_random_seed()
 513                 sage: random_eja().one().is_nilpotent()
 514                 False
 515
 516             The additive identity is always nilpotent::
 517
 518                 sage: set_random_seed()
 519                 sage: random_eja().zero().is_nilpotent()
 520                 True
 521
 522             """
 523             # The element we're going to call "is_nilpotent()" on.
 524             # Either myself, interpreted as an element of a finite-
 525             # dimensional algebra, or an element of an associative
 526             # subalgebra.
 527             elt = None
 528
 529             if self.parent().is_associative():
 530                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 531             else:
 532                 V = self.span_of_powers()
 533                 assoc_subalg = self.subalgebra_generated_by()
 534                 # Mis-design warning: the basis used for span_of_powers()
 535                 # and subalgebra_generated_by() must be the same, and in
 536                 # the same order!
 537                 elt = assoc_subalg(V.coordinates(self.vector()))
 538
 539             # Recursive call, but should work since elt lives in an
 540             # associative algebra.
 541             return elt.is_nilpotent()
 542
 543
 544         def is_regular(self):
 545             """
 546             Return whether or not this is a regular element.
 547
 548             EXAMPLES:
 549
 550             The identity element always has degree one, but any element
 551             linearly-independent from it is regular::
 552
 553                 sage: J = JordanSpinEJA(5)
 554                 sage: J.one().is_regular()
 555                 False
 556                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 557                 sage: for x in J.gens():
 558                 ....:     (J.one() + x).is_regular()
 559                 False
 560                 True
 561                 True
 562                 True
 563                 True
 564
 565             """
 566             return self.degree() == self.parent().rank()
 567
 568
 569         def degree(self):
 570             """
 571             Compute the degree of this element the straightforward way
 572             according to the definition; by appending powers to a list
 573             and figuring out its dimension (that is, whether or not
 574             they're linearly dependent).
 575
 576             EXAMPLES::
 577
 578                 sage: J = JordanSpinEJA(4)
 579                 sage: J.one().degree()
 580                 1
 581                 sage: e0,e1,e2,e3 = J.gens()
 582                 sage: (e0 - e1).degree()
 583                 2
 584
 585             In the spin factor algebra (of rank two), all elements that
 586             aren't multiples of the identity are regular::
 587
 588                 sage: set_random_seed()
 589                 sage: n = ZZ.random_element(1,10)
 590                 sage: J = JordanSpinEJA(n)
 591                 sage: x = J.random_element()
 592                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 593                 True
 594
 595             """
 596             return self.span_of_powers().dimension()
 597
 598
 599         def minimal_polynomial(self):
 600             """
 601             EXAMPLES::
 602
 603                 sage: set_random_seed()
 604                 sage: x = random_eja().random_element()
 605                 sage: x.degree() == x.minimal_polynomial().degree()
 606                 True
 607
 608             ::
 609
 610                 sage: set_random_seed()
 611                 sage: x = random_eja().random_element()
 612                 sage: x.degree() == x.minimal_polynomial().degree()
 613                 True
 614
 615             The minimal polynomial and the characteristic polynomial coincide
 616             and are known (see Alizadeh, Example 11.11) for all elements of
 617             the spin factor algebra that aren't scalar multiples of the
 618             identity::
 619
 620                 sage: set_random_seed()
 621                 sage: n = ZZ.random_element(2,10)
 622                 sage: J = JordanSpinEJA(n)
 623                 sage: y = J.random_element()
 624                 sage: while y == y.coefficient(0)*J.one():
 625                 ....:     y = J.random_element()
 626                 sage: y0 = y.vector()[0]
 627                 sage: y_bar = y.vector()[1:]
 628                 sage: actual = y.minimal_polynomial()
 629                 sage: x = SR.symbol('x', domain='real')
 630                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 631                 sage: bool(actual == expected)
 632                 True
 633
 634             """
 635             # The element we're going to call "minimal_polynomial()" on.
 636             # Either myself, interpreted as an element of a finite-
 637             # dimensional algebra, or an element of an associative
 638             # subalgebra.
 639             elt = None
 640
 641             if self.parent().is_associative():
 642                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 643             else:
 644                 V = self.span_of_powers()
 645                 assoc_subalg = self.subalgebra_generated_by()
 646                 # Mis-design warning: the basis used for span_of_powers()
 647                 # and subalgebra_generated_by() must be the same, and in
 648                 # the same order!
 649                 elt = assoc_subalg(V.coordinates(self.vector()))
 650
 651             # Recursive call, but should work since elt lives in an
 652             # associative algebra.
 653             return elt.minimal_polynomial()
 654
 655
 656         def natural_representation(self):
 657             """
 658             Return a more-natural representation of this element.
 659
 660             Every finite-dimensional Euclidean Jordan Algebra is a
 661             direct sum of five simple algebras, four of which comprise
 662             Hermitian matrices. This method returns the original
 663             "natural" representation of this element as a Hermitian
 664             matrix, if it has one. If not, you get the usual representation.
 665
 666             EXAMPLES::
 667
 668                 sage: J = ComplexHermitianEJA(3)
 669                 sage: J.one()
 670                 e0 + e5 + e8
 671                 sage: J.one().natural_representation()
 672                 [1 0 0 0 0 0]
 673                 [0 1 0 0 0 0]
 674                 [0 0 1 0 0 0]
 675                 [0 0 0 1 0 0]
 676                 [0 0 0 0 1 0]
 677                 [0 0 0 0 0 1]
 678
 679             ::
 680
 681                 sage: J = QuaternionHermitianEJA(3)
 682                 sage: J.one()
 683                 e0 + e9 + e14
 684                 sage: J.one().natural_representation()
 685                 [1 0 0 0 0 0 0 0 0 0 0 0]
 686                 [0 1 0 0 0 0 0 0 0 0 0 0]
 687                 [0 0 1 0 0 0 0 0 0 0 0 0]
 688                 [0 0 0 1 0 0 0 0 0 0 0 0]
 689                 [0 0 0 0 1 0 0 0 0 0 0 0]
 690                 [0 0 0 0 0 1 0 0 0 0 0 0]
 691                 [0 0 0 0 0 0 1 0 0 0 0 0]
 692                 [0 0 0 0 0 0 0 1 0 0 0 0]
 693                 [0 0 0 0 0 0 0 0 1 0 0 0]
 694                 [0 0 0 0 0 0 0 0 0 1 0 0]
 695                 [0 0 0 0 0 0 0 0 0 0 1 0]
 696                 [0 0 0 0 0 0 0 0 0 0 0 1]
 697
 698             """
 699             B = self.parent().natural_basis()
 700             W = B[0].matrix_space()
 701             return W.linear_combination(zip(self.vector(), B))
 702
 703
 704         def operator_matrix(self):
 705             """
 706             Return the matrix that represents left- (or right-)
 707             multiplication by this element in the parent algebra.
 708
 709             We have to override this because the superclass method
 710             returns a matrix that acts on row vectors (that is, on
 711             the right).
 712
 713             EXAMPLES:
 714
 715             Test the first polarization identity from my notes, Koecher Chapter
 716             III, or from Baes (2.3)::
 717
 718                 sage: set_random_seed()
 719                 sage: J = random_eja()
 720                 sage: x = J.random_element()
 721                 sage: y = J.random_element()
 722                 sage: Lx = x.operator_matrix()
 723                 sage: Ly = y.operator_matrix()
 724                 sage: Lxx = (x*x).operator_matrix()
 725                 sage: Lxy = (x*y).operator_matrix()
 726                 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
 727                 True
 728
 729             Test the second polarization identity from my notes or from
 730             Baes (2.4)::
 731
 732                 sage: set_random_seed()
 733                 sage: J = random_eja()
 734                 sage: x = J.random_element()
 735                 sage: y = J.random_element()
 736                 sage: z = J.random_element()
 737                 sage: Lx = x.operator_matrix()
 738                 sage: Ly = y.operator_matrix()
 739                 sage: Lz = z.operator_matrix()
 740                 sage: Lzy = (z*y).operator_matrix()
 741                 sage: Lxy = (x*y).operator_matrix()
 742                 sage: Lxz = (x*z).operator_matrix()
 743                 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
 744                 True
 745
 746             Test the third polarization identity from my notes or from
 747             Baes (2.5)::
 748
 749                 sage: set_random_seed()
 750                 sage: J = random_eja()
 751                 sage: u = J.random_element()
 752                 sage: y = J.random_element()
 753                 sage: z = J.random_element()
 754                 sage: Lu = u.operator_matrix()
 755                 sage: Ly = y.operator_matrix()
 756                 sage: Lz = z.operator_matrix()
 757                 sage: Lzy = (z*y).operator_matrix()
 758                 sage: Luy = (u*y).operator_matrix()
 759                 sage: Luz = (u*z).operator_matrix()
 760                 sage: Luyz = (u*(y*z)).operator_matrix()
 761                 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
 762                 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
 763                 sage: bool(lhs == rhs)
 764                 True
 765
 766             """
 767             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 768             return fda_elt.matrix().transpose()
 769
 770
 771         def quadratic_representation(self, other=None):
 772             """
 773             Return the quadratic representation of this element.
 774
 775             EXAMPLES:
 776
 777             The explicit form in the spin factor algebra is given by
 778             Alizadeh's Example 11.12::
 779
 780                 sage: set_random_seed()
 781                 sage: n = ZZ.random_element(1,10)
 782                 sage: J = JordanSpinEJA(n)
 783                 sage: x = J.random_element()
 784                 sage: x_vec = x.vector()
 785                 sage: x0 = x_vec[0]
 786                 sage: x_bar = x_vec[1:]
 787                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 788                 sage: B = 2*x0*x_bar.row()
 789                 sage: C = 2*x0*x_bar.column()
 790                 sage: D = identity_matrix(QQ, n-1)
 791                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 792                 sage: D = D + 2*x_bar.tensor_product(x_bar)
 793                 sage: Q = block_matrix(2,2,[A,B,C,D])
 794                 sage: Q == x.quadratic_representation()
 795                 True
 796
 797             Test all of the properties from Theorem 11.2 in Alizadeh::
 798
 799                 sage: set_random_seed()
 800                 sage: J = random_eja()
 801                 sage: x = J.random_element()
 802                 sage: y = J.random_element()
 803
 804             Property 1:
 805
 806                 sage: actual = x.quadratic_representation(y)
 807                 sage: expected = ( (x+y).quadratic_representation()
 808                 ....:              -x.quadratic_representation()
 809                 ....:              -y.quadratic_representation() ) / 2
 810                 sage: actual == expected
 811                 True
 812
 813             Property 2:
 814
 815                 sage: alpha = QQ.random_element()
 816                 sage: actual = (alpha*x).quadratic_representation()
 817                 sage: expected = (alpha^2)*x.quadratic_representation()
 818                 sage: actual == expected
 819                 True
 820
 821             Property 5:
 822
 823                 sage: Qy = y.quadratic_representation()
 824                 sage: actual = J(Qy*x.vector()).quadratic_representation()
 825                 sage: expected = Qy*x.quadratic_representation()*Qy
 826                 sage: actual == expected
 827                 True
 828
 829             Property 6:
 830
 831                 sage: k = ZZ.random_element(1,10)
 832                 sage: actual = (x^k).quadratic_representation()
 833                 sage: expected = (x.quadratic_representation())^k
 834                 sage: actual == expected
 835                 True
 836
 837             """
 838             if other is None:
 839                 other=self
 840             elif not other in self.parent():
 841                 raise TypeError("'other' must live in the same algebra")
 842
 843             L = self.operator_matrix()
 844             M = other.operator_matrix()
 845             return ( L*M + M*L - (self*other).operator_matrix() )
 846
 847
 848         def span_of_powers(self):
 849             """
 850             Return the vector space spanned by successive powers of
 851             this element.
 852             """
 853             # The dimension of the subalgebra can't be greater than
 854             # the big algebra, so just put everything into a list
 855             # and let span() get rid of the excess.
 856             V = self.vector().parent()
 857             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 858
 859
 860         def subalgebra_generated_by(self):
 861             """
 862             Return the associative subalgebra of the parent EJA generated
 863             by this element.
 864
 865             TESTS::
 866
 867                 sage: set_random_seed()
 868                 sage: x = random_eja().random_element()
 869                 sage: x.subalgebra_generated_by().is_associative()
 870                 True
 871
 872             Squaring in the subalgebra should be the same thing as
 873             squaring in the superalgebra::
 874
 875                 sage: set_random_seed()
 876                 sage: x = random_eja().random_element()
 877                 sage: u = x.subalgebra_generated_by().random_element()
 878                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
 879                 True
 880
 881             """
 882             # First get the subspace spanned by the powers of myself...
 883             V = self.span_of_powers()
 884             F = self.base_ring()
 885
 886             # Now figure out the entries of the right-multiplication
 887             # matrix for the successive basis elements b0, b1,... of
 888             # that subspace.
 889             mats = []
 890             for b_right in V.basis():
 891                 eja_b_right = self.parent()(b_right)
 892                 b_right_rows = []
 893                 # The first row of the right-multiplication matrix by
 894                 # b1 is what we get if we apply that matrix to b1. The
 895                 # second row of the right multiplication matrix by b1
 896                 # is what we get when we apply that matrix to b2...
 897                 #
 898                 # IMPORTANT: this assumes that all vectors are COLUMN
 899                 # vectors, unlike our superclass (which uses row vectors).
 900                 for b_left in V.basis():
 901                     eja_b_left = self.parent()(b_left)
 902                     # Multiply in the original EJA, but then get the
 903                     # coordinates from the subalgebra in terms of its
 904                     # basis.
 905                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 906                     b_right_rows.append(this_row)
 907                 b_right_matrix = matrix(F, b_right_rows)
 908                 mats.append(b_right_matrix)
 909
 910             # It's an algebra of polynomials in one element, and EJAs
 911             # are power-associative.
 912             #
 913             # TODO: choose generator names intelligently.
 914             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 915
 916
 917         def subalgebra_idempotent(self):
 918             """
 919             Find an idempotent in the associative subalgebra I generate
 920             using Proposition 2.3.5 in Baes.
 921
 922             TESTS::
 923
 924                 sage: set_random_seed()
 925                 sage: J = eja_rn(5)
 926                 sage: c = J.random_element().subalgebra_idempotent()
 927                 sage: c^2 == c
 928                 True
 929                 sage: J = JordanSpinEJA(5)
 930                 sage: c = J.random_element().subalgebra_idempotent()
 931                 sage: c^2 == c
 932                 True
 933
 934             """
 935             if self.is_nilpotent():
 936                 raise ValueError("this only works with non-nilpotent elements!")
 937
 938             V = self.span_of_powers()
 939             J = self.subalgebra_generated_by()
 940             # Mis-design warning: the basis used for span_of_powers()
 941             # and subalgebra_generated_by() must be the same, and in
 942             # the same order!
 943             u = J(V.coordinates(self.vector()))
 944
 945             # The image of the matrix of left-u^m-multiplication
 946             # will be minimal for some natural number s...
 947             s = 0
 948             minimal_dim = V.dimension()
 949             for i in xrange(1, V.dimension()):
 950                 this_dim = (u**i).operator_matrix().image().dimension()
 951                 if this_dim < minimal_dim:
 952                     minimal_dim = this_dim
 953                     s = i
 954
 955             # Now minimal_matrix should correspond to the smallest
 956             # non-zero subspace in Baes's (or really, Koecher's)
 957             # proposition.
 958             #
 959             # However, we need to restrict the matrix to work on the
 960             # subspace... or do we? Can't we just solve, knowing that
 961             # A(c) = u^(s+1) should have a solution in the big space,
 962             # too?
 963             #
 964             # Beware, solve_right() means that we're using COLUMN vectors.
 965             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 966             u_next = u**(s+1)
 967             A = u_next.operator_matrix()
 968             c_coordinates = A.solve_right(u_next.vector())
 969
 970             # Now c_coordinates is the idempotent we want, but it's in
 971             # the coordinate system of the subalgebra.
 972             #
 973             # We need the basis for J, but as elements of the parent algebra.
 974             #
 975             basis = [self.parent(v) for v in V.basis()]
 976             return self.parent().linear_combination(zip(c_coordinates, basis))
 977
 978
 979         def trace(self):
 980             """
 981             Return my trace, the sum of my eigenvalues.
 982
 983             EXAMPLES::
 984
 985                 sage: J = JordanSpinEJA(3)
 986                 sage: e0,e1,e2 = J.gens()
 987                 sage: x = e0 + e1 + e2
 988                 sage: x.trace()
 989                 2
 990
 991             """
 992             cs = self.characteristic_polynomial().coefficients(sparse=False)
 993             if len(cs) >= 2:
 994                 return -1*cs[-2]
 995             else:
 996                 raise ValueError('charpoly had fewer than 2 coefficients')
 997
 998
 999         def trace_inner_product(self, other):
1000             """
1001             Return the trace inner product of myself and ``other``.
1002             """
1003             if not other in self.parent():
1004                 raise TypeError("'other' must live in the same algebra")
1005
1006             return (self*other).trace()
1007
1008
1009 def eja_rn(dimension, field=QQ):
1010     """
1011     Return the Euclidean Jordan Algebra corresponding to the set
1012     `R^n` under the Hadamard product.
1013
1014     EXAMPLES:
1015
1016     This multiplication table can be verified by hand::
1017
1018         sage: J = eja_rn(3)
1019         sage: e0,e1,e2 = J.gens()
1020         sage: e0*e0
1021         e0
1022         sage: e0*e1
1023         0
1024         sage: e0*e2
1025         0
1026         sage: e1*e1
1027         e1
1028         sage: e1*e2
1029         0
1030         sage: e2*e2
1031         e2
1032
1033     """
1034     # The FiniteDimensionalAlgebra constructor takes a list of
1035     # matrices, the ith representing right multiplication by the ith
1036     # basis element in the vector space. So if e_1 = (1,0,0), then
1037     # right (Hadamard) multiplication of x by e_1 picks out the first
1038     # component of x; and likewise for the ith basis element e_i.
1039     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
1040            for i in xrange(dimension) ]
1041
1042     return FiniteDimensionalEuclideanJordanAlgebra(field,
1043                                                    Qs,
1044                                                    rank=dimension)
1045
1046
1047
1048 def random_eja():
1049     """
1050     Return a "random" finite-dimensional Euclidean Jordan Algebra.
1051
1052     ALGORITHM:
1053
1054     For now, we choose a random natural number ``n`` (greater than zero)
1055     and then give you back one of the following:
1056
1057       * The cartesian product of the rational numbers ``n`` times; this is
1058         ``QQ^n`` with the Hadamard product.
1059
1060       * The Jordan spin algebra on ``QQ^n``.
1061
1062       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1063         product.
1064
1065       * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1066         in the space of ``2n``-by-``2n`` real symmetric matrices.
1067
1068       * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1069         in the space of ``4n``-by-``4n`` real symmetric matrices.
1070
1071     Later this might be extended to return Cartesian products of the
1072     EJAs above.
1073
1074     TESTS::
1075
1076         sage: random_eja()
1077         Euclidean Jordan algebra of degree...
1078
1079     """
1080     n = ZZ.random_element(1,5)
1081     constructor = choice([eja_rn,
1082                           JordanSpinEJA,
1083                           RealSymmetricEJA,
1084                           ComplexHermitianEJA,
1085                           QuaternionHermitianEJA])
1086     return constructor(n, field=QQ)
1087
1088
1089
1090 def _real_symmetric_basis(n, field=QQ):
1091     """
1092     Return a basis for the space of real symmetric n-by-n matrices.
1093     """
1094     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1095     # coordinates.
1096     S = []
1097     for i in xrange(n):
1098         for j in xrange(i+1):
1099             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1100             if i == j:
1101                 Sij = Eij
1102             else:
1103                 # Beware, orthogonal but not normalized!
1104                 Sij = Eij + Eij.transpose()
1105             S.append(Sij)
1106     return tuple(S)
1107
1108
1109 def _complex_hermitian_basis(n, field=QQ):
1110     """
1111     Returns a basis for the space of complex Hermitian n-by-n matrices.
1112
1113     TESTS::
1114
1115         sage: set_random_seed()
1116         sage: n = ZZ.random_element(1,5)
1117         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1118         True
1119
1120     """
1121     F = QuadraticField(-1, 'I')
1122     I = F.gen()
1123
1124     # This is like the symmetric case, but we need to be careful:
1125     #
1126     #   * We want conjugate-symmetry, not just symmetry.
1127     #   * The diagonal will (as a result) be real.
1128     #
1129     S = []
1130     for i in xrange(n):
1131         for j in xrange(i+1):
1132             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1133             if i == j:
1134                 Sij = _embed_complex_matrix(Eij)
1135                 S.append(Sij)
1136             else:
1137                 # Beware, orthogonal but not normalized! The second one
1138                 # has a minus because it's conjugated.
1139                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1140                 S.append(Sij_real)
1141                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1142                 S.append(Sij_imag)
1143     return tuple(S)
1144
1145
1146 def _quaternion_hermitian_basis(n, field=QQ):
1147     """
1148     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1149
1150     TESTS::
1151
1152         sage: set_random_seed()
1153         sage: n = ZZ.random_element(1,5)
1154         sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1155         True
1156
1157     """
1158     Q = QuaternionAlgebra(QQ,-1,-1)
1159     I,J,K = Q.gens()
1160
1161     # This is like the symmetric case, but we need to be careful:
1162     #
1163     #   * We want conjugate-symmetry, not just symmetry.
1164     #   * The diagonal will (as a result) be real.
1165     #
1166     S = []
1167     for i in xrange(n):
1168         for j in xrange(i+1):
1169             Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1170             if i == j:
1171                 Sij = _embed_quaternion_matrix(Eij)
1172                 S.append(Sij)
1173             else:
1174                 # Beware, orthogonal but not normalized! The second,
1175                 # third, and fourth ones have a minus because they're
1176                 # conjugated.
1177                 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1178                 S.append(Sij_real)
1179                 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1180                 S.append(Sij_I)
1181                 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1182                 S.append(Sij_J)
1183                 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1184                 S.append(Sij_K)
1185     return tuple(S)
1186
1187
1188 def _mat2vec(m):
1189         return vector(m.base_ring(), m.list())
1190
1191 def _vec2mat(v):
1192         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1193
1194 def _multiplication_table_from_matrix_basis(basis):
1195     """
1196     At least three of the five simple Euclidean Jordan algebras have the
1197     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1198     multiplication on the right is matrix multiplication. Given a basis
1199     for the underlying matrix space, this function returns a
1200     multiplication table (obtained by looping through the basis
1201     elements) for an algebra of those matrices. A reordered copy
1202     of the basis is also returned to work around the fact that
1203     the ``span()`` in this function will change the order of the basis
1204     from what we think it is, to... something else.
1205     """
1206     # In S^2, for example, we nominally have four coordinates even
1207     # though the space is of dimension three only. The vector space V
1208     # is supposed to hold the entire long vector, and the subspace W
1209     # of V will be spanned by the vectors that arise from symmetric
1210     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1211     field = basis[0].base_ring()
1212     dimension = basis[0].nrows()
1213
1214     V = VectorSpace(field, dimension**2)
1215     W = V.span( _mat2vec(s) for s in basis )
1216
1217     # Taking the span above reorders our basis (thanks, jerk!) so we
1218     # need to put our "matrix basis" in the same order as the
1219     # (reordered) vector basis.
1220     S = tuple( _vec2mat(b) for b in W.basis() )
1221
1222     Qs = []
1223     for s in S:
1224         # Brute force the multiplication-by-s matrix by looping
1225         # through all elements of the basis and doing the computation
1226         # to find out what the corresponding row should be. BEWARE:
1227         # these multiplication tables won't be symmetric! It therefore
1228         # becomes REALLY IMPORTANT that the underlying algebra
1229         # constructor uses ROW vectors and not COLUMN vectors. That's
1230         # why we're computing rows here and not columns.
1231         Q_rows = []
1232         for t in S:
1233             this_row = _mat2vec((s*t + t*s)/2)
1234             Q_rows.append(W.coordinates(this_row))
1235         Q = matrix(field, W.dimension(), Q_rows)
1236         Qs.append(Q)
1237
1238     return (Qs, S)
1239
1240
1241 def _embed_complex_matrix(M):
1242     """
1243     Embed the n-by-n complex matrix ``M`` into the space of real
1244     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1245     bi` to the block matrix ``[[a,b],[-b,a]]``.
1246
1247     EXAMPLES::
1248
1249         sage: F = QuadraticField(-1,'i')
1250         sage: x1 = F(4 - 2*i)
1251         sage: x2 = F(1 + 2*i)
1252         sage: x3 = F(-i)
1253         sage: x4 = F(6)
1254         sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1255         sage: _embed_complex_matrix(M)
1256         [ 4 -2| 1  2]
1257         [ 2  4|-2  1]
1258         [-----+-----]
1259         [ 0 -1| 6  0]
1260         [ 1  0| 0  6]
1261
1262     TESTS:
1263
1264     Embedding is a homomorphism (isomorphism, in fact)::
1265
1266         sage: set_random_seed()
1267         sage: n = ZZ.random_element(5)
1268         sage: F = QuadraticField(-1, 'i')
1269         sage: X = random_matrix(F, n)
1270         sage: Y = random_matrix(F, n)
1271         sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1272         sage: expected = _embed_complex_matrix(X*Y)
1273         sage: actual == expected
1274         True
1275
1276     """
1277     n = M.nrows()
1278     if M.ncols() != n:
1279         raise ValueError("the matrix 'M' must be square")
1280     field = M.base_ring()
1281     blocks = []
1282     for z in M.list():
1283         a = z.real()
1284         b = z.imag()
1285         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1286
1287     # We can drop the imaginaries here.
1288     return block_matrix(field.base_ring(), n, blocks)
1289
1290
1291 def _unembed_complex_matrix(M):
1292     """
1293     The inverse of _embed_complex_matrix().
1294
1295     EXAMPLES::
1296
1297         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1298         ....:                 [-2,  1,  -4,  3],
1299         ....:                 [ 9,  10, 11, 12],
1300         ....:                 [-10, 9, -12, 11] ])
1301         sage: _unembed_complex_matrix(A)
1302         [  2*i + 1   4*i + 3]
1303         [ 10*i + 9 12*i + 11]
1304
1305     TESTS:
1306
1307     Unembedding is the inverse of embedding::
1308
1309         sage: set_random_seed()
1310         sage: F = QuadraticField(-1, 'i')
1311         sage: M = random_matrix(F, 3)
1312         sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1313         True
1314
1315     """
1316     n = ZZ(M.nrows())
1317     if M.ncols() != n:
1318         raise ValueError("the matrix 'M' must be square")
1319     if not n.mod(2).is_zero():
1320         raise ValueError("the matrix 'M' must be a complex embedding")
1321
1322     F = QuadraticField(-1, 'i')
1323     i = F.gen()
1324
1325     # Go top-left to bottom-right (reading order), converting every
1326     # 2-by-2 block we see to a single complex element.
1327     elements = []
1328     for k in xrange(n/2):
1329         for j in xrange(n/2):
1330             submat = M[2*k:2*k+2,2*j:2*j+2]
1331             if submat[0,0] != submat[1,1]:
1332                 raise ValueError('bad on-diagonal submatrix')
1333             if submat[0,1] != -submat[1,0]:
1334                 raise ValueError('bad off-diagonal submatrix')
1335             z = submat[0,0] + submat[0,1]*i
1336             elements.append(z)
1337
1338     return matrix(F, n/2, elements)
1339
1340
1341 def _embed_quaternion_matrix(M):
1342     """
1343     Embed the n-by-n quaternion matrix ``M`` into the space of real
1344     matrices of size 4n-by-4n by first sending each quaternion entry
1345     `z = a + bi + cj + dk` to the block-complex matrix
1346     ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1347     a real matrix.
1348
1349     EXAMPLES::
1350
1351         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1352         sage: i,j,k = Q.gens()
1353         sage: x = 1 + 2*i + 3*j + 4*k
1354         sage: M = matrix(Q, 1, [[x]])
1355         sage: _embed_quaternion_matrix(M)
1356         [ 1  2  3  4]
1357         [-2  1 -4  3]
1358         [-3  4  1 -2]
1359         [-4 -3  2  1]
1360
1361     Embedding is a homomorphism (isomorphism, in fact)::
1362
1363         sage: set_random_seed()
1364         sage: n = ZZ.random_element(5)
1365         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1366         sage: X = random_matrix(Q, n)
1367         sage: Y = random_matrix(Q, n)
1368         sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1369         sage: expected = _embed_quaternion_matrix(X*Y)
1370         sage: actual == expected
1371         True
1372
1373     """
1374     quaternions = M.base_ring()
1375     n = M.nrows()
1376     if M.ncols() != n:
1377         raise ValueError("the matrix 'M' must be square")
1378
1379     F = QuadraticField(-1, 'i')
1380     i = F.gen()
1381
1382     blocks = []
1383     for z in M.list():
1384         t = z.coefficient_tuple()
1385         a = t[0]
1386         b = t[1]
1387         c = t[2]
1388         d = t[3]
1389         cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1390                                     [-c + d*i, a - b*i]])
1391         blocks.append(_embed_complex_matrix(cplx_matrix))
1392
1393     # We should have real entries by now, so use the realest field
1394     # we've got for the return value.
1395     return block_matrix(quaternions.base_ring(), n, blocks)
1396
1397
1398 def _unembed_quaternion_matrix(M):
1399     """
1400     The inverse of _embed_quaternion_matrix().
1401
1402     EXAMPLES::
1403
1404         sage: M = matrix(QQ, [[ 1,  2,  3,  4],
1405         ....:                 [-2,  1, -4,  3],
1406         ....:                 [-3,  4,  1, -2],
1407         ....:                 [-4, -3,  2,  1]])
1408         sage: _unembed_quaternion_matrix(M)
1409         [1 + 2*i + 3*j + 4*k]
1410
1411     TESTS:
1412
1413     Unembedding is the inverse of embedding::
1414
1415         sage: set_random_seed()
1416         sage: Q = QuaternionAlgebra(QQ, -1, -1)
1417         sage: M = random_matrix(Q, 3)
1418         sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1419         True
1420
1421     """
1422     n = ZZ(M.nrows())
1423     if M.ncols() != n:
1424         raise ValueError("the matrix 'M' must be square")
1425     if not n.mod(4).is_zero():
1426         raise ValueError("the matrix 'M' must be a complex embedding")
1427
1428     Q = QuaternionAlgebra(QQ,-1,-1)
1429     i,j,k = Q.gens()
1430
1431     # Go top-left to bottom-right (reading order), converting every
1432     # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1433     # quaternion block.
1434     elements = []
1435     for l in xrange(n/4):
1436         for m in xrange(n/4):
1437             submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1438             if submat[0,0] != submat[1,1].conjugate():
1439                 raise ValueError('bad on-diagonal submatrix')
1440             if submat[0,1] != -submat[1,0].conjugate():
1441                 raise ValueError('bad off-diagonal submatrix')
1442             z  = submat[0,0].real() + submat[0,0].imag()*i
1443             z += submat[0,1].real()*j + submat[0,1].imag()*k
1444             elements.append(z)
1445
1446     return matrix(Q, n/4, elements)
1447
1448
1449 # The usual inner product on R^n.
1450 def _usual_ip(x,y):
1451     return x.vector().inner_product(y.vector())
1452
1453 # The inner product used for the real symmetric simple EJA.
1454 # We keep it as a separate function because e.g. the complex
1455 # algebra uses the same inner product, except divided by 2.
1456 def _matrix_ip(X,Y):
1457     X_mat = X.natural_representation()
1458     Y_mat = Y.natural_representation()
1459     return (X_mat*Y_mat).trace()
1460
1461
1462 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
1463     """
1464     The rank-n simple EJA consisting of real symmetric n-by-n
1465     matrices, the usual symmetric Jordan product, and the trace inner
1466     product. It has dimension `(n^2 + n)/2` over the reals.
1467
1468     EXAMPLES::
1469
1470         sage: J = RealSymmetricEJA(2)
1471         sage: e0, e1, e2 = J.gens()
1472         sage: e0*e0
1473         e0
1474         sage: e1*e1
1475         e0 + e2
1476         sage: e2*e2
1477         e2
1478
1479     TESTS:
1480
1481     The degree of this algebra is `(n^2 + n) / 2`::
1482
1483         sage: set_random_seed()
1484         sage: n = ZZ.random_element(1,5)
1485         sage: J = RealSymmetricEJA(n)
1486         sage: J.degree() == (n^2 + n)/2
1487         True
1488
1489     The Jordan multiplication is what we think it is::
1490
1491         sage: set_random_seed()
1492         sage: n = ZZ.random_element(1,5)
1493         sage: J = RealSymmetricEJA(n)
1494         sage: x = J.random_element()
1495         sage: y = J.random_element()
1496         sage: actual = (x*y).natural_representation()
1497         sage: X = x.natural_representation()
1498         sage: Y = y.natural_representation()
1499         sage: expected = (X*Y + Y*X)/2
1500         sage: actual == expected
1501         True
1502         sage: J(expected) == x*y
1503         True
1504
1505     """
1506     @staticmethod
1507     def __classcall_private__(cls, n, field=QQ):
1508         S = _real_symmetric_basis(n, field=field)
1509         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1510
1511         fdeja = super(RealSymmetricEJA, cls)
1512         return fdeja.__classcall_private__(cls,
1513                                            field,
1514                                            Qs,
1515                                            rank=n,
1516                                            natural_basis=T)
1517
1518     def inner_product(self, x, y):
1519         return _matrix_ip(x,y)
1520
1521
1522 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1523     """
1524     The rank-n simple EJA consisting of complex Hermitian n-by-n
1525     matrices over the real numbers, the usual symmetric Jordan product,
1526     and the real-part-of-trace inner product. It has dimension `n^2` over
1527     the reals.
1528
1529     TESTS:
1530
1531     The degree of this algebra is `n^2`::
1532
1533         sage: set_random_seed()
1534         sage: n = ZZ.random_element(1,5)
1535         sage: J = ComplexHermitianEJA(n)
1536         sage: J.degree() == n^2
1537         True
1538
1539     The Jordan multiplication is what we think it is::
1540
1541         sage: set_random_seed()
1542         sage: n = ZZ.random_element(1,5)
1543         sage: J = ComplexHermitianEJA(n)
1544         sage: x = J.random_element()
1545         sage: y = J.random_element()
1546         sage: actual = (x*y).natural_representation()
1547         sage: X = x.natural_representation()
1548         sage: Y = y.natural_representation()
1549         sage: expected = (X*Y + Y*X)/2
1550         sage: actual == expected
1551         True
1552         sage: J(expected) == x*y
1553         True
1554
1555     """
1556     @staticmethod
1557     def __classcall_private__(cls, n, field=QQ):
1558         S = _complex_hermitian_basis(n)
1559         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1560
1561         fdeja = super(ComplexHermitianEJA, cls)
1562         return fdeja.__classcall_private__(cls,
1563                                            field,
1564                                            Qs,
1565                                            rank=n,
1566                                            natural_basis=T)
1567
1568     def inner_product(self, x, y):
1569         # Since a+bi on the diagonal is represented as
1570         #
1571         #   a + bi  = [  a  b  ]
1572         #             [ -b  a  ],
1573         #
1574         # we'll double-count the "a" entries if we take the trace of
1575         # the embedding.
1576         return _matrix_ip(x,y)/2
1577
1578
1579 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1580     """
1581     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1582     matrices, the usual symmetric Jordan product, and the
1583     real-part-of-trace inner product. It has dimension `2n^2 - n` over
1584     the reals.
1585
1586     TESTS:
1587
1588     The degree of this algebra is `n^2`::
1589
1590         sage: set_random_seed()
1591         sage: n = ZZ.random_element(1,5)
1592         sage: J = QuaternionHermitianEJA(n)
1593         sage: J.degree() == 2*(n^2) - n
1594         True
1595
1596     The Jordan multiplication is what we think it is::
1597
1598         sage: set_random_seed()
1599         sage: n = ZZ.random_element(1,5)
1600         sage: J = QuaternionHermitianEJA(n)
1601         sage: x = J.random_element()
1602         sage: y = J.random_element()
1603         sage: actual = (x*y).natural_representation()
1604         sage: X = x.natural_representation()
1605         sage: Y = y.natural_representation()
1606         sage: expected = (X*Y + Y*X)/2
1607         sage: actual == expected
1608         True
1609         sage: J(expected) == x*y
1610         True
1611
1612     """
1613     @staticmethod
1614     def __classcall_private__(cls, n, field=QQ):
1615         S = _quaternion_hermitian_basis(n)
1616         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1617
1618         fdeja = super(QuaternionHermitianEJA, cls)
1619         return fdeja.__classcall_private__(cls,
1620                                            field,
1621                                            Qs,
1622                                            rank=n,
1623                                            natural_basis=T)
1624
1625     def inner_product(self, x, y):
1626         # Since a+bi+cj+dk on the diagonal is represented as
1627         #
1628         #   a + bi +cj + dk = [  a  b  c  d]
1629         #                     [ -b  a -d  c]
1630         #                     [ -c  d  a -b]
1631         #                     [ -d -c  b  a],
1632         #
1633         # we'll quadruple-count the "a" entries if we take the trace of
1634         # the embedding.
1635         return _matrix_ip(x,y)/4
1636
1637
1638 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
1639     """
1640     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1641     with the usual inner product and jordan product ``x*y =
1642     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1643     the reals.
1644
1645     EXAMPLES:
1646
1647     This multiplication table can be verified by hand::
1648
1649         sage: J = JordanSpinEJA(4)
1650         sage: e0,e1,e2,e3 = J.gens()
1651         sage: e0*e0
1652         e0
1653         sage: e0*e1
1654         e1
1655         sage: e0*e2
1656         e2
1657         sage: e0*e3
1658         e3
1659         sage: e1*e2
1660         0
1661         sage: e1*e3
1662         0
1663         sage: e2*e3
1664         0
1665
1666     """
1667     @staticmethod
1668     def __classcall_private__(cls, n, field=QQ):
1669         Qs = []
1670         id_matrix = identity_matrix(field, n)
1671         for i in xrange(n):
1672             ei = id_matrix.column(i)
1673             Qi = zero_matrix(field, n)
1674             Qi.set_row(0, ei)
1675             Qi.set_column(0, ei)
1676             Qi += diagonal_matrix(n, [ei[0]]*n)
1677             # The addition of the diagonal matrix adds an extra ei[0] in the
1678             # upper-left corner of the matrix.
1679             Qi[0,0] = Qi[0,0] * ~field(2)
1680             Qs.append(Qi)
1681
1682         fdeja = super(JordanSpinEJA, cls)
1683         return fdeja.__classcall_private__(cls, field, Qs)
1684
1685     def rank(self):
1686         """
1687         Return the rank of this Jordan Spin Algebra.
1688
1689         The rank of the spin algebra is two, unless we're in a
1690         one-dimensional ambient space (because the rank is bounded by
1691         the ambient dimension).
1692         """
1693         return min(self.dimension(),2)
1694
1695     def inner_product(self, x, y):
1696         return _usual_ip(x,y)