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             ALGORITHM:
 497
 498             The usual way to do this is to check if the determinant is
 499             zero, but we need the characteristic polynomial for the
 500             determinant. The minimal polynomial is a lot easier to get,
 501             so we use Corollary 2 in Chapter V of Koecher to check
 502             whether or not the paren't algebra's zero element is a root
 503             of this element's minimal polynomial.
 504
 505             TESTS:
 506
 507             The identity element is always invertible::
 508
 509                 sage: set_random_seed()
 510                 sage: J = random_eja()
 511                 sage: J.one().is_invertible()
 512                 True
 513
 514             The zero element is never invertible::
 515
 516                 sage: set_random_seed()
 517                 sage: J = random_eja()
 518                 sage: J.zero().is_invertible()
 519                 False
 520
 521             """
 522             zero = self.parent().zero()
 523             p = self.minimal_polynomial()
 524             return not (p(zero) == zero)
 525
 526
 527         def is_nilpotent(self):
 528             """
 529             Return whether or not some power of this element is zero.
 530
 531             The superclass method won't work unless we're in an
 532             associative algebra, and we aren't. However, we generate
 533             an assocoative subalgebra and we're nilpotent there if and
 534             only if we're nilpotent here (probably).
 535
 536             TESTS:
 537
 538             The identity element is never nilpotent::
 539
 540                 sage: set_random_seed()
 541                 sage: random_eja().one().is_nilpotent()
 542                 False
 543
 544             The additive identity is always nilpotent::
 545
 546                 sage: set_random_seed()
 547                 sage: random_eja().zero().is_nilpotent()
 548                 True
 549
 550             """
 551             # The element we're going to call "is_nilpotent()" on.
 552             # Either myself, interpreted as an element of a finite-
 553             # dimensional algebra, or an element of an associative
 554             # subalgebra.
 555             elt = None
 556
 557             if self.parent().is_associative():
 558                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 559             else:
 560                 V = self.span_of_powers()
 561                 assoc_subalg = self.subalgebra_generated_by()
 562                 # Mis-design warning: the basis used for span_of_powers()
 563                 # and subalgebra_generated_by() must be the same, and in
 564                 # the same order!
 565                 elt = assoc_subalg(V.coordinates(self.vector()))
 566
 567             # Recursive call, but should work since elt lives in an
 568             # associative algebra.
 569             return elt.is_nilpotent()
 570
 571
 572         def is_regular(self):
 573             """
 574             Return whether or not this is a regular element.
 575
 576             EXAMPLES:
 577
 578             The identity element always has degree one, but any element
 579             linearly-independent from it is regular::
 580
 581                 sage: J = JordanSpinEJA(5)
 582                 sage: J.one().is_regular()
 583                 False
 584                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 585                 sage: for x in J.gens():
 586                 ....:     (J.one() + x).is_regular()
 587                 False
 588                 True
 589                 True
 590                 True
 591                 True
 592
 593             """
 594             return self.degree() == self.parent().rank()
 595
 596
 597         def degree(self):
 598             """
 599             Compute the degree of this element the straightforward way
 600             according to the definition; by appending powers to a list
 601             and figuring out its dimension (that is, whether or not
 602             they're linearly dependent).
 603
 604             EXAMPLES::
 605
 606                 sage: J = JordanSpinEJA(4)
 607                 sage: J.one().degree()
 608                 1
 609                 sage: e0,e1,e2,e3 = J.gens()
 610                 sage: (e0 - e1).degree()
 611                 2
 612
 613             In the spin factor algebra (of rank two), all elements that
 614             aren't multiples of the identity are regular::
 615
 616                 sage: set_random_seed()
 617                 sage: n = ZZ.random_element(1,10)
 618                 sage: J = JordanSpinEJA(n)
 619                 sage: x = J.random_element()
 620                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 621                 True
 622
 623             """
 624             return self.span_of_powers().dimension()
 625
 626
 627         def minimal_polynomial(self):
 628             """
 629             ALGORITHM:
 630
 631             We restrict ourselves to the associative subalgebra
 632             generated by this element, and then return the minimal
 633             polynomial of this element's operator matrix (in that
 634             subalgebra). This works by Baes Proposition 2.3.16.
 635
 636             EXAMPLES::
 637
 638                 sage: set_random_seed()
 639                 sage: x = random_eja().random_element()
 640                 sage: x.degree() == x.minimal_polynomial().degree()
 641                 True
 642
 643             ::
 644
 645                 sage: set_random_seed()
 646                 sage: x = random_eja().random_element()
 647                 sage: x.degree() == x.minimal_polynomial().degree()
 648                 True
 649
 650             The minimal polynomial and the characteristic polynomial coincide
 651             and are known (see Alizadeh, Example 11.11) for all elements of
 652             the spin factor algebra that aren't scalar multiples of the
 653             identity::
 654
 655                 sage: set_random_seed()
 656                 sage: n = ZZ.random_element(2,10)
 657                 sage: J = JordanSpinEJA(n)
 658                 sage: y = J.random_element()
 659                 sage: while y == y.coefficient(0)*J.one():
 660                 ....:     y = J.random_element()
 661                 sage: y0 = y.vector()[0]
 662                 sage: y_bar = y.vector()[1:]
 663                 sage: actual = y.minimal_polynomial()
 664                 sage: x = SR.symbol('x', domain='real')
 665                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 666                 sage: bool(actual == expected)
 667                 True
 668
 669             """
 670             V = self.span_of_powers()
 671             assoc_subalg = self.subalgebra_generated_by()
 672             # Mis-design warning: the basis used for span_of_powers()
 673             # and subalgebra_generated_by() must be the same, and in
 674             # the same order!
 675             elt = assoc_subalg(V.coordinates(self.vector()))
 676             return elt.operator_matrix().minimal_polynomial()
 677
 678
 679         def natural_representation(self):
 680             """
 681             Return a more-natural representation of this element.
 682
 683             Every finite-dimensional Euclidean Jordan Algebra is a
 684             direct sum of five simple algebras, four of which comprise
 685             Hermitian matrices. This method returns the original
 686             "natural" representation of this element as a Hermitian
 687             matrix, if it has one. If not, you get the usual representation.
 688
 689             EXAMPLES::
 690
 691                 sage: J = ComplexHermitianEJA(3)
 692                 sage: J.one()
 693                 e0 + e5 + e8
 694                 sage: J.one().natural_representation()
 695                 [1 0 0 0 0 0]
 696                 [0 1 0 0 0 0]
 697                 [0 0 1 0 0 0]
 698                 [0 0 0 1 0 0]
 699                 [0 0 0 0 1 0]
 700                 [0 0 0 0 0 1]
 701
 702             ::
 703
 704                 sage: J = QuaternionHermitianEJA(3)
 705                 sage: J.one()
 706                 e0 + e9 + e14
 707                 sage: J.one().natural_representation()
 708                 [1 0 0 0 0 0 0 0 0 0 0 0]
 709                 [0 1 0 0 0 0 0 0 0 0 0 0]
 710                 [0 0 1 0 0 0 0 0 0 0 0 0]
 711                 [0 0 0 1 0 0 0 0 0 0 0 0]
 712                 [0 0 0 0 1 0 0 0 0 0 0 0]
 713                 [0 0 0 0 0 1 0 0 0 0 0 0]
 714                 [0 0 0 0 0 0 1 0 0 0 0 0]
 715                 [0 0 0 0 0 0 0 1 0 0 0 0]
 716                 [0 0 0 0 0 0 0 0 1 0 0 0]
 717                 [0 0 0 0 0 0 0 0 0 1 0 0]
 718                 [0 0 0 0 0 0 0 0 0 0 1 0]
 719                 [0 0 0 0 0 0 0 0 0 0 0 1]
 720
 721             """
 722             B = self.parent().natural_basis()
 723             W = B[0].matrix_space()
 724             return W.linear_combination(zip(self.vector(), B))
 725
 726
 727         def operator_matrix(self):
 728             """
 729             Return the matrix that represents left- (or right-)
 730             multiplication by this element in the parent algebra.
 731
 732             We have to override this because the superclass method
 733             returns a matrix that acts on row vectors (that is, on
 734             the right).
 735
 736             EXAMPLES:
 737
 738             Test the first polarization identity from my notes, Koecher Chapter
 739             III, or from Baes (2.3)::
 740
 741                 sage: set_random_seed()
 742                 sage: J = random_eja()
 743                 sage: x = J.random_element()
 744                 sage: y = J.random_element()
 745                 sage: Lx = x.operator_matrix()
 746                 sage: Ly = y.operator_matrix()
 747                 sage: Lxx = (x*x).operator_matrix()
 748                 sage: Lxy = (x*y).operator_matrix()
 749                 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
 750                 True
 751
 752             Test the second polarization identity from my notes or from
 753             Baes (2.4)::
 754
 755                 sage: set_random_seed()
 756                 sage: J = random_eja()
 757                 sage: x = J.random_element()
 758                 sage: y = J.random_element()
 759                 sage: z = J.random_element()
 760                 sage: Lx = x.operator_matrix()
 761                 sage: Ly = y.operator_matrix()
 762                 sage: Lz = z.operator_matrix()
 763                 sage: Lzy = (z*y).operator_matrix()
 764                 sage: Lxy = (x*y).operator_matrix()
 765                 sage: Lxz = (x*z).operator_matrix()
 766                 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
 767                 True
 768
 769             Test the third polarization identity from my notes or from
 770             Baes (2.5)::
 771
 772                 sage: set_random_seed()
 773                 sage: J = random_eja()
 774                 sage: u = J.random_element()
 775                 sage: y = J.random_element()
 776                 sage: z = J.random_element()
 777                 sage: Lu = u.operator_matrix()
 778                 sage: Ly = y.operator_matrix()
 779                 sage: Lz = z.operator_matrix()
 780                 sage: Lzy = (z*y).operator_matrix()
 781                 sage: Luy = (u*y).operator_matrix()
 782                 sage: Luz = (u*z).operator_matrix()
 783                 sage: Luyz = (u*(y*z)).operator_matrix()
 784                 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
 785                 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
 786                 sage: bool(lhs == rhs)
 787                 True
 788
 789             """
 790             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 791             return fda_elt.matrix().transpose()
 792
 793
 794         def quadratic_representation(self, other=None):
 795             """
 796             Return the quadratic representation of this element.
 797
 798             EXAMPLES:
 799
 800             The explicit form in the spin factor algebra is given by
 801             Alizadeh's Example 11.12::
 802
 803                 sage: set_random_seed()
 804                 sage: n = ZZ.random_element(1,10)
 805                 sage: J = JordanSpinEJA(n)
 806                 sage: x = J.random_element()
 807                 sage: x_vec = x.vector()
 808                 sage: x0 = x_vec[0]
 809                 sage: x_bar = x_vec[1:]
 810                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 811                 sage: B = 2*x0*x_bar.row()
 812                 sage: C = 2*x0*x_bar.column()
 813                 sage: D = identity_matrix(QQ, n-1)
 814                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 815                 sage: D = D + 2*x_bar.tensor_product(x_bar)
 816                 sage: Q = block_matrix(2,2,[A,B,C,D])
 817                 sage: Q == x.quadratic_representation()
 818                 True
 819
 820             Test all of the properties from Theorem 11.2 in Alizadeh::
 821
 822                 sage: set_random_seed()
 823                 sage: J = random_eja()
 824                 sage: x = J.random_element()
 825                 sage: y = J.random_element()
 826
 827             Property 1:
 828
 829                 sage: actual = x.quadratic_representation(y)
 830                 sage: expected = ( (x+y).quadratic_representation()
 831                 ....:              -x.quadratic_representation()
 832                 ....:              -y.quadratic_representation() ) / 2
 833                 sage: actual == expected
 834                 True
 835
 836             Property 2:
 837
 838                 sage: alpha = QQ.random_element()
 839                 sage: actual = (alpha*x).quadratic_representation()
 840                 sage: expected = (alpha^2)*x.quadratic_representation()
 841                 sage: actual == expected
 842                 True
 843
 844             Property 5:
 845
 846                 sage: Qy = y.quadratic_representation()
 847                 sage: actual = J(Qy*x.vector()).quadratic_representation()
 848                 sage: expected = Qy*x.quadratic_representation()*Qy
 849                 sage: actual == expected
 850                 True
 851
 852             Property 6:
 853
 854                 sage: k = ZZ.random_element(1,10)
 855                 sage: actual = (x^k).quadratic_representation()
 856                 sage: expected = (x.quadratic_representation())^k
 857                 sage: actual == expected
 858                 True
 859
 860             """
 861             if other is None:
 862                 other=self
 863             elif not other in self.parent():
 864                 raise TypeError("'other' must live in the same algebra")
 865
 866             L = self.operator_matrix()
 867             M = other.operator_matrix()
 868             return ( L*M + M*L - (self*other).operator_matrix() )
 869
 870
 871         def span_of_powers(self):
 872             """
 873             Return the vector space spanned by successive powers of
 874             this element.
 875             """
 876             # The dimension of the subalgebra can't be greater than
 877             # the big algebra, so just put everything into a list
 878             # and let span() get rid of the excess.
 879             #
 880             # We do the extra ambient_vector_space() in case we're messing
 881             # with polynomials and the direct parent is a module.
 882             V = self.vector().parent().ambient_vector_space()
 883             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 884
 885
 886         def subalgebra_generated_by(self):
 887             """
 888             Return the associative subalgebra of the parent EJA generated
 889             by this element.
 890
 891             TESTS::
 892
 893                 sage: set_random_seed()
 894                 sage: x = random_eja().random_element()
 895                 sage: x.subalgebra_generated_by().is_associative()
 896                 True
 897
 898             Squaring in the subalgebra should be the same thing as
 899             squaring in the superalgebra::
 900
 901                 sage: set_random_seed()
 902                 sage: x = random_eja().random_element()
 903                 sage: u = x.subalgebra_generated_by().random_element()
 904                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
 905                 True
 906
 907             """
 908             # First get the subspace spanned by the powers of myself...
 909             V = self.span_of_powers()
 910             F = self.base_ring()
 911
 912             # Now figure out the entries of the right-multiplication
 913             # matrix for the successive basis elements b0, b1,... of
 914             # that subspace.
 915             mats = []
 916             for b_right in V.basis():
 917                 eja_b_right = self.parent()(b_right)
 918                 b_right_rows = []
 919                 # The first row of the right-multiplication matrix by
 920                 # b1 is what we get if we apply that matrix to b1. The
 921                 # second row of the right multiplication matrix by b1
 922                 # is what we get when we apply that matrix to b2...
 923                 #
 924                 # IMPORTANT: this assumes that all vectors are COLUMN
 925                 # vectors, unlike our superclass (which uses row vectors).
 926                 for b_left in V.basis():
 927                     eja_b_left = self.parent()(b_left)
 928                     # Multiply in the original EJA, but then get the
 929                     # coordinates from the subalgebra in terms of its
 930                     # basis.
 931                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 932                     b_right_rows.append(this_row)
 933                 b_right_matrix = matrix(F, b_right_rows)
 934                 mats.append(b_right_matrix)
 935
 936             # It's an algebra of polynomials in one element, and EJAs
 937             # are power-associative.
 938             #
 939             # TODO: choose generator names intelligently.
 940             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 941
 942
 943         def subalgebra_idempotent(self):
 944             """
 945             Find an idempotent in the associative subalgebra I generate
 946             using Proposition 2.3.5 in Baes.
 947
 948             TESTS::
 949
 950                 sage: set_random_seed()
 951                 sage: J = RealCartesianProductEJA(5)
 952                 sage: c = J.random_element().subalgebra_idempotent()
 953                 sage: c^2 == c
 954                 True
 955                 sage: J = JordanSpinEJA(5)
 956                 sage: c = J.random_element().subalgebra_idempotent()
 957                 sage: c^2 == c
 958                 True
 959
 960             """
 961             if self.is_nilpotent():
 962                 raise ValueError("this only works with non-nilpotent elements!")
 963
 964             V = self.span_of_powers()
 965             J = self.subalgebra_generated_by()
 966             # Mis-design warning: the basis used for span_of_powers()
 967             # and subalgebra_generated_by() must be the same, and in
 968             # the same order!
 969             u = J(V.coordinates(self.vector()))
 970
 971             # The image of the matrix of left-u^m-multiplication
 972             # will be minimal for some natural number s...
 973             s = 0
 974             minimal_dim = V.dimension()
 975             for i in xrange(1, V.dimension()):
 976                 this_dim = (u**i).operator_matrix().image().dimension()
 977                 if this_dim < minimal_dim:
 978                     minimal_dim = this_dim
 979                     s = i
 980
 981             # Now minimal_matrix should correspond to the smallest
 982             # non-zero subspace in Baes's (or really, Koecher's)
 983             # proposition.
 984             #
 985             # However, we need to restrict the matrix to work on the
 986             # subspace... or do we? Can't we just solve, knowing that
 987             # A(c) = u^(s+1) should have a solution in the big space,
 988             # too?
 989             #
 990             # Beware, solve_right() means that we're using COLUMN vectors.
 991             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 992             u_next = u**(s+1)
 993             A = u_next.operator_matrix()
 994             c_coordinates = A.solve_right(u_next.vector())
 995
 996             # Now c_coordinates is the idempotent we want, but it's in
 997             # the coordinate system of the subalgebra.
 998             #
 999             # We need the basis for J, but as elements of the parent algebra.
1000             #
1001             basis = [self.parent(v) for v in V.basis()]
1002             return self.parent().linear_combination(zip(c_coordinates, basis))
1003
1004
1005         def trace(self):
1006             """
1007             Return my trace, the sum of my eigenvalues.
1008
1009             EXAMPLES::
1010
1011                 sage: J = JordanSpinEJA(3)
1012                 sage: e0,e1,e2 = J.gens()
1013                 sage: x = e0 + e1 + e2
1014                 sage: x.trace()
1015                 2
1016
1017             """
1018             cs = self.characteristic_polynomial().coefficients(sparse=False)
1019             if len(cs) >= 2:
1020                 return -1*cs[-2]
1021             else:
1022                 raise ValueError('charpoly had fewer than 2 coefficients')
1023
1024
1025         def trace_inner_product(self, other):
1026             """
1027             Return the trace inner product of myself and ``other``.
1028             """
1029             if not other in self.parent():
1030                 raise TypeError("'other' must live in the same algebra")
1031
1032             return (self*other).trace()
1033
1034
1035 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
1036     """
1037     Return the Euclidean Jordan Algebra corresponding to the set
1038     `R^n` under the Hadamard product.
1039
1040     Note: this is nothing more than the Cartesian product of ``n``
1041     copies of the spin algebra. Once Cartesian product algebras
1042     are implemented, this can go.
1043
1044     EXAMPLES:
1045
1046     This multiplication table can be verified by hand::
1047
1048         sage: J = RealCartesianProductEJA(3)
1049         sage: e0,e1,e2 = J.gens()
1050         sage: e0*e0
1051         e0
1052         sage: e0*e1
1053         0
1054         sage: e0*e2
1055         0
1056         sage: e1*e1
1057         e1
1058         sage: e1*e2
1059         0
1060         sage: e2*e2
1061         e2
1062
1063     """
1064     @staticmethod
1065     def __classcall_private__(cls, n, field=QQ):
1066         # The FiniteDimensionalAlgebra constructor takes a list of
1067         # matrices, the ith representing right multiplication by the ith
1068         # basis element in the vector space. So if e_1 = (1,0,0), then
1069         # right (Hadamard) multiplication of x by e_1 picks out the first
1070         # component of x; and likewise for the ith basis element e_i.
1071         Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i))
1072                for i in xrange(n) ]
1073
1074         fdeja = super(RealCartesianProductEJA, cls)
1075         return fdeja.__classcall_private__(cls, field, Qs, rank=n)
1076
1077     def inner_product(self, x, y):
1078         return _usual_ip(x,y)
1079
1080
1081 def random_eja():
1082     """
1083     Return a "random" finite-dimensional Euclidean Jordan Algebra.
1084
1085     ALGORITHM:
1086
1087     For now, we choose a random natural number ``n`` (greater than zero)
1088     and then give you back one of the following:
1089
1090       * The cartesian product of the rational numbers ``n`` times; this is
1091         ``QQ^n`` with the Hadamard product.
1092
1093       * The Jordan spin algebra on ``QQ^n``.
1094
1095       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1096         product.
1097
1098       * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1099         in the space of ``2n``-by-``2n`` real symmetric matrices.
1100
1101       * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1102         in the space of ``4n``-by-``4n`` real symmetric matrices.
1103
1104     Later this might be extended to return Cartesian products of the
1105     EJAs above.
1106
1107     TESTS::
1108
1109         sage: random_eja()
1110         Euclidean Jordan algebra of degree...
1111
1112     """
1113     n = ZZ.random_element(1,5)
1114     constructor = choice([RealCartesianProductEJA,
1115                           JordanSpinEJA,
1116                           RealSymmetricEJA,
1117                           ComplexHermitianEJA,
1118                           QuaternionHermitianEJA])
1119     return constructor(n, field=QQ)
1120
1121
1122
1123 def _real_symmetric_basis(n, field=QQ):
1124     """
1125     Return a basis for the space of real symmetric n-by-n matrices.
1126     """
1127     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1128     # coordinates.
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 = Eij
1135             else:
1136                 # Beware, orthogonal but not normalized!
1137                 Sij = Eij + Eij.transpose()
1138             S.append(Sij)
1139     return tuple(S)
1140
1141
1142 def _complex_hermitian_basis(n, field=QQ):
1143     """
1144     Returns a basis for the space of complex Hermitian n-by-n matrices.
1145
1146     TESTS::
1147
1148         sage: set_random_seed()
1149         sage: n = ZZ.random_element(1,5)
1150         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1151         True
1152
1153     """
1154     F = QuadraticField(-1, 'I')
1155     I = F.gen()
1156
1157     # This is like the symmetric case, but we need to be careful:
1158     #
1159     #   * We want conjugate-symmetry, not just symmetry.
1160     #   * The diagonal will (as a result) be real.
1161     #
1162     S = []
1163     for i in xrange(n):
1164         for j in xrange(i+1):
1165             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1166             if i == j:
1167                 Sij = _embed_complex_matrix(Eij)
1168                 S.append(Sij)
1169             else:
1170                 # Beware, orthogonal but not normalized! The second one
1171                 # has a minus because it's conjugated.
1172                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1173                 S.append(Sij_real)
1174                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1175                 S.append(Sij_imag)
1176     return tuple(S)
1177
1178
1179 def _quaternion_hermitian_basis(n, field=QQ):
1180     """
1181     Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1182
1183     TESTS::
1184
1185         sage: set_random_seed()
1186         sage: n = ZZ.random_element(1,5)
1187         sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1188         True
1189
1190     """
1191     Q = QuaternionAlgebra(QQ,-1,-1)
1192     I,J,K = Q.gens()
1193
1194     # This is like the symmetric case, but we need to be careful:
1195     #
1196     #   * We want conjugate-symmetry, not just symmetry.
1197     #   * The diagonal will (as a result) be real.
1198     #
1199     S = []
1200     for i in xrange(n):
1201         for j in xrange(i+1):
1202             Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1203             if i == j:
1204                 Sij = _embed_quaternion_matrix(Eij)
1205                 S.append(Sij)
1206             else:
1207                 # Beware, orthogonal but not normalized! The second,
1208                 # third, and fourth ones have a minus because they're
1209                 # conjugated.
1210                 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1211                 S.append(Sij_real)
1212                 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1213                 S.append(Sij_I)
1214                 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1215                 S.append(Sij_J)
1216                 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1217                 S.append(Sij_K)
1218     return tuple(S)
1219
1220
1221 def _mat2vec(m):
1222         return vector(m.base_ring(), m.list())
1223
1224 def _vec2mat(v):
1225         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1226
1227 def _multiplication_table_from_matrix_basis(basis):
1228     """
1229     At least three of the five simple Euclidean Jordan algebras have the
1230     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1231     multiplication on the right is matrix multiplication. Given a basis
1232     for the underlying matrix space, this function returns a
1233     multiplication table (obtained by looping through the basis
1234     elements) for an algebra of those matrices. A reordered copy
1235     of the basis is also returned to work around the fact that
1236     the ``span()`` in this function will change the order of the basis
1237     from what we think it is, to... something else.
1238     """
1239     # In S^2, for example, we nominally have four coordinates even
1240     # though the space is of dimension three only. The vector space V
1241     # is supposed to hold the entire long vector, and the subspace W
1242     # of V will be spanned by the vectors that arise from symmetric
1243     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1244     field = basis[0].base_ring()
1245     dimension = basis[0].nrows()
1246
1247     V = VectorSpace(field, dimension**2)
1248     W = V.span( _mat2vec(s) for s in basis )
1249
1250     # Taking the span above reorders our basis (thanks, jerk!) so we
1251     # need to put our "matrix basis" in the same order as the
1252     # (reordered) vector basis.
1253     S = tuple( _vec2mat(b) for b in W.basis() )
1254
1255     Qs = []
1256     for s in S:
1257         # Brute force the multiplication-by-s matrix by looping
1258         # through all elements of the basis and doing the computation
1259         # to find out what the corresponding row should be. BEWARE:
1260         # these multiplication tables won't be symmetric! It therefore
1261         # becomes REALLY IMPORTANT that the underlying algebra
1262         # constructor uses ROW vectors and not COLUMN vectors. That's
1263         # why we're computing rows here and not columns.
1264         Q_rows = []
1265         for t in S:
1266             this_row = _mat2vec((s*t + t*s)/2)
1267             Q_rows.append(W.coordinates(this_row))
1268         Q = matrix(field, W.dimension(), Q_rows)
1269         Qs.append(Q)
1270
1271     return (Qs, S)
1272
1273
1274 def _embed_complex_matrix(M):
1275     """
1276     Embed the n-by-n complex matrix ``M`` into the space of real
1277     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1278     bi` to the block matrix ``[[a,b],[-b,a]]``.
1279
1280     EXAMPLES::
1281
1282         sage: F = QuadraticField(-1,'i')
1283         sage: x1 = F(4 - 2*i)
1284         sage: x2 = F(1 + 2*i)
1285         sage: x3 = F(-i)
1286         sage: x4 = F(6)
1287         sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1288         sage: _embed_complex_matrix(M)
1289         [ 4 -2| 1  2]
1290         [ 2  4|-2  1]
1291         [-----+-----]
1292         [ 0 -1| 6  0]
1293         [ 1  0| 0  6]
1294
1295     TESTS:
1296
1297     Embedding is a homomorphism (isomorphism, in fact)::
1298
1299         sage: set_random_seed()
1300         sage: n = ZZ.random_element(5)
1301         sage: F = QuadraticField(-1, 'i')
1302         sage: X = random_matrix(F, n)
1303         sage: Y = random_matrix(F, n)
1304         sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1305         sage: expected = _embed_complex_matrix(X*Y)
1306         sage: actual == expected
1307         True
1308
1309     """
1310     n = M.nrows()
1311     if M.ncols() != n:
1312         raise ValueError("the matrix 'M' must be square")
1313     field = M.base_ring()
1314     blocks = []
1315     for z in M.list():
1316         a = z.real()
1317         b = z.imag()
1318         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1319
1320     # We can drop the imaginaries here.
1321     return block_matrix(field.base_ring(), n, blocks)
1322
1323
1324 def _unembed_complex_matrix(M):
1325     """
1326     The inverse of _embed_complex_matrix().
1327
1328     EXAMPLES::
1329
1330         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1331         ....:                 [-2,  1,  -4,  3],
1332         ....:                 [ 9,  10, 11, 12],
1333         ....:                 [-10, 9, -12, 11] ])
1334         sage: _unembed_complex_matrix(A)
1335         [  2*i + 1   4*i + 3]
1336         [ 10*i + 9 12*i + 11]
1337
1338     TESTS:
1339
1340     Unembedding is the inverse of embedding::
1341
1342         sage: set_random_seed()
1343         sage: F = QuadraticField(-1, 'i')
1344         sage: M = random_matrix(F, 3)
1345         sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1346         True
1347
1348     """
1349     n = ZZ(M.nrows())
1350     if M.ncols() != n:
1351         raise ValueError("the matrix 'M' must be square")
1352     if not n.mod(2).is_zero():
1353         raise ValueError("the matrix 'M' must be a complex embedding")
1354
1355     F = QuadraticField(-1, 'i')
1356     i = F.gen()
1357
1358     # Go top-left to bottom-right (reading order), converting every
1359     # 2-by-2 block we see to a single complex element.
1360     elements = []
1361     for k in xrange(n/2):
1362         for j in xrange(n/2):
1363             submat = M[2*k:2*k+2,2*j:2*j+2]
1364             if submat[0,0] != submat[1,1]:
1365                 raise ValueError('bad on-diagonal submatrix')
1366             if submat[0,1] != -submat[1,0]:
1367                 raise ValueError('bad off-diagonal submatrix')
1368             z = submat[0,0] + submat[0,1]*i
1369             elements.append(z)
1370
1371     return matrix(F, n/2, elements)
1372
1373
1374 def _embed_quaternion_matrix(M):
1375     """
1376     Embed the n-by-n quaternion matrix ``M`` into the space of real
1377     matrices of size 4n-by-4n by first sending each quaternion entry
1378     `z = a + bi + cj + dk` to the block-complex matrix
1379     ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1380     a real matrix.
1381
1382     EXAMPLES::
1383
1384         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1385         sage: i,j,k = Q.gens()
1386         sage: x = 1 + 2*i + 3*j + 4*k
1387         sage: M = matrix(Q, 1, [[x]])
1388         sage: _embed_quaternion_matrix(M)
1389         [ 1  2  3  4]
1390         [-2  1 -4  3]
1391         [-3  4  1 -2]
1392         [-4 -3  2  1]
1393
1394     Embedding is a homomorphism (isomorphism, in fact)::
1395
1396         sage: set_random_seed()
1397         sage: n = ZZ.random_element(5)
1398         sage: Q = QuaternionAlgebra(QQ,-1,-1)
1399         sage: X = random_matrix(Q, n)
1400         sage: Y = random_matrix(Q, n)
1401         sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1402         sage: expected = _embed_quaternion_matrix(X*Y)
1403         sage: actual == expected
1404         True
1405
1406     """
1407     quaternions = M.base_ring()
1408     n = M.nrows()
1409     if M.ncols() != n:
1410         raise ValueError("the matrix 'M' must be square")
1411
1412     F = QuadraticField(-1, 'i')
1413     i = F.gen()
1414
1415     blocks = []
1416     for z in M.list():
1417         t = z.coefficient_tuple()
1418         a = t[0]
1419         b = t[1]
1420         c = t[2]
1421         d = t[3]
1422         cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1423                                     [-c + d*i, a - b*i]])
1424         blocks.append(_embed_complex_matrix(cplx_matrix))
1425
1426     # We should have real entries by now, so use the realest field
1427     # we've got for the return value.
1428     return block_matrix(quaternions.base_ring(), n, blocks)
1429
1430
1431 def _unembed_quaternion_matrix(M):
1432     """
1433     The inverse of _embed_quaternion_matrix().
1434
1435     EXAMPLES::
1436
1437         sage: M = matrix(QQ, [[ 1,  2,  3,  4],
1438         ....:                 [-2,  1, -4,  3],
1439         ....:                 [-3,  4,  1, -2],
1440         ....:                 [-4, -3,  2,  1]])
1441         sage: _unembed_quaternion_matrix(M)
1442         [1 + 2*i + 3*j + 4*k]
1443
1444     TESTS:
1445
1446     Unembedding is the inverse of embedding::
1447
1448         sage: set_random_seed()
1449         sage: Q = QuaternionAlgebra(QQ, -1, -1)
1450         sage: M = random_matrix(Q, 3)
1451         sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1452         True
1453
1454     """
1455     n = ZZ(M.nrows())
1456     if M.ncols() != n:
1457         raise ValueError("the matrix 'M' must be square")
1458     if not n.mod(4).is_zero():
1459         raise ValueError("the matrix 'M' must be a complex embedding")
1460
1461     Q = QuaternionAlgebra(QQ,-1,-1)
1462     i,j,k = Q.gens()
1463
1464     # Go top-left to bottom-right (reading order), converting every
1465     # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1466     # quaternion block.
1467     elements = []
1468     for l in xrange(n/4):
1469         for m in xrange(n/4):
1470             submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1471             if submat[0,0] != submat[1,1].conjugate():
1472                 raise ValueError('bad on-diagonal submatrix')
1473             if submat[0,1] != -submat[1,0].conjugate():
1474                 raise ValueError('bad off-diagonal submatrix')
1475             z  = submat[0,0].real() + submat[0,0].imag()*i
1476             z += submat[0,1].real()*j + submat[0,1].imag()*k
1477             elements.append(z)
1478
1479     return matrix(Q, n/4, elements)
1480
1481
1482 # The usual inner product on R^n.
1483 def _usual_ip(x,y):
1484     return x.vector().inner_product(y.vector())
1485
1486 # The inner product used for the real symmetric simple EJA.
1487 # We keep it as a separate function because e.g. the complex
1488 # algebra uses the same inner product, except divided by 2.
1489 def _matrix_ip(X,Y):
1490     X_mat = X.natural_representation()
1491     Y_mat = Y.natural_representation()
1492     return (X_mat*Y_mat).trace()
1493
1494
1495 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
1496     """
1497     The rank-n simple EJA consisting of real symmetric n-by-n
1498     matrices, the usual symmetric Jordan product, and the trace inner
1499     product. It has dimension `(n^2 + n)/2` over the reals.
1500
1501     EXAMPLES::
1502
1503         sage: J = RealSymmetricEJA(2)
1504         sage: e0, e1, e2 = J.gens()
1505         sage: e0*e0
1506         e0
1507         sage: e1*e1
1508         e0 + e2
1509         sage: e2*e2
1510         e2
1511
1512     TESTS:
1513
1514     The degree of this algebra is `(n^2 + n) / 2`::
1515
1516         sage: set_random_seed()
1517         sage: n = ZZ.random_element(1,5)
1518         sage: J = RealSymmetricEJA(n)
1519         sage: J.degree() == (n^2 + n)/2
1520         True
1521
1522     The Jordan multiplication is what we think it is::
1523
1524         sage: set_random_seed()
1525         sage: n = ZZ.random_element(1,5)
1526         sage: J = RealSymmetricEJA(n)
1527         sage: x = J.random_element()
1528         sage: y = J.random_element()
1529         sage: actual = (x*y).natural_representation()
1530         sage: X = x.natural_representation()
1531         sage: Y = y.natural_representation()
1532         sage: expected = (X*Y + Y*X)/2
1533         sage: actual == expected
1534         True
1535         sage: J(expected) == x*y
1536         True
1537
1538     """
1539     @staticmethod
1540     def __classcall_private__(cls, n, field=QQ):
1541         S = _real_symmetric_basis(n, field=field)
1542         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1543
1544         fdeja = super(RealSymmetricEJA, cls)
1545         return fdeja.__classcall_private__(cls,
1546                                            field,
1547                                            Qs,
1548                                            rank=n,
1549                                            natural_basis=T)
1550
1551     def inner_product(self, x, y):
1552         return _matrix_ip(x,y)
1553
1554
1555 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1556     """
1557     The rank-n simple EJA consisting of complex Hermitian n-by-n
1558     matrices over the real numbers, the usual symmetric Jordan product,
1559     and the real-part-of-trace inner product. It has dimension `n^2` over
1560     the reals.
1561
1562     TESTS:
1563
1564     The degree of this algebra is `n^2`::
1565
1566         sage: set_random_seed()
1567         sage: n = ZZ.random_element(1,5)
1568         sage: J = ComplexHermitianEJA(n)
1569         sage: J.degree() == n^2
1570         True
1571
1572     The Jordan multiplication is what we think it is::
1573
1574         sage: set_random_seed()
1575         sage: n = ZZ.random_element(1,5)
1576         sage: J = ComplexHermitianEJA(n)
1577         sage: x = J.random_element()
1578         sage: y = J.random_element()
1579         sage: actual = (x*y).natural_representation()
1580         sage: X = x.natural_representation()
1581         sage: Y = y.natural_representation()
1582         sage: expected = (X*Y + Y*X)/2
1583         sage: actual == expected
1584         True
1585         sage: J(expected) == x*y
1586         True
1587
1588     """
1589     @staticmethod
1590     def __classcall_private__(cls, n, field=QQ):
1591         S = _complex_hermitian_basis(n)
1592         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1593
1594         fdeja = super(ComplexHermitianEJA, cls)
1595         return fdeja.__classcall_private__(cls,
1596                                            field,
1597                                            Qs,
1598                                            rank=n,
1599                                            natural_basis=T)
1600
1601     def inner_product(self, x, y):
1602         # Since a+bi on the diagonal is represented as
1603         #
1604         #   a + bi  = [  a  b  ]
1605         #             [ -b  a  ],
1606         #
1607         # we'll double-count the "a" entries if we take the trace of
1608         # the embedding.
1609         return _matrix_ip(x,y)/2
1610
1611
1612 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1613     """
1614     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1615     matrices, the usual symmetric Jordan product, and the
1616     real-part-of-trace inner product. It has dimension `2n^2 - n` over
1617     the reals.
1618
1619     TESTS:
1620
1621     The degree of this algebra is `n^2`::
1622
1623         sage: set_random_seed()
1624         sage: n = ZZ.random_element(1,5)
1625         sage: J = QuaternionHermitianEJA(n)
1626         sage: J.degree() == 2*(n^2) - n
1627         True
1628
1629     The Jordan multiplication is what we think it is::
1630
1631         sage: set_random_seed()
1632         sage: n = ZZ.random_element(1,5)
1633         sage: J = QuaternionHermitianEJA(n)
1634         sage: x = J.random_element()
1635         sage: y = J.random_element()
1636         sage: actual = (x*y).natural_representation()
1637         sage: X = x.natural_representation()
1638         sage: Y = y.natural_representation()
1639         sage: expected = (X*Y + Y*X)/2
1640         sage: actual == expected
1641         True
1642         sage: J(expected) == x*y
1643         True
1644
1645     """
1646     @staticmethod
1647     def __classcall_private__(cls, n, field=QQ):
1648         S = _quaternion_hermitian_basis(n)
1649         (Qs, T) = _multiplication_table_from_matrix_basis(S)
1650
1651         fdeja = super(QuaternionHermitianEJA, cls)
1652         return fdeja.__classcall_private__(cls,
1653                                            field,
1654                                            Qs,
1655                                            rank=n,
1656                                            natural_basis=T)
1657
1658     def inner_product(self, x, y):
1659         # Since a+bi+cj+dk on the diagonal is represented as
1660         #
1661         #   a + bi +cj + dk = [  a  b  c  d]
1662         #                     [ -b  a -d  c]
1663         #                     [ -c  d  a -b]
1664         #                     [ -d -c  b  a],
1665         #
1666         # we'll quadruple-count the "a" entries if we take the trace of
1667         # the embedding.
1668         return _matrix_ip(x,y)/4
1669
1670
1671 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
1672     """
1673     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1674     with the usual inner product and jordan product ``x*y =
1675     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1676     the reals.
1677
1678     EXAMPLES:
1679
1680     This multiplication table can be verified by hand::
1681
1682         sage: J = JordanSpinEJA(4)
1683         sage: e0,e1,e2,e3 = J.gens()
1684         sage: e0*e0
1685         e0
1686         sage: e0*e1
1687         e1
1688         sage: e0*e2
1689         e2
1690         sage: e0*e3
1691         e3
1692         sage: e1*e2
1693         0
1694         sage: e1*e3
1695         0
1696         sage: e2*e3
1697         0
1698
1699     """
1700     @staticmethod
1701     def __classcall_private__(cls, n, field=QQ):
1702         Qs = []
1703         id_matrix = identity_matrix(field, n)
1704         for i in xrange(n):
1705             ei = id_matrix.column(i)
1706             Qi = zero_matrix(field, n)
1707             Qi.set_row(0, ei)
1708             Qi.set_column(0, ei)
1709             Qi += diagonal_matrix(n, [ei[0]]*n)
1710             # The addition of the diagonal matrix adds an extra ei[0] in the
1711             # upper-left corner of the matrix.
1712             Qi[0,0] = Qi[0,0] * ~field(2)
1713             Qs.append(Qi)
1714
1715         # The rank of the spin algebra is two, unless we're in a
1716         # one-dimensional ambient space (because the rank is bounded by
1717         # the ambient dimension).
1718         fdeja = super(JordanSpinEJA, cls)
1719         return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2))
1720
1721     def inner_product(self, x, y):
1722         return _usual_ip(x,y)