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