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 = RealCartesianProductEJA(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 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
1010     """
1011     Return the Euclidean Jordan Algebra corresponding to the set
1012     `R^n` under the Hadamard product.
1013
1014     Note: this is nothing more than the Cartesian product of ``n``
1015     copies of the spin algebra. Once Cartesian product algebras
1016     are implemented, this can go.
1017
1018     EXAMPLES:
1019
1020     This multiplication table can be verified by hand::
1021
1022         sage: J = RealCartesianProductEJA(3)
1023         sage: e0,e1,e2 = J.gens()
1024         sage: e0*e0
1025         e0
1026         sage: e0*e1
1027         0
1028         sage: e0*e2
1029         0
1030         sage: e1*e1
1031         e1
1032         sage: e1*e2
1033         0
1034         sage: e2*e2
1035         e2
1036
1037     """
1038     @staticmethod
1039     def __classcall_private__(cls, n, field=QQ):
1040         # The FiniteDimensionalAlgebra constructor takes a list of
1041         # matrices, the ith representing right multiplication by the ith
1042         # basis element in the vector space. So if e_1 = (1,0,0), then
1043         # right (Hadamard) multiplication of x by e_1 picks out the first
1044         # component of x; and likewise for the ith basis element e_i.
1045         Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i))
1046                for i in xrange(n) ]
1047
1048         fdeja = super(RealCartesianProductEJA, cls)
1049         return fdeja.__classcall_private__(cls, field, Qs, rank=n)
1050
1051     def inner_product(self, x, y):
1052         return _usual_ip(x,y)
1053
1054
1055 def random_eja():
1056     """
1057     Return a "random" finite-dimensional Euclidean Jordan Algebra.
1058
1059     ALGORITHM:
1060
1061     For now, we choose a random natural number ``n`` (greater than zero)
1062     and then give you back one of the following:
1063
1064       * The cartesian product of the rational numbers ``n`` times; this is
1065         ``QQ^n`` with the Hadamard product.
1066
1067       * The Jordan spin algebra on ``QQ^n``.
1068
1069       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1070         product.
1071
1072       * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1073         in the space of ``2n``-by-``2n`` real symmetric matrices.
1074
1075       * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1076         in the space of ``4n``-by-``4n`` real symmetric matrices.
1077
1078     Later this might be extended to return Cartesian products of the
1079     EJAs above.
1080
1081     TESTS::
1082
1083         sage: random_eja()
1084         Euclidean Jordan algebra of degree...
1085
1086     """
1087     n = ZZ.random_element(1,5)
1088     constructor = choice([RealCartesianProductEJA,
1089                           JordanSpinEJA,
1090                           RealSymmetricEJA,
1091                           ComplexHermitianEJA,
1092                           QuaternionHermitianEJA])
1093     return constructor(n, field=QQ)
1094
1095
1096
1097 def _real_symmetric_basis(n, field=QQ):
1098     """
1099     Return a basis for the space of real symmetric n-by-n matrices.
1100     """
1101     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1102     # coordinates.
1103     S = []
1104     for i in xrange(n):
1105         for j in xrange(i+1):
1106             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1107             if i == j:
1108                 Sij = Eij
1109             else:
1110                 # Beware, orthogonal but not normalized!
1111                 Sij = Eij + Eij.transpose()
1112             S.append(Sij)
1113     return tuple(S)
1114
1115
1116 def _complex_hermitian_basis(n, field=QQ):
1117     """
1118     Returns a basis for the space of complex Hermitian n-by-n matrices.
1119
1120     TESTS::
1121
1122         sage: set_random_seed()
1123         sage: n = ZZ.random_element(1,5)
1124         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1125         True
1126
1127     """
1128     F = QuadraticField(-1, 'I')
1129     I = F.gen()
1130
1131     # This is like the symmetric case, but we need to be careful:
1132     #
1133     #   * We want conjugate-symmetry, not just symmetry.
1134     #   * The diagonal will (as a result) be real.
1135     #
1136     S = []
1137     for i in xrange(n):
1138         for j in xrange(i+1):
1139             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1140             if i == j:
1141                 Sij = _embed_complex_matrix(Eij)
1142                 S.append(Sij)
1143             else:
1144                 # Beware, orthogonal but not normalized! The second one
1145                 # has a minus because it's conjugated.
1146                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1147                 S.append(Sij_real)
1148                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1149                 S.append(Sij_imag)
1150     return tuple(S)
1151
1152
1153 def _quaternion_hermitian_basis(n, field=QQ):
1154     """
1155     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1156
1157     TESTS::
1158
1159         sage: set_random_seed()
1160         sage: n = ZZ.random_element(1,5)
1161         sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1162         True
1163
1164     """
1165     Q = QuaternionAlgebra(QQ,-1,-1)
1166     I,J,K = Q.gens()
1167
1168     # This is like the symmetric case, but we need to be careful:
1169     #
1170     #   * We want conjugate-symmetry, not just symmetry.
1171     #   * The diagonal will (as a result) be real.
1172     #
1173     S = []
1174     for i in xrange(n):
1175         for j in xrange(i+1):
1176             Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1177             if i == j:
1178                 Sij = _embed_quaternion_matrix(Eij)
1179                 S.append(Sij)
1180             else:
1181                 # Beware, orthogonal but not normalized! The second,
1182                 # third, and fourth ones have a minus because they're
1183                 # conjugated.
1184                 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1185                 S.append(Sij_real)
1186                 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1187                 S.append(Sij_I)
1188                 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1189                 S.append(Sij_J)
1190                 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1191                 S.append(Sij_K)
1192     return tuple(S)
1193
1194
1195 def _mat2vec(m):
1196         return vector(m.base_ring(), m.list())
1197
1198 def _vec2mat(v):
1199         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1200
1201 def _multiplication_table_from_matrix_basis(basis):
1202     """
1203     At least three of the five simple Euclidean Jordan algebras have the
1204     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1205     multiplication on the right is matrix multiplication. Given a basis
1206     for the underlying matrix space, this function returns a
1207     multiplication table (obtained by looping through the basis
1208     elements) for an algebra of those matrices. A reordered copy
1209     of the basis is also returned to work around the fact that
1210     the ``span()`` in this function will change the order of the basis
1211     from what we think it is, to... something else.
1212     """
1213     # In S^2, for example, we nominally have four coordinates even
1214     # though the space is of dimension three only. The vector space V
1215     # is supposed to hold the entire long vector, and the subspace W
1216     # of V will be spanned by the vectors that arise from symmetric
1217     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1218     field = basis[0].base_ring()
1219     dimension = basis[0].nrows()
1220
1221     V = VectorSpace(field, dimension**2)
1222     W = V.span( _mat2vec(s) for s in basis )
1223
1224     # Taking the span above reorders our basis (thanks, jerk!) so we
1225     # need to put our "matrix basis" in the same order as the
1226     # (reordered) vector basis.
1227     S = tuple( _vec2mat(b) for b in W.basis() )
1228
1229     Qs = []
1230     for s in S:
1231         # Brute force the multiplication-by-s matrix by looping
1232         # through all elements of the basis and doing the computation
1233         # to find out what the corresponding row should be. BEWARE:
1234         # these multiplication tables won't be symmetric! It therefore
1235         # becomes REALLY IMPORTANT that the underlying algebra
1236         # constructor uses ROW vectors and not COLUMN vectors. That's
1237         # why we're computing rows here and not columns.
1238         Q_rows = []
1239         for t in S:
1240             this_row = _mat2vec((s*t + t*s)/2)
1241             Q_rows.append(W.coordinates(this_row))
1242         Q = matrix(field, W.dimension(), Q_rows)
1243         Qs.append(Q)
1244
1245     return (Qs, S)
1246
1247
1248 def _embed_complex_matrix(M):
1249     """
1250     Embed the n-by-n complex matrix ``M`` into the space of real
1251     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1252     bi` to the block matrix ``[[a,b],[-b,a]]``.
1253
1254     EXAMPLES::
1255
1256         sage: F = QuadraticField(-1,'i')
1257         sage: x1 = F(4 - 2*i)
1258         sage: x2 = F(1 + 2*i)
1259         sage: x3 = F(-i)
1260         sage: x4 = F(6)
1261         sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1262         sage: _embed_complex_matrix(M)
1263         [ 4 -2| 1  2]
1264         [ 2  4|-2  1]
1265         [-----+-----]
1266         [ 0 -1| 6  0]
1267         [ 1  0| 0  6]
1268
1269     TESTS:
1270
1271     Embedding is a homomorphism (isomorphism, in fact)::
1272
1273         sage: set_random_seed()
1274         sage: n = ZZ.random_element(5)
1275         sage: F = QuadraticField(-1, 'i')
1276         sage: X = random_matrix(F, n)
1277         sage: Y = random_matrix(F, n)
1278         sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1279         sage: expected = _embed_complex_matrix(X*Y)
1280         sage: actual == expected
1281         True
1282
1283     """
1284     n = M.nrows()
1285     if M.ncols() != n:
1286         raise ValueError("the matrix 'M' must be square")
1287     field = M.base_ring()
1288     blocks = []
1289     for z in M.list():
1290         a = z.real()
1291         b = z.imag()
1292         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1293
1294     # We can drop the imaginaries here.
1295     return block_matrix(field.base_ring(), n, blocks)
1296
1297
1298 def _unembed_complex_matrix(M):
1299     """
1300     The inverse of _embed_complex_matrix().
1301
1302     EXAMPLES::
1303
1304         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1305         ....:                 [-2,  1,  -4,  3],
1306         ....:                 [ 9,  10, 11, 12],
1307         ....:                 [-10, 9, -12, 11] ])
1308         sage: _unembed_complex_matrix(A)
1309         [  2*i + 1   4*i + 3]
1310         [ 10*i + 9 12*i + 11]
1311
1312     TESTS:
1313
1314     Unembedding is the inverse of embedding::
1315
1316         sage: set_random_seed()
1317         sage: F = QuadraticField(-1, 'i')
1318         sage: M = random_matrix(F, 3)
1319         sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1320         True
1321
1322     """
1323     n = ZZ(M.nrows())
1324     if M.ncols() != n:
1325         raise ValueError("the matrix 'M' must be square")
1326     if not n.mod(2).is_zero():
1327         raise ValueError("the matrix 'M' must be a complex embedding")
1328
1329     F = QuadraticField(-1, 'i')
1330     i = F.gen()
1331
1332     # Go top-left to bottom-right (reading order), converting every
1333     # 2-by-2 block we see to a single complex element.
1334     elements = []
1335     for k in xrange(n/2):
1336         for j in xrange(n/2):
1337             submat = M[2*k:2*k+2,2*j:2*j+2]
1338             if submat[0,0] != submat[1,1]:
1339                 raise ValueError('bad on-diagonal submatrix')
1340             if submat[0,1] != -submat[1,0]:
1341                 raise ValueError('bad off-diagonal submatrix')
1342             z = submat[0,0] + submat[0,1]*i
1343             elements.append(z)
1344
1345     return matrix(F, n/2, elements)
1346
1347
1348 def _embed_quaternion_matrix(M):
1349     """
1350     Embed the n-by-n quaternion matrix ``M`` into the space of real
1351     matrices of size 4n-by-4n by first sending each quaternion entry
1352     `z = a + bi + cj + dk` to the block-complex matrix
1353     ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1354     a real matrix.
1355
1356     EXAMPLES::
1357
1358         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1359         sage: i,j,k = Q.gens()
1360         sage: x = 1 + 2*i + 3*j + 4*k
1361         sage: M = matrix(Q, 1, [[x]])
1362         sage: _embed_quaternion_matrix(M)
1363         [ 1  2  3  4]
1364         [-2  1 -4  3]
1365         [-3  4  1 -2]
1366         [-4 -3  2  1]
1367
1368     Embedding is a homomorphism (isomorphism, in fact)::
1369
1370         sage: set_random_seed()
1371         sage: n = ZZ.random_element(5)
1372         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1373         sage: X = random_matrix(Q, n)
1374         sage: Y = random_matrix(Q, n)
1375         sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1376         sage: expected = _embed_quaternion_matrix(X*Y)
1377         sage: actual == expected
1378         True
1379
1380     """
1381     quaternions = M.base_ring()
1382     n = M.nrows()
1383     if M.ncols() != n:
1384         raise ValueError("the matrix 'M' must be square")
1385
1386     F = QuadraticField(-1, 'i')
1387     i = F.gen()
1388
1389     blocks = []
1390     for z in M.list():
1391         t = z.coefficient_tuple()
1392         a = t[0]
1393         b = t[1]
1394         c = t[2]
1395         d = t[3]
1396         cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1397                                     [-c + d*i, a - b*i]])
1398         blocks.append(_embed_complex_matrix(cplx_matrix))
1399
1400     # We should have real entries by now, so use the realest field
1401     # we've got for the return value.
1402     return block_matrix(quaternions.base_ring(), n, blocks)
1403
1404
1405 def _unembed_quaternion_matrix(M):
1406     """
1407     The inverse of _embed_quaternion_matrix().
1408
1409     EXAMPLES::
1410
1411         sage: M = matrix(QQ, [[ 1,  2,  3,  4],
1412         ....:                 [-2,  1, -4,  3],
1413         ....:                 [-3,  4,  1, -2],
1414         ....:                 [-4, -3,  2,  1]])
1415         sage: _unembed_quaternion_matrix(M)
1416         [1 + 2*i + 3*j + 4*k]
1417
1418     TESTS:
1419
1420     Unembedding is the inverse of embedding::
1421
1422         sage: set_random_seed()
1423         sage: Q = QuaternionAlgebra(QQ, -1, -1)
1424         sage: M = random_matrix(Q, 3)
1425         sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1426         True
1427
1428     """
1429     n = ZZ(M.nrows())
1430     if M.ncols() != n:
1431         raise ValueError("the matrix 'M' must be square")
1432     if not n.mod(4).is_zero():
1433         raise ValueError("the matrix 'M' must be a complex embedding")
1434
1435     Q = QuaternionAlgebra(QQ,-1,-1)
1436     i,j,k = Q.gens()
1437
1438     # Go top-left to bottom-right (reading order), converting every
1439     # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1440     # quaternion block.
1441     elements = []
1442     for l in xrange(n/4):
1443         for m in xrange(n/4):
1444             submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1445             if submat[0,0] != submat[1,1].conjugate():
1446                 raise ValueError('bad on-diagonal submatrix')
1447             if submat[0,1] != -submat[1,0].conjugate():
1448                 raise ValueError('bad off-diagonal submatrix')
1449             z  = submat[0,0].real() + submat[0,0].imag()*i
1450             z += submat[0,1].real()*j + submat[0,1].imag()*k
1451             elements.append(z)
1452
1453     return matrix(Q, n/4, elements)
1454
1455
1456 # The usual inner product on R^n.
1457 def _usual_ip(x,y):
1458     return x.vector().inner_product(y.vector())
1459
1460 # The inner product used for the real symmetric simple EJA.
1461 # We keep it as a separate function because e.g. the complex
1462 # algebra uses the same inner product, except divided by 2.
1463 def _matrix_ip(X,Y):
1464     X_mat = X.natural_representation()
1465     Y_mat = Y.natural_representation()
1466     return (X_mat*Y_mat).trace()
1467
1468
1469 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
1470     """
1471     The rank-n simple EJA consisting of real symmetric n-by-n
1472     matrices, the usual symmetric Jordan product, and the trace inner
1473     product. It has dimension `(n^2 + n)/2` over the reals.
1474
1475     EXAMPLES::
1476
1477         sage: J = RealSymmetricEJA(2)
1478         sage: e0, e1, e2 = J.gens()
1479         sage: e0*e0
1480         e0
1481         sage: e1*e1
1482         e0 + e2
1483         sage: e2*e2
1484         e2
1485
1486     TESTS:
1487
1488     The degree of this algebra is `(n^2 + n) / 2`::
1489
1490         sage: set_random_seed()
1491         sage: n = ZZ.random_element(1,5)
1492         sage: J = RealSymmetricEJA(n)
1493         sage: J.degree() == (n^2 + n)/2
1494         True
1495
1496     The Jordan multiplication is what we think it is::
1497
1498         sage: set_random_seed()
1499         sage: n = ZZ.random_element(1,5)
1500         sage: J = RealSymmetricEJA(n)
1501         sage: x = J.random_element()
1502         sage: y = J.random_element()
1503         sage: actual = (x*y).natural_representation()
1504         sage: X = x.natural_representation()
1505         sage: Y = y.natural_representation()
1506         sage: expected = (X*Y + Y*X)/2
1507         sage: actual == expected
1508         True
1509         sage: J(expected) == x*y
1510         True
1511
1512     """
1513     @staticmethod
1514     def __classcall_private__(cls, n, field=QQ):
1515         S = _real_symmetric_basis(n, field=field)
1516         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1517
1518         fdeja = super(RealSymmetricEJA, cls)
1519         return fdeja.__classcall_private__(cls,
1520                                            field,
1521                                            Qs,
1522                                            rank=n,
1523                                            natural_basis=T)
1524
1525     def inner_product(self, x, y):
1526         return _matrix_ip(x,y)
1527
1528
1529 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1530     """
1531     The rank-n simple EJA consisting of complex Hermitian n-by-n
1532     matrices over the real numbers, the usual symmetric Jordan product,
1533     and the real-part-of-trace inner product. It has dimension `n^2` over
1534     the reals.
1535
1536     TESTS:
1537
1538     The degree of this algebra is `n^2`::
1539
1540         sage: set_random_seed()
1541         sage: n = ZZ.random_element(1,5)
1542         sage: J = ComplexHermitianEJA(n)
1543         sage: J.degree() == n^2
1544         True
1545
1546     The Jordan multiplication is what we think it is::
1547
1548         sage: set_random_seed()
1549         sage: n = ZZ.random_element(1,5)
1550         sage: J = ComplexHermitianEJA(n)
1551         sage: x = J.random_element()
1552         sage: y = J.random_element()
1553         sage: actual = (x*y).natural_representation()
1554         sage: X = x.natural_representation()
1555         sage: Y = y.natural_representation()
1556         sage: expected = (X*Y + Y*X)/2
1557         sage: actual == expected
1558         True
1559         sage: J(expected) == x*y
1560         True
1561
1562     """
1563     @staticmethod
1564     def __classcall_private__(cls, n, field=QQ):
1565         S = _complex_hermitian_basis(n)
1566         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1567
1568         fdeja = super(ComplexHermitianEJA, cls)
1569         return fdeja.__classcall_private__(cls,
1570                                            field,
1571                                            Qs,
1572                                            rank=n,
1573                                            natural_basis=T)
1574
1575     def inner_product(self, x, y):
1576         # Since a+bi on the diagonal is represented as
1577         #
1578         #   a + bi  = [  a  b  ]
1579         #             [ -b  a  ],
1580         #
1581         # we'll double-count the "a" entries if we take the trace of
1582         # the embedding.
1583         return _matrix_ip(x,y)/2
1584
1585
1586 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1587     """
1588     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1589     matrices, the usual symmetric Jordan product, and the
1590     real-part-of-trace inner product. It has dimension `2n^2 - n` over
1591     the reals.
1592
1593     TESTS:
1594
1595     The degree of this algebra is `n^2`::
1596
1597         sage: set_random_seed()
1598         sage: n = ZZ.random_element(1,5)
1599         sage: J = QuaternionHermitianEJA(n)
1600         sage: J.degree() == 2*(n^2) - n
1601         True
1602
1603     The Jordan multiplication is what we think it is::
1604
1605         sage: set_random_seed()
1606         sage: n = ZZ.random_element(1,5)
1607         sage: J = QuaternionHermitianEJA(n)
1608         sage: x = J.random_element()
1609         sage: y = J.random_element()
1610         sage: actual = (x*y).natural_representation()
1611         sage: X = x.natural_representation()
1612         sage: Y = y.natural_representation()
1613         sage: expected = (X*Y + Y*X)/2
1614         sage: actual == expected
1615         True
1616         sage: J(expected) == x*y
1617         True
1618
1619     """
1620     @staticmethod
1621     def __classcall_private__(cls, n, field=QQ):
1622         S = _quaternion_hermitian_basis(n)
1623         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1624
1625         fdeja = super(QuaternionHermitianEJA, cls)
1626         return fdeja.__classcall_private__(cls,
1627                                            field,
1628                                            Qs,
1629                                            rank=n,
1630                                            natural_basis=T)
1631
1632     def inner_product(self, x, y):
1633         # Since a+bi+cj+dk on the diagonal is represented as
1634         #
1635         #   a + bi +cj + dk = [  a  b  c  d]
1636         #                     [ -b  a -d  c]
1637         #                     [ -c  d  a -b]
1638         #                     [ -d -c  b  a],
1639         #
1640         # we'll quadruple-count the "a" entries if we take the trace of
1641         # the embedding.
1642         return _matrix_ip(x,y)/4
1643
1644
1645 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
1646     """
1647     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1648     with the usual inner product and jordan product ``x*y =
1649     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1650     the reals.
1651
1652     EXAMPLES:
1653
1654     This multiplication table can be verified by hand::
1655
1656         sage: J = JordanSpinEJA(4)
1657         sage: e0,e1,e2,e3 = J.gens()
1658         sage: e0*e0
1659         e0
1660         sage: e0*e1
1661         e1
1662         sage: e0*e2
1663         e2
1664         sage: e0*e3
1665         e3
1666         sage: e1*e2
1667         0
1668         sage: e1*e3
1669         0
1670         sage: e2*e3
1671         0
1672
1673     """
1674     @staticmethod
1675     def __classcall_private__(cls, n, field=QQ):
1676         Qs = []
1677         id_matrix = identity_matrix(field, n)
1678         for i in xrange(n):
1679             ei = id_matrix.column(i)
1680             Qi = zero_matrix(field, n)
1681             Qi.set_row(0, ei)
1682             Qi.set_column(0, ei)
1683             Qi += diagonal_matrix(n, [ei[0]]*n)
1684             # The addition of the diagonal matrix adds an extra ei[0] in the
1685             # upper-left corner of the matrix.
1686             Qi[0,0] = Qi[0,0] * ~field(2)
1687             Qs.append(Qi)
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         fdeja = super(JordanSpinEJA, cls)
1693         return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2))
1694
1695     def inner_product(self, x, y):
1696         return _usual_ip(x,y)