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