mjo/eja/eja_algebra.py

   1 """
   2 Representations and constructions for Euclidean Jordan algebras.
   3
   4 A Euclidean Jordan algebra is a Jordan algebra that has some
   5 additional properties:
   6
   7   1.   It is finite-dimensional.
   8   2.   Its scalar field is the real numbers.
   9   3a.  An inner product is defined on it, and...
  10   3b.  That inner product is compatible with the Jordan product
  11        in the sense that `<x*y,z> = <y,x*z>` for all elements
  12        `x,y,z` in the algebra.
  13
  14 Every Euclidean Jordan algebra is formally-real: for any two elements
  15 `x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
  16 0`. Conversely, every finite-dimensional formally-real Jordan algebra
  17 can be made into a Euclidean Jordan algebra with an appropriate choice
  18 of inner-product.
  19
  20 Formally-real Jordan algebras were originally studied as a framework
  21 for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
  22 symmetric cone optimization, since every symmetric cone arises as the
  23 cone of squares in some Euclidean Jordan algebra.
  24
  25 It is known that every Euclidean Jordan algebra decomposes into an
  26 orthogonal direct sum (essentially, a Cartesian product) of simple
  27 algebras, and that moreover, up to Jordan-algebra isomorphism, there
  28 are only five families of simple algebras. We provide constructions
  29 for these simple algebras:
  30
  31   * :class:`BilinearFormEJA`
  32   * :class:`RealSymmetricEJA`
  33   * :class:`ComplexHermitianEJA`
  34   * :class:`QuaternionHermitianEJA`
  35
  36 Missing from this list is the algebra of three-by-three octononion
  37 Hermitian matrices, as there is (as of yet) no implementation of the
  38 octonions in SageMath. In addition to these, we provide two other
  39 example constructions,
  40
  41   * :class:`HadamardEJA`
  42   * :class:`TrivialEJA`
  43
  44 The Jordan spin algebra is a bilinear form algebra where the bilinear
  45 form is the identity. The Hadamard EJA is simply a Cartesian product
  46 of one-dimensional spin algebras. And last but not least, the trivial
  47 EJA is exactly what you think. Cartesian products of these are also
  48 supported using the usual ``cartesian_product()`` function; as a
  49 result, we support (up to isomorphism) all Euclidean Jordan algebras
  50 that don't involve octonions.
  51
  52 SETUP::
  53
  54     sage: from mjo.eja.eja_algebra import random_eja
  55
  56 EXAMPLES::
  57
  58     sage: random_eja()
  59     Euclidean Jordan algebra of dimension...
  60 """
  61
  62 from itertools import repeat
  63
  64 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
  65 from sage.categories.magmatic_algebras import MagmaticAlgebras
  66 from sage.categories.sets_cat import cartesian_product
  67 from sage.combinat.free_module import CombinatorialFreeModule
  68 from sage.matrix.constructor import matrix
  69 from sage.matrix.matrix_space import MatrixSpace
  70 from sage.misc.cachefunc import cached_method
  71 from sage.misc.table import table
  72 from sage.modules.free_module import FreeModule, VectorSpace
  73 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
  74                             PolynomialRing,
  75                             QuadraticField)
  76 from mjo.eja.eja_element import FiniteDimensionalEJAElement
  77 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
  78 from mjo.eja.eja_utils import _all2list, _mat2vec
  79
  80 class FiniteDimensionalEJA(CombinatorialFreeModule):
  81     r"""
  82     A finite-dimensional Euclidean Jordan algebra.
  83
  84     INPUT:
  85
  86       - ``basis`` -- a tuple; a tuple of basis elements in "matrix
  87         form," which must be the same form as the arguments to
  88         ``jordan_product`` and ``inner_product``. In reality, "matrix
  89         form" can be either vectors, matrices, or a Cartesian product
  90         (ordered tuple) of vectors or matrices. All of these would
  91         ideally be vector spaces in sage with no special-casing
  92         needed; but in reality we turn vectors into column-matrices
  93         and Cartesian products `(a,b)` into column matrices
  94         `(a,b)^{T}` after converting `a` and `b` themselves.
  95
  96       - ``jordan_product`` -- a function; afunction of two ``basis``
  97         elements (in matrix form) that returns their jordan product,
  98         also in matrix form; this will be applied to ``basis`` to
  99         compute a multiplication table for the algebra.
 100
 101       - ``inner_product`` -- a function; a function of two ``basis``
 102         elements (in matrix form) that returns their inner
 103         product. This will be applied to ``basis`` to compute an
 104         inner-product table (basically a matrix) for this algebra.
 105
 106       - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
 107         field for the algebra.
 108
 109       - ``orthonormalize`` -- boolean (default: ``True``); whether or
 110         not to orthonormalize the basis. Doing so is expensive and
 111         generally rules out using the rationals as your ``field``, but
 112         is required for spectral decompositions.
 113
 114     SETUP::
 115
 116         sage: from mjo.eja.eja_algebra import random_eja
 117
 118     TESTS:
 119
 120     We should compute that an element subalgebra is associative even
 121     if we circumvent the element method::
 122
 123         sage: set_random_seed()
 124         sage: J = random_eja(field=QQ,orthonormalize=False)
 125         sage: x = J.random_element()
 126         sage: A = x.subalgebra_generated_by(orthonormalize=False)
 127         sage: basis = tuple(b.superalgebra_element() for b in A.basis())
 128         sage: J.subalgebra(basis, orthonormalize=False).is_associative()
 129         True
 130
 131     """
 132     Element = FiniteDimensionalEJAElement
 133
 134     def __init__(self,
 135                  basis,
 136                  jordan_product,
 137                  inner_product,
 138                  field=AA,
 139                  orthonormalize=True,
 140                  associative=None,
 141                  cartesian_product=False,
 142                  check_field=True,
 143                  check_axioms=True,
 144                  prefix="b"):
 145
 146         n = len(basis)
 147
 148         if check_field:
 149             if not field.is_subring(RR):
 150                 # Note: this does return true for the real algebraic
 151                 # field, the rationals, and any quadratic field where
 152                 # we've specified a real embedding.
 153                 raise ValueError("scalar field is not real")
 154
 155         if check_axioms:
 156             # Check commutativity of the Jordan and inner-products.
 157             # This has to be done before we build the multiplication
 158             # and inner-product tables/matrices, because we take
 159             # advantage of symmetry in the process.
 160             if not all( jordan_product(bi,bj) == jordan_product(bj,bi)
 161                         for bi in basis
 162                         for bj in basis ):
 163                 raise ValueError("Jordan product is not commutative")
 164
 165             if not all( inner_product(bi,bj) == inner_product(bj,bi)
 166                         for bi in basis
 167                         for bj in basis ):
 168                 raise ValueError("inner-product is not commutative")
 169
 170
 171         category = MagmaticAlgebras(field).FiniteDimensional()
 172         category = category.WithBasis().Unital().Commutative()
 173
 174         if associative is None:
 175             # We should figure it out. As with check_axioms, we have to do
 176             # this without the help of the _jordan_product_is_associative()
 177             # method because we need to know the category before we
 178             # initialize the algebra.
 179             associative = all( jordan_product(jordan_product(bi,bj),bk)
 180                                ==
 181                                jordan_product(bi,jordan_product(bj,bk))
 182                                for bi in basis
 183                                for bj in basis
 184                                for bk in basis)
 185
 186         if associative:
 187             # Element subalgebras can take advantage of this.
 188             category = category.Associative()
 189         if cartesian_product:
 190             # Use join() here because otherwise we only get the
 191             # "Cartesian product of..." and not the things themselves.
 192             category = category.join([category,
 193                                       category.CartesianProducts()])
 194
 195         # Call the superclass constructor so that we can use its from_vector()
 196         # method to build our multiplication table.
 197         CombinatorialFreeModule.__init__(self,
 198                                          field,
 199                                          range(n),
 200                                          prefix=prefix,
 201                                          category=category,
 202                                          bracket=False)
 203
 204         # Now comes all of the hard work. We'll be constructing an
 205         # ambient vector space V that our (vectorized) basis lives in,
 206         # as well as a subspace W of V spanned by those (vectorized)
 207         # basis elements. The W-coordinates are the coefficients that
 208         # we see in things like x = 1*b1 + 2*b2.
 209         vector_basis = basis
 210
 211         degree = 0
 212         if n > 0:
 213             degree = len(_all2list(basis[0]))
 214
 215         # Build an ambient space that fits our matrix basis when
 216         # written out as "long vectors."
 217         V = VectorSpace(field, degree)
 218
 219         # The matrix that will hole the orthonormal -> unorthonormal
 220         # coordinate transformation.
 221         self._deortho_matrix = None
 222
 223         if orthonormalize:
 224             # Save a copy of the un-orthonormalized basis for later.
 225             # Convert it to ambient V (vector) coordinates while we're
 226             # at it, because we'd have to do it later anyway.
 227             deortho_vector_basis = tuple( V(_all2list(b)) for b in basis )
 228
 229             from mjo.eja.eja_utils import gram_schmidt
 230             basis = tuple(gram_schmidt(basis, inner_product))
 231
 232         # Save the (possibly orthonormalized) matrix basis for
 233         # later...
 234         self._matrix_basis = basis
 235
 236         # Now create the vector space for the algebra, which will have
 237         # its own set of non-ambient coordinates (in terms of the
 238         # supplied basis).
 239         vector_basis = tuple( V(_all2list(b)) for b in basis )
 240         W = V.span_of_basis( vector_basis, check=check_axioms)
 241
 242         if orthonormalize:
 243             # Now "W" is the vector space of our algebra coordinates. The
 244             # variables "X1", "X2",...  refer to the entries of vectors in
 245             # W. Thus to convert back and forth between the orthonormal
 246             # coordinates and the given ones, we need to stick the original
 247             # basis in W.
 248             U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
 249             self._deortho_matrix = matrix( U.coordinate_vector(q)
 250                                            for q in vector_basis )
 251
 252
 253         # Now we actually compute the multiplication and inner-product
 254         # tables/matrices using the possibly-orthonormalized basis.
 255         self._inner_product_matrix = matrix.identity(field, n)
 256         self._multiplication_table = [ [0 for j in range(i+1)]
 257                                        for i in range(n) ]
 258
 259         # Note: the Jordan and inner-products are defined in terms
 260         # of the ambient basis. It's important that their arguments
 261         # are in ambient coordinates as well.
 262         for i in range(n):
 263             for j in range(i+1):
 264                 # ortho basis w.r.t. ambient coords
 265                 q_i = basis[i]
 266                 q_j = basis[j]
 267
 268                 # The jordan product returns a matrixy answer, so we
 269                 # have to convert it to the algebra coordinates.
 270                 elt = jordan_product(q_i, q_j)
 271                 elt = W.coordinate_vector(V(_all2list(elt)))
 272                 self._multiplication_table[i][j] = self.from_vector(elt)
 273
 274                 if not orthonormalize:
 275                     # If we're orthonormalizing the basis with respect
 276                     # to an inner-product, then the inner-product
 277                     # matrix with respect to the resulting basis is
 278                     # just going to be the identity.
 279                     ip = inner_product(q_i, q_j)
 280                     self._inner_product_matrix[i,j] = ip
 281                     self._inner_product_matrix[j,i] = ip
 282
 283         self._inner_product_matrix._cache = {'hermitian': True}
 284         self._inner_product_matrix.set_immutable()
 285
 286         if check_axioms:
 287             if not self._is_jordanian():
 288                 raise ValueError("Jordan identity does not hold")
 289             if not self._inner_product_is_associative():
 290                 raise ValueError("inner product is not associative")
 291
 292
 293     def _coerce_map_from_base_ring(self):
 294         """
 295         Disable the map from the base ring into the algebra.
 296
 297         Performing a nonsense conversion like this automatically
 298         is counterpedagogical. The fallback is to try the usual
 299         element constructor, which should also fail.
 300
 301         SETUP::
 302
 303             sage: from mjo.eja.eja_algebra import random_eja
 304
 305         TESTS::
 306
 307             sage: set_random_seed()
 308             sage: J = random_eja()
 309             sage: J(1)
 310             Traceback (most recent call last):
 311             ...
 312             ValueError: not an element of this algebra
 313
 314         """
 315         return None
 316
 317
 318     def product_on_basis(self, i, j):
 319         r"""
 320         Returns the Jordan product of the `i` and `j`th basis elements.
 321
 322         This completely defines the Jordan product on the algebra, and
 323         is used direclty by our superclass machinery to implement
 324         :meth:`product`.
 325
 326         SETUP::
 327
 328             sage: from mjo.eja.eja_algebra import random_eja
 329
 330         TESTS::
 331
 332             sage: set_random_seed()
 333             sage: J = random_eja()
 334             sage: n = J.dimension()
 335             sage: bi = J.zero()
 336             sage: bj = J.zero()
 337             sage: bi_bj = J.zero()*J.zero()
 338             sage: if n > 0:
 339             ....:     i = ZZ.random_element(n)
 340             ....:     j = ZZ.random_element(n)
 341             ....:     bi = J.monomial(i)
 342             ....:     bj = J.monomial(j)
 343             ....:     bi_bj = J.product_on_basis(i,j)
 344             sage: bi*bj == bi_bj
 345             True
 346
 347         """
 348         # We only stored the lower-triangular portion of the
 349         # multiplication table.
 350         if j <= i:
 351             return self._multiplication_table[i][j]
 352         else:
 353             return self._multiplication_table[j][i]
 354
 355     def inner_product(self, x, y):
 356         """
 357         The inner product associated with this Euclidean Jordan algebra.
 358
 359         Defaults to the trace inner product, but can be overridden by
 360         subclasses if they are sure that the necessary properties are
 361         satisfied.
 362
 363         SETUP::
 364
 365             sage: from mjo.eja.eja_algebra import (random_eja,
 366             ....:                                  HadamardEJA,
 367             ....:                                  BilinearFormEJA)
 368
 369         EXAMPLES:
 370
 371         Our inner product is "associative," which means the following for
 372         a symmetric bilinear form::
 373
 374             sage: set_random_seed()
 375             sage: J = random_eja()
 376             sage: x,y,z = J.random_elements(3)
 377             sage: (x*y).inner_product(z) == y.inner_product(x*z)
 378             True
 379
 380         TESTS:
 381
 382         Ensure that this is the usual inner product for the algebras
 383         over `R^n`::
 384
 385             sage: set_random_seed()
 386             sage: J = HadamardEJA.random_instance()
 387             sage: x,y = J.random_elements(2)
 388             sage: actual = x.inner_product(y)
 389             sage: expected = x.to_vector().inner_product(y.to_vector())
 390             sage: actual == expected
 391             True
 392
 393         Ensure that this is one-half of the trace inner-product in a
 394         BilinearFormEJA that isn't just the reals (when ``n`` isn't
 395         one). This is in Faraut and Koranyi, and also my "On the
 396         symmetry..." paper::
 397
 398             sage: set_random_seed()
 399             sage: J = BilinearFormEJA.random_instance()
 400             sage: n = J.dimension()
 401             sage: x = J.random_element()
 402             sage: y = J.random_element()
 403             sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
 404             True
 405
 406         """
 407         B = self._inner_product_matrix
 408         return (B*x.to_vector()).inner_product(y.to_vector())
 409
 410
 411     def is_associative(self):
 412         r"""
 413         Return whether or not this algebra's Jordan product is associative.
 414
 415         SETUP::
 416
 417             sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
 418
 419         EXAMPLES::
 420
 421             sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
 422             sage: J.is_associative()
 423             False
 424             sage: x = sum(J.gens())
 425             sage: A = x.subalgebra_generated_by(orthonormalize=False)
 426             sage: A.is_associative()
 427             True
 428
 429         """
 430         return "Associative" in self.category().axioms()
 431
 432     def _is_commutative(self):
 433         r"""
 434         Whether or not this algebra's multiplication table is commutative.
 435
 436         This method should of course always return ``True``, unless
 437         this algebra was constructed with ``check_axioms=False`` and
 438         passed an invalid multiplication table.
 439         """
 440         return all( x*y == y*x for x in self.gens() for y in self.gens() )
 441
 442     def _is_jordanian(self):
 443         r"""
 444         Whether or not this algebra's multiplication table respects the
 445         Jordan identity `(x^{2})(xy) = x(x^{2}y)`.
 446
 447         We only check one arrangement of `x` and `y`, so for a
 448         ``True`` result to be truly true, you should also check
 449         :meth:`_is_commutative`. This method should of course always
 450         return ``True``, unless this algebra was constructed with
 451         ``check_axioms=False`` and passed an invalid multiplication table.
 452         """
 453         return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
 454                     ==
 455                     (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
 456                     for i in range(self.dimension())
 457                     for j in range(self.dimension()) )
 458
 459     def _jordan_product_is_associative(self):
 460         r"""
 461         Return whether or not this algebra's Jordan product is
 462         associative; that is, whether or not `x*(y*z) = (x*y)*z`
 463         for all `x,y,x`.
 464
 465         This method should agree with :meth:`is_associative` unless
 466         you lied about the value of the ``associative`` parameter
 467         when you constructed the algebra.
 468
 469         SETUP::
 470
 471             sage: from mjo.eja.eja_algebra import (random_eja,
 472             ....:                                  RealSymmetricEJA,
 473             ....:                                  ComplexHermitianEJA,
 474             ....:                                  QuaternionHermitianEJA)
 475
 476         EXAMPLES::
 477
 478             sage: J = RealSymmetricEJA(4, orthonormalize=False)
 479             sage: J._jordan_product_is_associative()
 480             False
 481             sage: x = sum(J.gens())
 482             sage: A = x.subalgebra_generated_by()
 483             sage: A._jordan_product_is_associative()
 484             True
 485
 486         ::
 487
 488             sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
 489             sage: J._jordan_product_is_associative()
 490             False
 491             sage: x = sum(J.gens())
 492             sage: A = x.subalgebra_generated_by(orthonormalize=False)
 493             sage: A._jordan_product_is_associative()
 494             True
 495
 496         ::
 497
 498             sage: J = QuaternionHermitianEJA(2)
 499             sage: J._jordan_product_is_associative()
 500             False
 501             sage: x = sum(J.gens())
 502             sage: A = x.subalgebra_generated_by()
 503             sage: A._jordan_product_is_associative()
 504             True
 505
 506         TESTS:
 507
 508         The values we've presupplied to the constructors agree with
 509         the computation::
 510
 511             sage: set_random_seed()
 512             sage: J = random_eja()
 513             sage: J.is_associative() == J._jordan_product_is_associative()
 514             True
 515
 516         """
 517         R = self.base_ring()
 518
 519         # Used to check whether or not something is zero.
 520         epsilon = R.zero()
 521         if not R.is_exact():
 522             # I don't know of any examples that make this magnitude
 523             # necessary because I don't know how to make an
 524             # associative algebra when the element subalgebra
 525             # construction is unreliable (as it is over RDF; we can't
 526             # find the degree of an element because we can't compute
 527             # the rank of a matrix). But even multiplication of floats
 528             # is non-associative, so *some* epsilon is needed... let's
 529             # just take the one from _inner_product_is_associative?
 530             epsilon = 1e-15
 531
 532         for i in range(self.dimension()):
 533             for j in range(self.dimension()):
 534                 for k in range(self.dimension()):
 535                     x = self.monomial(i)
 536                     y = self.monomial(j)
 537                     z = self.monomial(k)
 538                     diff = (x*y)*z - x*(y*z)
 539
 540                     if diff.norm() > epsilon:
 541                         return False
 542
 543         return True
 544
 545     def _inner_product_is_associative(self):
 546         r"""
 547         Return whether or not this algebra's inner product `B` is
 548         associative; that is, whether or not `B(xy,z) = B(x,yz)`.
 549
 550         This method should of course always return ``True``, unless
 551         this algebra was constructed with ``check_axioms=False`` and
 552         passed an invalid Jordan or inner-product.
 553         """
 554         R = self.base_ring()
 555
 556         # Used to check whether or not something is zero.
 557         epsilon = R.zero()
 558         if not R.is_exact():
 559             # This choice is sufficient to allow the construction of
 560             # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
 561             epsilon = 1e-15
 562
 563         for i in range(self.dimension()):
 564             for j in range(self.dimension()):
 565                 for k in range(self.dimension()):
 566                     x = self.monomial(i)
 567                     y = self.monomial(j)
 568                     z = self.monomial(k)
 569                     diff = (x*y).inner_product(z) - x.inner_product(y*z)
 570
 571                     if diff.abs() > epsilon:
 572                         return False
 573
 574         return True
 575
 576     def _element_constructor_(self, elt):
 577         """
 578         Construct an element of this algebra from its vector or matrix
 579         representation.
 580
 581         This gets called only after the parent element _call_ method
 582         fails to find a coercion for the argument.
 583
 584         SETUP::
 585
 586             sage: from mjo.eja.eja_algebra import (random_eja,
 587             ....:                                  JordanSpinEJA,
 588             ....:                                  HadamardEJA,
 589             ....:                                  RealSymmetricEJA)
 590
 591         EXAMPLES:
 592
 593         The identity in `S^n` is converted to the identity in the EJA::
 594
 595             sage: J = RealSymmetricEJA(3)
 596             sage: I = matrix.identity(QQ,3)
 597             sage: J(I) == J.one()
 598             True
 599
 600         This skew-symmetric matrix can't be represented in the EJA::
 601
 602             sage: J = RealSymmetricEJA(3)
 603             sage: A = matrix(QQ,3, lambda i,j: i-j)
 604             sage: J(A)
 605             Traceback (most recent call last):
 606             ...
 607             ValueError: not an element of this algebra
 608
 609         Tuples work as well, provided that the matrix basis for the
 610         algebra consists of them::
 611
 612             sage: J1 = HadamardEJA(3)
 613             sage: J2 = RealSymmetricEJA(2)
 614             sage: J = cartesian_product([J1,J2])
 615             sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
 616             b1 + b5
 617
 618         TESTS:
 619
 620         Ensure that we can convert any element back and forth
 621         faithfully between its matrix and algebra representations::
 622
 623             sage: set_random_seed()
 624             sage: J = random_eja()
 625             sage: x = J.random_element()
 626             sage: J(x.to_matrix()) == x
 627             True
 628
 629         We cannot coerce elements between algebras just because their
 630         matrix representations are compatible::
 631
 632             sage: J1 = HadamardEJA(3)
 633             sage: J2 = JordanSpinEJA(3)
 634             sage: J2(J1.one())
 635             Traceback (most recent call last):
 636             ...
 637             ValueError: not an element of this algebra
 638             sage: J1(J2.zero())
 639             Traceback (most recent call last):
 640             ...
 641             ValueError: not an element of this algebra
 642         """
 643         msg = "not an element of this algebra"
 644         if elt in self.base_ring():
 645             # Ensure that no base ring -> algebra coercion is performed
 646             # by this method. There's some stupidity in sage that would
 647             # otherwise propagate to this method; for example, sage thinks
 648             # that the integer 3 belongs to the space of 2-by-2 matrices.
 649             raise ValueError(msg)
 650
 651         try:
 652             # Try to convert a vector into a column-matrix...
 653             elt = elt.column()
 654         except (AttributeError, TypeError):
 655             # and ignore failure, because we weren't really expecting
 656             # a vector as an argument anyway.
 657             pass
 658
 659         if elt not in self.matrix_space():
 660             raise ValueError(msg)
 661
 662         # Thanks for nothing! Matrix spaces aren't vector spaces in
 663         # Sage, so we have to figure out its matrix-basis coordinates
 664         # ourselves. We use the basis space's ring instead of the
 665         # element's ring because the basis space might be an algebraic
 666         # closure whereas the base ring of the 3-by-3 identity matrix
 667         # could be QQ instead of QQbar.
 668         #
 669         # And, we also have to handle Cartesian product bases (when
 670         # the matrix basis consists of tuples) here. The "good news"
 671         # is that we're already converting everything to long vectors,
 672         # and that strategy works for tuples as well.
 673         #
 674         # We pass check=False because the matrix basis is "guaranteed"
 675         # to be linearly independent... right? Ha ha.
 676         elt = _all2list(elt)
 677         V = VectorSpace(self.base_ring(), len(elt))
 678         W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()),
 679                              check=False)
 680
 681         try:
 682             coords = W.coordinate_vector(V(elt))
 683         except ArithmeticError:  # vector is not in free module
 684             raise ValueError(msg)
 685
 686         return self.from_vector(coords)
 687
 688     def _repr_(self):
 689         """
 690         Return a string representation of ``self``.
 691
 692         SETUP::
 693
 694             sage: from mjo.eja.eja_algebra import JordanSpinEJA
 695
 696         TESTS:
 697
 698         Ensure that it says what we think it says::
 699
 700             sage: JordanSpinEJA(2, field=AA)
 701             Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
 702             sage: JordanSpinEJA(3, field=RDF)
 703             Euclidean Jordan algebra of dimension 3 over Real Double Field
 704
 705         """
 706         fmt = "Euclidean Jordan algebra of dimension {} over {}"
 707         return fmt.format(self.dimension(), self.base_ring())
 708
 709
 710     @cached_method
 711     def characteristic_polynomial_of(self):
 712         """
 713         Return the algebra's "characteristic polynomial of" function,
 714         which is itself a multivariate polynomial that, when evaluated
 715         at the coordinates of some algebra element, returns that
 716         element's characteristic polynomial.
 717
 718         The resulting polynomial has `n+1` variables, where `n` is the
 719         dimension of this algebra. The first `n` variables correspond to
 720         the coordinates of an algebra element: when evaluated at the
 721         coordinates of an algebra element with respect to a certain
 722         basis, the result is a univariate polynomial (in the one
 723         remaining variable ``t``), namely the characteristic polynomial
 724         of that element.
 725
 726         SETUP::
 727
 728             sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA
 729
 730         EXAMPLES:
 731
 732         The characteristic polynomial in the spin algebra is given in
 733         Alizadeh, Example 11.11::
 734
 735             sage: J = JordanSpinEJA(3)
 736             sage: p = J.characteristic_polynomial_of(); p
 737             X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
 738             sage: xvec = J.one().to_vector()
 739             sage: p(*xvec)
 740             t^2 - 2*t + 1
 741
 742         By definition, the characteristic polynomial is a monic
 743         degree-zero polynomial in a rank-zero algebra. Note that
 744         Cayley-Hamilton is indeed satisfied since the polynomial
 745         ``1`` evaluates to the identity element of the algebra on
 746         any argument::
 747
 748             sage: J = TrivialEJA()
 749             sage: J.characteristic_polynomial_of()
 750             1
 751
 752         """
 753         r = self.rank()
 754         n = self.dimension()
 755
 756         # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1).
 757         a = self._charpoly_coefficients()
 758
 759         # We go to a bit of trouble here to reorder the
 760         # indeterminates, so that it's easier to evaluate the
 761         # characteristic polynomial at x's coordinates and get back
 762         # something in terms of t, which is what we want.
 763         S = PolynomialRing(self.base_ring(),'t')
 764         t = S.gen(0)
 765         if r > 0:
 766             R = a[0].parent()
 767             S = PolynomialRing(S, R.variable_names())
 768             t = S(t)
 769
 770         return (t**r + sum( a[k]*(t**k) for k in range(r) ))
 771
 772     def coordinate_polynomial_ring(self):
 773         r"""
 774         The multivariate polynomial ring in which this algebra's
 775         :meth:`characteristic_polynomial_of` lives.
 776
 777         SETUP::
 778
 779             sage: from mjo.eja.eja_algebra import (HadamardEJA,
 780             ....:                                  RealSymmetricEJA)
 781
 782         EXAMPLES::
 783
 784             sage: J = HadamardEJA(2)
 785             sage: J.coordinate_polynomial_ring()
 786             Multivariate Polynomial Ring in X1, X2...
 787             sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
 788             sage: J.coordinate_polynomial_ring()
 789             Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6...
 790
 791         """
 792         var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) )
 793         return PolynomialRing(self.base_ring(), var_names)
 794
 795     def inner_product(self, x, y):
 796         """
 797         The inner product associated with this Euclidean Jordan algebra.
 798
 799         Defaults to the trace inner product, but can be overridden by
 800         subclasses if they are sure that the necessary properties are
 801         satisfied.
 802
 803         SETUP::
 804
 805             sage: from mjo.eja.eja_algebra import (random_eja,
 806             ....:                                  HadamardEJA,
 807             ....:                                  BilinearFormEJA)
 808
 809         EXAMPLES:
 810
 811         Our inner product is "associative," which means the following for
 812         a symmetric bilinear form::
 813
 814             sage: set_random_seed()
 815             sage: J = random_eja()
 816             sage: x,y,z = J.random_elements(3)
 817             sage: (x*y).inner_product(z) == y.inner_product(x*z)
 818             True
 819
 820         TESTS:
 821
 822         Ensure that this is the usual inner product for the algebras
 823         over `R^n`::
 824
 825             sage: set_random_seed()
 826             sage: J = HadamardEJA.random_instance()
 827             sage: x,y = J.random_elements(2)
 828             sage: actual = x.inner_product(y)
 829             sage: expected = x.to_vector().inner_product(y.to_vector())
 830             sage: actual == expected
 831             True
 832
 833         Ensure that this is one-half of the trace inner-product in a
 834         BilinearFormEJA that isn't just the reals (when ``n`` isn't
 835         one). This is in Faraut and Koranyi, and also my "On the
 836         symmetry..." paper::
 837
 838             sage: set_random_seed()
 839             sage: J = BilinearFormEJA.random_instance()
 840             sage: n = J.dimension()
 841             sage: x = J.random_element()
 842             sage: y = J.random_element()
 843             sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
 844             True
 845         """
 846         B = self._inner_product_matrix
 847         return (B*x.to_vector()).inner_product(y.to_vector())
 848
 849
 850     def is_trivial(self):
 851         """
 852         Return whether or not this algebra is trivial.
 853
 854         A trivial algebra contains only the zero element.
 855
 856         SETUP::
 857
 858             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 859             ....:                                  TrivialEJA)
 860
 861         EXAMPLES::
 862
 863             sage: J = ComplexHermitianEJA(3)
 864             sage: J.is_trivial()
 865             False
 866
 867         ::
 868
 869             sage: J = TrivialEJA()
 870             sage: J.is_trivial()
 871             True
 872
 873         """
 874         return self.dimension() == 0
 875
 876
 877     def multiplication_table(self):
 878         """
 879         Return a visual representation of this algebra's multiplication
 880         table (on basis elements).
 881
 882         SETUP::
 883
 884             sage: from mjo.eja.eja_algebra import JordanSpinEJA
 885
 886         EXAMPLES::
 887
 888             sage: J = JordanSpinEJA(4)
 889             sage: J.multiplication_table()
 890             +----++----+----+----+----+
 891             | *  || b0 | b1 | b2 | b3 |
 892             +====++====+====+====+====+
 893             | b0 || b0 | b1 | b2 | b3 |
 894             +----++----+----+----+----+
 895             | b1 || b1 | b0 | 0  | 0  |
 896             +----++----+----+----+----+
 897             | b2 || b2 | 0  | b0 | 0  |
 898             +----++----+----+----+----+
 899             | b3 || b3 | 0  | 0  | b0 |
 900             +----++----+----+----+----+
 901
 902         """
 903         n = self.dimension()
 904         # Prepend the header row.
 905         M = [["*"] + list(self.gens())]
 906
 907         # And to each subsequent row, prepend an entry that belongs to
 908         # the left-side "header column."
 909         M += [ [self.monomial(i)] + [ self.monomial(i)*self.monomial(j)
 910                                     for j in range(n) ]
 911                for i in range(n) ]
 912
 913         return table(M, header_row=True, header_column=True, frame=True)
 914
 915
 916     def matrix_basis(self):
 917         """
 918         Return an (often more natural) representation of this algebras
 919         basis as an ordered tuple of matrices.
 920
 921         Every finite-dimensional Euclidean Jordan Algebra is a, up to
 922         Jordan isomorphism, a direct sum of five simple
 923         algebras---four of which comprise Hermitian matrices. And the
 924         last type of algebra can of course be thought of as `n`-by-`1`
 925         column matrices (ambiguusly called column vectors) to avoid
 926         special cases. As a result, matrices (and column vectors) are
 927         a natural representation format for Euclidean Jordan algebra
 928         elements.
 929
 930         But, when we construct an algebra from a basis of matrices,
 931         those matrix representations are lost in favor of coordinate
 932         vectors *with respect to* that basis. We could eventually
 933         convert back if we tried hard enough, but having the original
 934         representations handy is valuable enough that we simply store
 935         them and return them from this method.
 936
 937         Why implement this for non-matrix algebras? Avoiding special
 938         cases for the :class:`BilinearFormEJA` pays with simplicity in
 939         its own right. But mainly, we would like to be able to assume
 940         that elements of a :class:`CartesianProductEJA` can be displayed
 941         nicely, without having to have special classes for direct sums
 942         one of whose components was a matrix algebra.
 943
 944         SETUP::
 945
 946             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
 947             ....:                                  RealSymmetricEJA)
 948
 949         EXAMPLES::
 950
 951             sage: J = RealSymmetricEJA(2)
 952             sage: J.basis()
 953             Finite family {0: b0, 1: b1, 2: b2}
 954             sage: J.matrix_basis()
 955             (
 956             [1 0]  [                  0 0.7071067811865475?]  [0 0]
 957             [0 0], [0.7071067811865475?                   0], [0 1]
 958             )
 959
 960         ::
 961
 962             sage: J = JordanSpinEJA(2)
 963             sage: J.basis()
 964             Finite family {0: b0, 1: b1}
 965             sage: J.matrix_basis()
 966             (
 967             [1]  [0]
 968             [0], [1]
 969             )
 970         """
 971         return self._matrix_basis
 972
 973
 974     def matrix_space(self):
 975         """
 976         Return the matrix space in which this algebra's elements live, if
 977         we think of them as matrices (including column vectors of the
 978         appropriate size).
 979
 980         "By default" this will be an `n`-by-`1` column-matrix space,
 981         except when the algebra is trivial. There it's `n`-by-`n`
 982         (where `n` is zero), to ensure that two elements of the matrix
 983         space (empty matrices) can be multiplied. For algebras of
 984         matrices, this returns the space in which their
 985         real embeddings live.
 986
 987         SETUP::
 988
 989             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
 990             ....:                                  JordanSpinEJA,
 991             ....:                                  QuaternionHermitianEJA,
 992             ....:                                  TrivialEJA)
 993
 994         EXAMPLES:
 995
 996         By default, the matrix representation is just a column-matrix
 997         equivalent to the vector representation::
 998
 999             sage: J = JordanSpinEJA(3)
1000             sage: J.matrix_space()
1001             Full MatrixSpace of 3 by 1 dense matrices over Algebraic
1002             Real Field
1003
1004         The matrix representation in the trivial algebra is
1005         zero-by-zero instead of the usual `n`-by-one::
1006
1007             sage: J = TrivialEJA()
1008             sage: J.matrix_space()
1009             Full MatrixSpace of 0 by 0 dense matrices over Algebraic
1010             Real Field
1011
1012         The matrix space for complex/quaternion Hermitian matrix EJA
1013         is the space in which their real-embeddings live, not the
1014         original complex/quaternion matrix space::
1015
1016             sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1017             sage: J.matrix_space()
1018             Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1019             sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False)
1020             sage: J.matrix_space()
1021             Full MatrixSpace of 4 by 4 dense matrices over Rational Field
1022
1023         """
1024         if self.is_trivial():
1025             return MatrixSpace(self.base_ring(), 0)
1026         else:
1027             return self.matrix_basis()[0].parent()
1028
1029
1030     @cached_method
1031     def one(self):
1032         """
1033         Return the unit element of this algebra.
1034
1035         SETUP::
1036
1037             sage: from mjo.eja.eja_algebra import (HadamardEJA,
1038             ....:                                  random_eja)
1039
1040         EXAMPLES:
1041
1042         We can compute unit element in the Hadamard EJA::
1043
1044             sage: J = HadamardEJA(5)
1045             sage: J.one()
1046             b0 + b1 + b2 + b3 + b4
1047
1048         The unit element in the Hadamard EJA is inherited in the
1049         subalgebras generated by its elements::
1050
1051             sage: J = HadamardEJA(5)
1052             sage: J.one()
1053             b0 + b1 + b2 + b3 + b4
1054             sage: x = sum(J.gens())
1055             sage: A = x.subalgebra_generated_by(orthonormalize=False)
1056             sage: A.one()
1057             c0
1058             sage: A.one().superalgebra_element()
1059             b0 + b1 + b2 + b3 + b4
1060
1061         TESTS:
1062
1063         The identity element acts like the identity, regardless of
1064         whether or not we orthonormalize::
1065
1066             sage: set_random_seed()
1067             sage: J = random_eja()
1068             sage: x = J.random_element()
1069             sage: J.one()*x == x and x*J.one() == x
1070             True
1071             sage: A = x.subalgebra_generated_by()
1072             sage: y = A.random_element()
1073             sage: A.one()*y == y and y*A.one() == y
1074             True
1075
1076         ::
1077
1078             sage: set_random_seed()
1079             sage: J = random_eja(field=QQ, orthonormalize=False)
1080             sage: x = J.random_element()
1081             sage: J.one()*x == x and x*J.one() == x
1082             True
1083             sage: A = x.subalgebra_generated_by(orthonormalize=False)
1084             sage: y = A.random_element()
1085             sage: A.one()*y == y and y*A.one() == y
1086             True
1087
1088         The matrix of the unit element's operator is the identity,
1089         regardless of the base field and whether or not we
1090         orthonormalize::
1091
1092             sage: set_random_seed()
1093             sage: J = random_eja()
1094             sage: actual = J.one().operator().matrix()
1095             sage: expected = matrix.identity(J.base_ring(), J.dimension())
1096             sage: actual == expected
1097             True
1098             sage: x = J.random_element()
1099             sage: A = x.subalgebra_generated_by()
1100             sage: actual = A.one().operator().matrix()
1101             sage: expected = matrix.identity(A.base_ring(), A.dimension())
1102             sage: actual == expected
1103             True
1104
1105         ::
1106
1107             sage: set_random_seed()
1108             sage: J = random_eja(field=QQ, orthonormalize=False)
1109             sage: actual = J.one().operator().matrix()
1110             sage: expected = matrix.identity(J.base_ring(), J.dimension())
1111             sage: actual == expected
1112             True
1113             sage: x = J.random_element()
1114             sage: A = x.subalgebra_generated_by(orthonormalize=False)
1115             sage: actual = A.one().operator().matrix()
1116             sage: expected = matrix.identity(A.base_ring(), A.dimension())
1117             sage: actual == expected
1118             True
1119
1120         Ensure that the cached unit element (often precomputed by
1121         hand) agrees with the computed one::
1122
1123             sage: set_random_seed()
1124             sage: J = random_eja()
1125             sage: cached = J.one()
1126             sage: J.one.clear_cache()
1127             sage: J.one() == cached
1128             True
1129
1130         ::
1131
1132             sage: set_random_seed()
1133             sage: J = random_eja(field=QQ, orthonormalize=False)
1134             sage: cached = J.one()
1135             sage: J.one.clear_cache()
1136             sage: J.one() == cached
1137             True
1138
1139         """
1140         # We can brute-force compute the matrices of the operators
1141         # that correspond to the basis elements of this algebra.
1142         # If some linear combination of those basis elements is the
1143         # algebra identity, then the same linear combination of
1144         # their matrices has to be the identity matrix.
1145         #
1146         # Of course, matrices aren't vectors in sage, so we have to
1147         # appeal to the "long vectors" isometry.
1148         oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
1149
1150         # Now we use basic linear algebra to find the coefficients,
1151         # of the matrices-as-vectors-linear-combination, which should
1152         # work for the original algebra basis too.
1153         A = matrix(self.base_ring(), oper_vecs)
1154
1155         # We used the isometry on the left-hand side already, but we
1156         # still need to do it for the right-hand side. Recall that we
1157         # wanted something that summed to the identity matrix.
1158         b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
1159
1160         # Now if there's an identity element in the algebra, this
1161         # should work. We solve on the left to avoid having to
1162         # transpose the matrix "A".
1163         return self.from_vector(A.solve_left(b))
1164
1165
1166     def peirce_decomposition(self, c):
1167         """
1168         The Peirce decomposition of this algebra relative to the
1169         idempotent ``c``.
1170
1171         In the future, this can be extended to a complete system of
1172         orthogonal idempotents.
1173
1174         INPUT:
1175
1176           - ``c`` -- an idempotent of this algebra.
1177
1178         OUTPUT:
1179
1180         A triple (J0, J5, J1) containing two subalgebras and one subspace
1181         of this algebra,
1182
1183           - ``J0`` -- the algebra on the eigenspace of ``c.operator()``
1184             corresponding to the eigenvalue zero.
1185
1186           - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()``
1187             corresponding to the eigenvalue one-half.
1188
1189           - ``J1`` -- the algebra on the eigenspace of ``c.operator()``
1190             corresponding to the eigenvalue one.
1191
1192         These are the only possible eigenspaces for that operator, and this
1193         algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are
1194         orthogonal, and are subalgebras of this algebra with the appropriate
1195         restrictions.
1196
1197         SETUP::
1198
1199             sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA
1200
1201         EXAMPLES:
1202
1203         The canonical example comes from the symmetric matrices, which
1204         decompose into diagonal and off-diagonal parts::
1205
1206             sage: J = RealSymmetricEJA(3)
1207             sage: C = matrix(QQ, [ [1,0,0],
1208             ....:                  [0,1,0],
1209             ....:                  [0,0,0] ])
1210             sage: c = J(C)
1211             sage: J0,J5,J1 = J.peirce_decomposition(c)
1212             sage: J0
1213             Euclidean Jordan algebra of dimension 1...
1214             sage: J5
1215             Vector space of degree 6 and dimension 2...
1216             sage: J1
1217             Euclidean Jordan algebra of dimension 3...
1218             sage: J0.one().to_matrix()
1219             [0 0 0]
1220             [0 0 0]
1221             [0 0 1]
1222             sage: orig_df = AA.options.display_format
1223             sage: AA.options.display_format = 'radical'
1224             sage: J.from_vector(J5.basis()[0]).to_matrix()
1225             [          0           0 1/2*sqrt(2)]
1226             [          0           0           0]
1227             [1/2*sqrt(2)           0           0]
1228             sage: J.from_vector(J5.basis()[1]).to_matrix()
1229             [          0           0           0]
1230             [          0           0 1/2*sqrt(2)]
1231             [          0 1/2*sqrt(2)           0]
1232             sage: AA.options.display_format = orig_df
1233             sage: J1.one().to_matrix()
1234             [1 0 0]
1235             [0 1 0]
1236             [0 0 0]
1237
1238         TESTS:
1239
1240         Every algebra decomposes trivially with respect to its identity
1241         element::
1242
1243             sage: set_random_seed()
1244             sage: J = random_eja()
1245             sage: J0,J5,J1 = J.peirce_decomposition(J.one())
1246             sage: J0.dimension() == 0 and J5.dimension() == 0
1247             True
1248             sage: J1.superalgebra() == J and J1.dimension() == J.dimension()
1249             True
1250
1251         The decomposition is into eigenspaces, and its components are
1252         therefore necessarily orthogonal. Moreover, the identity
1253         elements in the two subalgebras are the projections onto their
1254         respective subspaces of the superalgebra's identity element::
1255
1256             sage: set_random_seed()
1257             sage: J = random_eja()
1258             sage: x = J.random_element()
1259             sage: if not J.is_trivial():
1260             ....:     while x.is_nilpotent():
1261             ....:         x = J.random_element()
1262             sage: c = x.subalgebra_idempotent()
1263             sage: J0,J5,J1 = J.peirce_decomposition(c)
1264             sage: ipsum = 0
1265             sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()):
1266             ....:     w = w.superalgebra_element()
1267             ....:     y = J.from_vector(y)
1268             ....:     z = z.superalgebra_element()
1269             ....:     ipsum += w.inner_product(y).abs()
1270             ....:     ipsum += w.inner_product(z).abs()
1271             ....:     ipsum += y.inner_product(z).abs()
1272             sage: ipsum
1273             0
1274             sage: J1(c) == J1.one()
1275             True
1276             sage: J0(J.one() - c) == J0.one()
1277             True
1278
1279         """
1280         if not c.is_idempotent():
1281             raise ValueError("element is not idempotent: %s" % c)
1282
1283         # Default these to what they should be if they turn out to be
1284         # trivial, because eigenspaces_left() won't return eigenvalues
1285         # corresponding to trivial spaces (e.g. it returns only the
1286         # eigenspace corresponding to lambda=1 if you take the
1287         # decomposition relative to the identity element).
1288         trivial = self.subalgebra(())
1289         J0 = trivial                          # eigenvalue zero
1290         J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
1291         J1 = trivial                          # eigenvalue one
1292
1293         for (eigval, eigspace) in c.operator().matrix().right_eigenspaces():
1294             if eigval == ~(self.base_ring()(2)):
1295                 J5 = eigspace
1296             else:
1297                 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
1298                 subalg = self.subalgebra(gens, check_axioms=False)
1299                 if eigval == 0:
1300                     J0 = subalg
1301                 elif eigval == 1:
1302                     J1 = subalg
1303                 else:
1304                     raise ValueError("unexpected eigenvalue: %s" % eigval)
1305
1306         return (J0, J5, J1)
1307
1308
1309     def random_element(self, thorough=False):
1310         r"""
1311         Return a random element of this algebra.
1312
1313         Our algebra superclass method only returns a linear
1314         combination of at most two basis elements. We instead
1315         want the vector space "random element" method that
1316         returns a more diverse selection.
1317
1318         INPUT:
1319
1320         - ``thorough`` -- (boolean; default False) whether or not we
1321           should generate irrational coefficients for the random
1322           element when our base ring is irrational; this slows the
1323           algebra operations to a crawl, but any truly random method
1324           should include them
1325
1326         """
1327         # For a general base ring... maybe we can trust this to do the
1328         # right thing? Unlikely, but.
1329         V = self.vector_space()
1330         v = V.random_element()
1331
1332         if self.base_ring() is AA:
1333             # The "random element" method of the algebraic reals is
1334             # stupid at the moment, and only returns integers between
1335             # -2 and 2, inclusive:
1336             #
1337             #   https://trac.sagemath.org/ticket/30875
1338             #
1339             # Instead, we implement our own "random vector" method,
1340             # and then coerce that into the algebra. We use the vector
1341             # space degree here instead of the dimension because a
1342             # subalgebra could (for example) be spanned by only two
1343             # vectors, each with five coordinates.  We need to
1344             # generate all five coordinates.
1345             if thorough:
1346                 v *= QQbar.random_element().real()
1347             else:
1348                 v *= QQ.random_element()
1349
1350         return self.from_vector(V.coordinate_vector(v))
1351
1352     def random_elements(self, count, thorough=False):
1353         """
1354         Return ``count`` random elements as a tuple.
1355
1356         INPUT:
1357
1358         - ``thorough`` -- (boolean; default False) whether or not we
1359           should generate irrational coefficients for the random
1360           elements when our base ring is irrational; this slows the
1361           algebra operations to a crawl, but any truly random method
1362           should include them
1363
1364         SETUP::
1365
1366             sage: from mjo.eja.eja_algebra import JordanSpinEJA
1367
1368         EXAMPLES::
1369
1370             sage: J = JordanSpinEJA(3)
1371             sage: x,y,z = J.random_elements(3)
1372             sage: all( [ x in J, y in J, z in J ])
1373             True
1374             sage: len( J.random_elements(10) ) == 10
1375             True
1376
1377         """
1378         return tuple( self.random_element(thorough)
1379                       for idx in range(count) )
1380
1381
1382     @cached_method
1383     def _charpoly_coefficients(self):
1384         r"""
1385         The `r` polynomial coefficients of the "characteristic polynomial
1386         of" function.
1387
1388         SETUP::
1389
1390             sage: from mjo.eja.eja_algebra import random_eja
1391
1392         TESTS:
1393
1394         The theory shows that these are all homogeneous polynomials of
1395         a known degree::
1396
1397             sage: set_random_seed()
1398             sage: J = random_eja()
1399             sage: all(p.is_homogeneous() for p in J._charpoly_coefficients())
1400             True
1401
1402         """
1403         n = self.dimension()
1404         R = self.coordinate_polynomial_ring()
1405         vars = R.gens()
1406         F = R.fraction_field()
1407
1408         def L_x_i_j(i,j):
1409             # From a result in my book, these are the entries of the
1410             # basis representation of L_x.
1411             return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
1412                         for k in range(n) )
1413
1414         L_x = matrix(F, n, n, L_x_i_j)
1415
1416         r = None
1417         if self.rank.is_in_cache():
1418             r = self.rank()
1419             # There's no need to pad the system with redundant
1420             # columns if we *know* they'll be redundant.
1421             n = r
1422
1423         # Compute an extra power in case the rank is equal to
1424         # the dimension (otherwise, we would stop at x^(r-1)).
1425         x_powers = [ (L_x**k)*self.one().to_vector()
1426                      for k in range(n+1) ]
1427         A = matrix.column(F, x_powers[:n])
1428         AE = A.extended_echelon_form()
1429         E = AE[:,n:]
1430         A_rref = AE[:,:n]
1431         if r is None:
1432             r = A_rref.rank()
1433         b = x_powers[r]
1434
1435         # The theory says that only the first "r" coefficients are
1436         # nonzero, and they actually live in the original polynomial
1437         # ring and not the fraction field. We negate them because in
1438         # the actual characteristic polynomial, they get moved to the
1439         # other side where x^r lives. We don't bother to trim A_rref
1440         # down to a square matrix and solve the resulting system,
1441         # because the upper-left r-by-r portion of A_rref is
1442         # guaranteed to be the identity matrix, so e.g.
1443         #
1444         #   A_rref.solve_right(Y)
1445         #
1446         # would just be returning Y.
1447         return (-E*b)[:r].change_ring(R)
1448
1449     @cached_method
1450     def rank(self):
1451         r"""
1452         Return the rank of this EJA.
1453
1454         This is a cached method because we know the rank a priori for
1455         all of the algebras we can construct. Thus we can avoid the
1456         expensive ``_charpoly_coefficients()`` call unless we truly
1457         need to compute the whole characteristic polynomial.
1458
1459         SETUP::
1460
1461             sage: from mjo.eja.eja_algebra import (HadamardEJA,
1462             ....:                                  JordanSpinEJA,
1463             ....:                                  RealSymmetricEJA,
1464             ....:                                  ComplexHermitianEJA,
1465             ....:                                  QuaternionHermitianEJA,
1466             ....:                                  random_eja)
1467
1468         EXAMPLES:
1469
1470         The rank of the Jordan spin algebra is always two::
1471
1472             sage: JordanSpinEJA(2).rank()
1473             2
1474             sage: JordanSpinEJA(3).rank()
1475             2
1476             sage: JordanSpinEJA(4).rank()
1477             2
1478
1479         The rank of the `n`-by-`n` Hermitian real, complex, or
1480         quaternion matrices is `n`::
1481
1482             sage: RealSymmetricEJA(4).rank()
1483             4
1484             sage: ComplexHermitianEJA(3).rank()
1485             3
1486             sage: QuaternionHermitianEJA(2).rank()
1487             2
1488
1489         TESTS:
1490
1491         Ensure that every EJA that we know how to construct has a
1492         positive integer rank, unless the algebra is trivial in
1493         which case its rank will be zero::
1494
1495             sage: set_random_seed()
1496             sage: J = random_eja()
1497             sage: r = J.rank()
1498             sage: r in ZZ
1499             True
1500             sage: r > 0 or (r == 0 and J.is_trivial())
1501             True
1502
1503         Ensure that computing the rank actually works, since the ranks
1504         of all simple algebras are known and will be cached by default::
1505
1506             sage: set_random_seed()    # long time
1507             sage: J = random_eja()     # long time
1508             sage: cached = J.rank()    # long time
1509             sage: J.rank.clear_cache() # long time
1510             sage: J.rank() == cached   # long time
1511             True
1512
1513         """
1514         return len(self._charpoly_coefficients())
1515
1516
1517     def subalgebra(self, basis, **kwargs):
1518         r"""
1519         Create a subalgebra of this algebra from the given basis.
1520         """
1521         from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
1522         return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
1523
1524
1525     def vector_space(self):
1526         """
1527         Return the vector space that underlies this algebra.
1528
1529         SETUP::
1530
1531             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1532
1533         EXAMPLES::
1534
1535             sage: J = RealSymmetricEJA(2)
1536             sage: J.vector_space()
1537             Vector space of dimension 3 over...
1538
1539         """
1540         return self.zero().to_vector().parent().ambient_vector_space()
1541
1542
1543
1544 class RationalBasisEJA(FiniteDimensionalEJA):
1545     r"""
1546     New class for algebras whose supplied basis elements have all rational entries.
1547
1548     SETUP::
1549
1550         sage: from mjo.eja.eja_algebra import BilinearFormEJA
1551
1552     EXAMPLES:
1553
1554     The supplied basis is orthonormalized by default::
1555
1556         sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]])
1557         sage: J = BilinearFormEJA(B)
1558         sage: J.matrix_basis()
1559         (
1560         [1]  [  0]  [   0]
1561         [0]  [1/5]  [32/5]
1562         [0], [  0], [   5]
1563         )
1564
1565     """
1566     def __init__(self,
1567                  basis,
1568                  jordan_product,
1569                  inner_product,
1570                  field=AA,
1571                  check_field=True,
1572                  **kwargs):
1573
1574         if check_field:
1575             # Abuse the check_field parameter to check that the entries of
1576             # out basis (in ambient coordinates) are in the field QQ.
1577             if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
1578                 raise TypeError("basis not rational")
1579
1580         super().__init__(basis,
1581                          jordan_product,
1582                          inner_product,
1583                          field=field,
1584                          check_field=check_field,
1585                          **kwargs)
1586
1587         self._rational_algebra = None
1588         if field is not QQ:
1589             # There's no point in constructing the extra algebra if this
1590             # one is already rational.
1591             #
1592             # Note: the same Jordan and inner-products work here,
1593             # because they are necessarily defined with respect to
1594             # ambient coordinates and not any particular basis.
1595             self._rational_algebra = FiniteDimensionalEJA(
1596                                        basis,
1597                                        jordan_product,
1598                                        inner_product,
1599                                        field=QQ,
1600                                        associative=self.is_associative(),
1601                                        orthonormalize=False,
1602                                        check_field=False,
1603                                        check_axioms=False)
1604
1605     @cached_method
1606     def _charpoly_coefficients(self):
1607         r"""
1608         SETUP::
1609
1610             sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
1611             ....:                                  JordanSpinEJA)
1612
1613         EXAMPLES:
1614
1615         The base ring of the resulting polynomial coefficients is what
1616         it should be, and not the rationals (unless the algebra was
1617         already over the rationals)::
1618
1619             sage: J = JordanSpinEJA(3)
1620             sage: J._charpoly_coefficients()
1621             (X1^2 - X2^2 - X3^2, -2*X1)
1622             sage: a0 = J._charpoly_coefficients()[0]
1623             sage: J.base_ring()
1624             Algebraic Real Field
1625             sage: a0.base_ring()
1626             Algebraic Real Field
1627
1628         """
1629         if self._rational_algebra is None:
1630             # There's no need to construct *another* algebra over the
1631             # rationals if this one is already over the
1632             # rationals. Likewise, if we never orthonormalized our
1633             # basis, we might as well just use the given one.
1634             return super()._charpoly_coefficients()
1635
1636         # Do the computation over the rationals. The answer will be
1637         # the same, because all we've done is a change of basis.
1638         # Then, change back from QQ to our real base ring
1639         a = ( a_i.change_ring(self.base_ring())
1640               for a_i in self._rational_algebra._charpoly_coefficients() )
1641
1642         if self._deortho_matrix is None:
1643             # This can happen if our base ring was, say, AA and we
1644             # chose not to (or didn't need to) orthonormalize. It's
1645             # still faster to do the computations over QQ even if
1646             # the numbers in the boxes stay the same.
1647             return tuple(a)
1648
1649         # Otherwise, convert the coordinate variables back to the
1650         # deorthonormalized ones.
1651         R = self.coordinate_polynomial_ring()
1652         from sage.modules.free_module_element import vector
1653         X = vector(R, R.gens())
1654         BX = self._deortho_matrix*X
1655
1656         subs_dict = { X[i]: BX[i] for i in range(len(X)) }
1657         return tuple( a_i.subs(subs_dict) for a_i in a )
1658
1659 class ConcreteEJA(RationalBasisEJA):
1660     r"""
1661     A class for the Euclidean Jordan algebras that we know by name.
1662
1663     These are the Jordan algebras whose basis, multiplication table,
1664     rank, and so on are known a priori. More to the point, they are
1665     the Euclidean Jordan algebras for which we are able to conjure up
1666     a "random instance."
1667
1668     SETUP::
1669
1670         sage: from mjo.eja.eja_algebra import ConcreteEJA
1671
1672     TESTS:
1673
1674     Our basis is normalized with respect to the algebra's inner
1675     product, unless we specify otherwise::
1676
1677         sage: set_random_seed()
1678         sage: J = ConcreteEJA.random_instance()
1679         sage: all( b.norm() == 1 for b in J.gens() )
1680         True
1681
1682     Since our basis is orthonormal with respect to the algebra's inner
1683     product, and since we know that this algebra is an EJA, any
1684     left-multiplication operator's matrix will be symmetric because
1685     natural->EJA basis representation is an isometry and within the
1686     EJA the operator is self-adjoint by the Jordan axiom::
1687
1688         sage: set_random_seed()
1689         sage: J = ConcreteEJA.random_instance()
1690         sage: x = J.random_element()
1691         sage: x.operator().is_self_adjoint()
1692         True
1693     """
1694
1695     @staticmethod
1696     def _max_random_instance_size():
1697         """
1698         Return an integer "size" that is an upper bound on the size of
1699         this algebra when it is used in a random test
1700         case. Unfortunately, the term "size" is ambiguous -- when
1701         dealing with `R^n` under either the Hadamard or Jordan spin
1702         product, the "size" refers to the dimension `n`. When dealing
1703         with a matrix algebra (real symmetric or complex/quaternion
1704         Hermitian), it refers to the size of the matrix, which is far
1705         less than the dimension of the underlying vector space.
1706
1707         This method must be implemented in each subclass.
1708         """
1709         raise NotImplementedError
1710
1711     @classmethod
1712     def random_instance(cls, *args, **kwargs):
1713         """
1714         Return a random instance of this type of algebra.
1715
1716         This method should be implemented in each subclass.
1717         """
1718         from sage.misc.prandom import choice
1719         eja_class = choice(cls.__subclasses__())
1720
1721         # These all bubble up to the RationalBasisEJA superclass
1722         # constructor, so any (kw)args valid there are also valid
1723         # here.
1724         return eja_class.random_instance(*args, **kwargs)
1725
1726
1727 class MatrixEJA:
1728     @staticmethod
1729     def jordan_product(X,Y):
1730         return (X*Y + Y*X)/2
1731
1732     @staticmethod
1733     def trace_inner_product(X,Y):
1734         r"""
1735         A trace inner-product for matrices that aren't embedded in the
1736         reals. It takes MATRICES as arguments, not EJA elements.
1737         """
1738         return (X*Y).trace().real()
1739
1740 class RealEmbeddedMatrixEJA(MatrixEJA):
1741     @staticmethod
1742     def dimension_over_reals():
1743         r"""
1744         The dimension of this matrix's base ring over the reals.
1745
1746         The reals are dimension one over themselves, obviously; that's
1747         just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
1748         have dimension two. Finally, the quaternions have dimension
1749         four over the reals.
1750
1751         This is used to determine the size of the matrix returned from
1752         :meth:`real_embed`, among other things.
1753         """
1754         raise NotImplementedError
1755
1756     @classmethod
1757     def real_embed(cls,M):
1758         """
1759         Embed the matrix ``M`` into a space of real matrices.
1760
1761         The matrix ``M`` can have entries in any field at the moment:
1762         the real numbers, complex numbers, or quaternions. And although
1763         they are not a field, we can probably support octonions at some
1764         point, too. This function returns a real matrix that "acts like"
1765         the original with respect to matrix multiplication; i.e.
1766
1767           real_embed(M*N) = real_embed(M)*real_embed(N)
1768
1769         """
1770         if M.ncols() != M.nrows():
1771             raise ValueError("the matrix 'M' must be square")
1772         return M
1773
1774
1775     @classmethod
1776     def real_unembed(cls,M):
1777         """
1778         The inverse of :meth:`real_embed`.
1779         """
1780         if M.ncols() != M.nrows():
1781             raise ValueError("the matrix 'M' must be square")
1782         if not ZZ(M.nrows()).mod(cls.dimension_over_reals()).is_zero():
1783             raise ValueError("the matrix 'M' must be a real embedding")
1784         return M
1785
1786
1787     @classmethod
1788     def trace_inner_product(cls,X,Y):
1789         r"""
1790         Compute the trace inner-product of two real-embeddings.
1791
1792         SETUP::
1793
1794             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
1795             ....:                                  QuaternionHermitianEJA)
1796
1797         EXAMPLES::
1798
1799             sage: set_random_seed()
1800             sage: J = ComplexHermitianEJA.random_instance()
1801             sage: x,y = J.random_elements(2)
1802             sage: Xe = x.to_matrix()
1803             sage: Ye = y.to_matrix()
1804             sage: X = J.real_unembed(Xe)
1805             sage: Y = J.real_unembed(Ye)
1806             sage: expected = (X*Y).trace().real()
1807             sage: actual = J.trace_inner_product(Xe,Ye)
1808             sage: actual == expected
1809             True
1810
1811         ::
1812
1813             sage: set_random_seed()
1814             sage: J = QuaternionHermitianEJA.random_instance()
1815             sage: x,y = J.random_elements(2)
1816             sage: Xe = x.to_matrix()
1817             sage: Ye = y.to_matrix()
1818             sage: X = J.real_unembed(Xe)
1819             sage: Y = J.real_unembed(Ye)
1820             sage: expected = (X*Y).trace().coefficient_tuple()[0]
1821             sage: actual = J.trace_inner_product(Xe,Ye)
1822             sage: actual == expected
1823             True
1824
1825         """
1826         # This does in fact compute the real part of the trace.
1827         # If we compute the trace of e.g. a complex matrix M,
1828         # then we do so by adding up its diagonal entries --
1829         # call them z_1 through z_n. The real embedding of z_1
1830         # will be a 2-by-2 REAL matrix [a, b; -b, a] whose trace
1831         # as a REAL matrix will be 2*a = 2*Re(z_1). And so forth.
1832         return (X*Y).trace()/cls.dimension_over_reals()
1833
1834 class RealSymmetricEJA(ConcreteEJA, MatrixEJA):
1835     """
1836     The rank-n simple EJA consisting of real symmetric n-by-n
1837     matrices, the usual symmetric Jordan product, and the trace inner
1838     product. It has dimension `(n^2 + n)/2` over the reals.
1839
1840     SETUP::
1841
1842         sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1843
1844     EXAMPLES::
1845
1846         sage: J = RealSymmetricEJA(2)
1847         sage: b0, b1, b2 = J.gens()
1848         sage: b0*b0
1849         b0
1850         sage: b1*b1
1851         1/2*b0 + 1/2*b2
1852         sage: b2*b2
1853         b2
1854
1855     In theory, our "field" can be any subfield of the reals::
1856
1857         sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1858         Euclidean Jordan algebra of dimension 3 over Real Double Field
1859         sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1860         Euclidean Jordan algebra of dimension 3 over Real Field with
1861         53 bits of precision
1862
1863     TESTS:
1864
1865     The dimension of this algebra is `(n^2 + n) / 2`::
1866
1867         sage: set_random_seed()
1868         sage: n_max = RealSymmetricEJA._max_random_instance_size()
1869         sage: n = ZZ.random_element(1, n_max)
1870         sage: J = RealSymmetricEJA(n)
1871         sage: J.dimension() == (n^2 + n)/2
1872         True
1873
1874     The Jordan multiplication is what we think it is::
1875
1876         sage: set_random_seed()
1877         sage: J = RealSymmetricEJA.random_instance()
1878         sage: x,y = J.random_elements(2)
1879         sage: actual = (x*y).to_matrix()
1880         sage: X = x.to_matrix()
1881         sage: Y = y.to_matrix()
1882         sage: expected = (X*Y + Y*X)/2
1883         sage: actual == expected
1884         True
1885         sage: J(expected) == x*y
1886         True
1887
1888     We can change the generator prefix::
1889
1890         sage: RealSymmetricEJA(3, prefix='q').gens()
1891         (q0, q1, q2, q3, q4, q5)
1892
1893     We can construct the (trivial) algebra of rank zero::
1894
1895         sage: RealSymmetricEJA(0)
1896         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1897
1898     """
1899     @classmethod
1900     def _denormalized_basis(cls, n, field):
1901         """
1902         Return a basis for the space of real symmetric n-by-n matrices.
1903
1904         SETUP::
1905
1906             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1907
1908         TESTS::
1909
1910             sage: set_random_seed()
1911             sage: n = ZZ.random_element(1,5)
1912             sage: B = RealSymmetricEJA._denormalized_basis(n,ZZ)
1913             sage: all( M.is_symmetric() for M in  B)
1914             True
1915
1916         """
1917         # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1918         # coordinates.
1919         S = []
1920         for i in range(n):
1921             for j in range(i+1):
1922                 Eij = matrix(field, n, lambda k,l: k==i and l==j)
1923                 if i == j:
1924                     Sij = Eij
1925                 else:
1926                     Sij = Eij + Eij.transpose()
1927                 S.append(Sij)
1928         return tuple(S)
1929
1930
1931     @staticmethod
1932     def _max_random_instance_size():
1933         return 4 # Dimension 10
1934
1935     @classmethod
1936     def random_instance(cls, **kwargs):
1937         """
1938         Return a random instance of this type of algebra.
1939         """
1940         n = ZZ.random_element(cls._max_random_instance_size() + 1)
1941         return cls(n, **kwargs)
1942
1943     def __init__(self, n, field=AA, **kwargs):
1944         # We know this is a valid EJA, but will double-check
1945         # if the user passes check_axioms=True.
1946         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
1947
1948         associative = False
1949         if n <= 1:
1950             associative = True
1951
1952         super().__init__(self._denormalized_basis(n,field),
1953                          self.jordan_product,
1954                          self.trace_inner_product,
1955                          field=field,
1956                          associative=associative,
1957                          **kwargs)
1958
1959         # TODO: this could be factored out somehow, but is left here
1960         # because the MatrixEJA is not presently a subclass of the
1961         # FDEJA class that defines rank() and one().
1962         self.rank.set_cache(n)
1963         idV = self.matrix_space().one()
1964         self.one.set_cache(self(idV))
1965
1966
1967
1968 class ComplexMatrixEJA(RealEmbeddedMatrixEJA):
1969     # A manual dictionary-cache for the complex_extension() method,
1970     # since apparently @classmethods can't also be @cached_methods.
1971     _complex_extension = {}
1972
1973     @classmethod
1974     def complex_extension(cls,field):
1975         r"""
1976         The complex field that we embed/unembed, as an extension
1977         of the given ``field``.
1978         """
1979         if field in cls._complex_extension:
1980             return cls._complex_extension[field]
1981
1982         # Sage doesn't know how to adjoin the complex "i" (the root of
1983         # x^2 + 1) to a field in a general way. Here, we just enumerate
1984         # all of the cases that I have cared to support so far.
1985         if field is AA:
1986             # Sage doesn't know how to embed AA into QQbar, i.e. how
1987             # to adjoin sqrt(-1) to AA.
1988             F = QQbar
1989         elif not field.is_exact():
1990             # RDF or RR
1991             F = field.complex_field()
1992         else:
1993             # Works for QQ and... maybe some other fields.
1994             R = PolynomialRing(field, 'z')
1995             z = R.gen()
1996             F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
1997
1998         cls._complex_extension[field] = F
1999         return F
2000
2001     @staticmethod
2002     def dimension_over_reals():
2003         return 2
2004
2005     @classmethod
2006     def real_embed(cls,M):
2007         """
2008         Embed the n-by-n complex matrix ``M`` into the space of real
2009         matrices of size 2n-by-2n via the map the sends each entry `z = a +
2010         bi` to the block matrix ``[[a,b],[-b,a]]``.
2011
2012         SETUP::
2013
2014             sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2015
2016         EXAMPLES::
2017
2018             sage: F = QuadraticField(-1, 'I')
2019             sage: x1 = F(4 - 2*i)
2020             sage: x2 = F(1 + 2*i)
2021             sage: x3 = F(-i)
2022             sage: x4 = F(6)
2023             sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
2024             sage: ComplexMatrixEJA.real_embed(M)
2025             [ 4 -2| 1  2]
2026             [ 2  4|-2  1]
2027             [-----+-----]
2028             [ 0 -1| 6  0]
2029             [ 1  0| 0  6]
2030
2031         TESTS:
2032
2033         Embedding is a homomorphism (isomorphism, in fact)::
2034
2035             sage: set_random_seed()
2036             sage: n = ZZ.random_element(3)
2037             sage: F = QuadraticField(-1, 'I')
2038             sage: X = random_matrix(F, n)
2039             sage: Y = random_matrix(F, n)
2040             sage: Xe = ComplexMatrixEJA.real_embed(X)
2041             sage: Ye = ComplexMatrixEJA.real_embed(Y)
2042             sage: XYe = ComplexMatrixEJA.real_embed(X*Y)
2043             sage: Xe*Ye == XYe
2044             True
2045
2046         """
2047         super().real_embed(M)
2048         n = M.nrows()
2049
2050         # We don't need any adjoined elements...
2051         field = M.base_ring().base_ring()
2052
2053         blocks = []
2054         for z in M.list():
2055             a = z.real()
2056             b = z.imag()
2057             blocks.append(matrix(field, 2, [ [ a, b],
2058                                              [-b, a] ]))
2059
2060         return matrix.block(field, n, blocks)
2061
2062
2063     @classmethod
2064     def real_unembed(cls,M):
2065         """
2066         The inverse of _embed_complex_matrix().
2067
2068         SETUP::
2069
2070             sage: from mjo.eja.eja_algebra import ComplexMatrixEJA
2071
2072         EXAMPLES::
2073
2074             sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
2075             ....:                 [-2,  1,  -4,  3],
2076             ....:                 [ 9,  10, 11, 12],
2077             ....:                 [-10, 9, -12, 11] ])
2078             sage: ComplexMatrixEJA.real_unembed(A)
2079             [  2*I + 1   4*I + 3]
2080             [ 10*I + 9 12*I + 11]
2081
2082         TESTS:
2083
2084         Unembedding is the inverse of embedding::
2085
2086             sage: set_random_seed()
2087             sage: F = QuadraticField(-1, 'I')
2088             sage: M = random_matrix(F, 3)
2089             sage: Me = ComplexMatrixEJA.real_embed(M)
2090             sage: ComplexMatrixEJA.real_unembed(Me) == M
2091             True
2092
2093         """
2094         super().real_unembed(M)
2095         n = ZZ(M.nrows())
2096         d = cls.dimension_over_reals()
2097         F = cls.complex_extension(M.base_ring())
2098         i = F.gen()
2099
2100         # Go top-left to bottom-right (reading order), converting every
2101         # 2-by-2 block we see to a single complex element.
2102         elements = []
2103         for k in range(n/d):
2104             for j in range(n/d):
2105                 submat = M[d*k:d*k+d,d*j:d*j+d]
2106                 if submat[0,0] != submat[1,1]:
2107                     raise ValueError('bad on-diagonal submatrix')
2108                 if submat[0,1] != -submat[1,0]:
2109                     raise ValueError('bad off-diagonal submatrix')
2110                 z = submat[0,0] + submat[0,1]*i
2111                 elements.append(z)
2112
2113         return matrix(F, n/d, elements)
2114
2115
2116 class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
2117     """
2118     The rank-n simple EJA consisting of complex Hermitian n-by-n
2119     matrices over the real numbers, the usual symmetric Jordan product,
2120     and the real-part-of-trace inner product. It has dimension `n^2` over
2121     the reals.
2122
2123     SETUP::
2124
2125         sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2126
2127     EXAMPLES:
2128
2129     In theory, our "field" can be any subfield of the reals::
2130
2131         sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
2132         Euclidean Jordan algebra of dimension 4 over Real Double Field
2133         sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
2134         Euclidean Jordan algebra of dimension 4 over Real Field with
2135         53 bits of precision
2136
2137     TESTS:
2138
2139     The dimension of this algebra is `n^2`::
2140
2141         sage: set_random_seed()
2142         sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2143         sage: n = ZZ.random_element(1, n_max)
2144         sage: J = ComplexHermitianEJA(n)
2145         sage: J.dimension() == n^2
2146         True
2147
2148     The Jordan multiplication is what we think it is::
2149
2150         sage: set_random_seed()
2151         sage: J = ComplexHermitianEJA.random_instance()
2152         sage: x,y = J.random_elements(2)
2153         sage: actual = (x*y).to_matrix()
2154         sage: X = x.to_matrix()
2155         sage: Y = y.to_matrix()
2156         sage: expected = (X*Y + Y*X)/2
2157         sage: actual == expected
2158         True
2159         sage: J(expected) == x*y
2160         True
2161
2162     We can change the generator prefix::
2163
2164         sage: ComplexHermitianEJA(2, prefix='z').gens()
2165         (z0, z1, z2, z3)
2166
2167     We can construct the (trivial) algebra of rank zero::
2168
2169         sage: ComplexHermitianEJA(0)
2170         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2171
2172     """
2173
2174     @classmethod
2175     def _denormalized_basis(cls, n, field):
2176         """
2177         Returns a basis for the space of complex Hermitian n-by-n matrices.
2178
2179         Why do we embed these? Basically, because all of numerical linear
2180         algebra assumes that you're working with vectors consisting of `n`
2181         entries from a field and scalars from the same field. There's no way
2182         to tell SageMath that (for example) the vectors contain complex
2183         numbers, while the scalar field is real.
2184
2185         SETUP::
2186
2187             sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
2188
2189         TESTS::
2190
2191             sage: set_random_seed()
2192             sage: n = ZZ.random_element(1,5)
2193             sage: B = ComplexHermitianEJA._denormalized_basis(n,ZZ)
2194             sage: all( M.is_symmetric() for M in  B)
2195             True
2196
2197         """
2198         R = PolynomialRing(ZZ, 'z')
2199         z = R.gen()
2200         F = ZZ.extension(z**2 + 1, 'I')
2201         I = F.gen(1)
2202
2203         # This is like the symmetric case, but we need to be careful:
2204         #
2205         #   * We want conjugate-symmetry, not just symmetry.
2206         #   * The diagonal will (as a result) be real.
2207         #
2208         S = []
2209         Eij = matrix.zero(F,n)
2210         for i in range(n):
2211             for j in range(i+1):
2212                 # "build" E_ij
2213                 Eij[i,j] = 1
2214                 if i == j:
2215                     Sij = cls.real_embed(Eij)
2216                     S.append(Sij)
2217                 else:
2218                     # The second one has a minus because it's conjugated.
2219                     Eij[j,i] = 1 # Eij = Eij + Eij.transpose()
2220                     Sij_real = cls.real_embed(Eij)
2221                     S.append(Sij_real)
2222                     # Eij = I*Eij - I*Eij.transpose()
2223                     Eij[i,j] = I
2224                     Eij[j,i] = -I
2225                     Sij_imag = cls.real_embed(Eij)
2226                     S.append(Sij_imag)
2227                     Eij[j,i] = 0
2228                 # "erase" E_ij
2229                 Eij[i,j] = 0
2230
2231         # Since we embedded the entries, we can drop back to the
2232         # desired real "field" instead of the extension "F".
2233         return tuple( s.change_ring(field) for s in S )
2234
2235
2236     def __init__(self, n, field=AA, **kwargs):
2237         # We know this is a valid EJA, but will double-check
2238         # if the user passes check_axioms=True.
2239         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2240
2241         associative = False
2242         if n <= 1:
2243             associative = True
2244
2245         super().__init__(self._denormalized_basis(n,field),
2246                          self.jordan_product,
2247                          self.trace_inner_product,
2248                          field=field,
2249                          associative=associative,
2250                          **kwargs)
2251         # TODO: this could be factored out somehow, but is left here
2252         # because the MatrixEJA is not presently a subclass of the
2253         # FDEJA class that defines rank() and one().
2254         self.rank.set_cache(n)
2255         idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
2256         self.one.set_cache(self(idV))
2257
2258     @staticmethod
2259     def _max_random_instance_size():
2260         return 3 # Dimension 9
2261
2262     @classmethod
2263     def random_instance(cls, **kwargs):
2264         """
2265         Return a random instance of this type of algebra.
2266         """
2267         n = ZZ.random_element(cls._max_random_instance_size() + 1)
2268         return cls(n, **kwargs)
2269
2270 class QuaternionMatrixEJA(RealEmbeddedMatrixEJA):
2271
2272     # A manual dictionary-cache for the quaternion_extension() method,
2273     # since apparently @classmethods can't also be @cached_methods.
2274     _quaternion_extension = {}
2275
2276     @classmethod
2277     def quaternion_extension(cls,field):
2278         r"""
2279         The quaternion field that we embed/unembed, as an extension
2280         of the given ``field``.
2281         """
2282         if field in cls._quaternion_extension:
2283             return cls._quaternion_extension[field]
2284
2285         Q = QuaternionAlgebra(field,-1,-1)
2286
2287         cls._quaternion_extension[field] = Q
2288         return Q
2289
2290     @staticmethod
2291     def dimension_over_reals():
2292         return 4
2293
2294     @classmethod
2295     def real_embed(cls,M):
2296         """
2297         Embed the n-by-n quaternion matrix ``M`` into the space of real
2298         matrices of size 4n-by-4n by first sending each quaternion entry `z
2299         = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
2300         c+di],[-c + di, a-bi]]`, and then embedding those into a real
2301         matrix.
2302
2303         SETUP::
2304
2305             sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2306
2307         EXAMPLES::
2308
2309             sage: Q = QuaternionAlgebra(QQ,-1,-1)
2310             sage: i,j,k = Q.gens()
2311             sage: x = 1 + 2*i + 3*j + 4*k
2312             sage: M = matrix(Q, 1, [[x]])
2313             sage: QuaternionMatrixEJA.real_embed(M)
2314             [ 1  2  3  4]
2315             [-2  1 -4  3]
2316             [-3  4  1 -2]
2317             [-4 -3  2  1]
2318
2319         Embedding is a homomorphism (isomorphism, in fact)::
2320
2321             sage: set_random_seed()
2322             sage: n = ZZ.random_element(2)
2323             sage: Q = QuaternionAlgebra(QQ,-1,-1)
2324             sage: X = random_matrix(Q, n)
2325             sage: Y = random_matrix(Q, n)
2326             sage: Xe = QuaternionMatrixEJA.real_embed(X)
2327             sage: Ye = QuaternionMatrixEJA.real_embed(Y)
2328             sage: XYe = QuaternionMatrixEJA.real_embed(X*Y)
2329             sage: Xe*Ye == XYe
2330             True
2331
2332         """
2333         super().real_embed(M)
2334         quaternions = M.base_ring()
2335         n = M.nrows()
2336
2337         F = QuadraticField(-1, 'I')
2338         i = F.gen()
2339
2340         blocks = []
2341         for z in M.list():
2342             t = z.coefficient_tuple()
2343             a = t[0]
2344             b = t[1]
2345             c = t[2]
2346             d = t[3]
2347             cplxM = matrix(F, 2, [[ a + b*i, c + d*i],
2348                                  [-c + d*i, a - b*i]])
2349             realM = ComplexMatrixEJA.real_embed(cplxM)
2350             blocks.append(realM)
2351
2352         # We should have real entries by now, so use the realest field
2353         # we've got for the return value.
2354         return matrix.block(quaternions.base_ring(), n, blocks)
2355
2356
2357
2358     @classmethod
2359     def real_unembed(cls,M):
2360         """
2361         The inverse of _embed_quaternion_matrix().
2362
2363         SETUP::
2364
2365             sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA
2366
2367         EXAMPLES::
2368
2369             sage: M = matrix(QQ, [[ 1,  2,  3,  4],
2370             ....:                 [-2,  1, -4,  3],
2371             ....:                 [-3,  4,  1, -2],
2372             ....:                 [-4, -3,  2,  1]])
2373             sage: QuaternionMatrixEJA.real_unembed(M)
2374             [1 + 2*i + 3*j + 4*k]
2375
2376         TESTS:
2377
2378         Unembedding is the inverse of embedding::
2379
2380             sage: set_random_seed()
2381             sage: Q = QuaternionAlgebra(QQ, -1, -1)
2382             sage: M = random_matrix(Q, 3)
2383             sage: Me = QuaternionMatrixEJA.real_embed(M)
2384             sage: QuaternionMatrixEJA.real_unembed(Me) == M
2385             True
2386
2387         """
2388         super().real_unembed(M)
2389         n = ZZ(M.nrows())
2390         d = cls.dimension_over_reals()
2391
2392         # Use the base ring of the matrix to ensure that its entries can be
2393         # multiplied by elements of the quaternion algebra.
2394         Q = cls.quaternion_extension(M.base_ring())
2395         i,j,k = Q.gens()
2396
2397         # Go top-left to bottom-right (reading order), converting every
2398         # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
2399         # quaternion block.
2400         elements = []
2401         for l in range(n/d):
2402             for m in range(n/d):
2403                 submat = ComplexMatrixEJA.real_unembed(
2404                     M[d*l:d*l+d,d*m:d*m+d] )
2405                 if submat[0,0] != submat[1,1].conjugate():
2406                     raise ValueError('bad on-diagonal submatrix')
2407                 if submat[0,1] != -submat[1,0].conjugate():
2408                     raise ValueError('bad off-diagonal submatrix')
2409                 z  = submat[0,0].real()
2410                 z += submat[0,0].imag()*i
2411                 z += submat[0,1].real()*j
2412                 z += submat[0,1].imag()*k
2413                 elements.append(z)
2414
2415         return matrix(Q, n/d, elements)
2416
2417
2418 class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
2419     r"""
2420     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2421     matrices, the usual symmetric Jordan product, and the
2422     real-part-of-trace inner product. It has dimension `2n^2 - n` over
2423     the reals.
2424
2425     SETUP::
2426
2427         sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2428
2429     EXAMPLES:
2430
2431     In theory, our "field" can be any subfield of the reals::
2432
2433         sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2434         Euclidean Jordan algebra of dimension 6 over Real Double Field
2435         sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2436         Euclidean Jordan algebra of dimension 6 over Real Field with
2437         53 bits of precision
2438
2439     TESTS:
2440
2441     The dimension of this algebra is `2*n^2 - n`::
2442
2443         sage: set_random_seed()
2444         sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2445         sage: n = ZZ.random_element(1, n_max)
2446         sage: J = QuaternionHermitianEJA(n)
2447         sage: J.dimension() == 2*(n^2) - n
2448         True
2449
2450     The Jordan multiplication is what we think it is::
2451
2452         sage: set_random_seed()
2453         sage: J = QuaternionHermitianEJA.random_instance()
2454         sage: x,y = J.random_elements(2)
2455         sage: actual = (x*y).to_matrix()
2456         sage: X = x.to_matrix()
2457         sage: Y = y.to_matrix()
2458         sage: expected = (X*Y + Y*X)/2
2459         sage: actual == expected
2460         True
2461         sage: J(expected) == x*y
2462         True
2463
2464     We can change the generator prefix::
2465
2466         sage: QuaternionHermitianEJA(2, prefix='a').gens()
2467         (a0, a1, a2, a3, a4, a5)
2468
2469     We can construct the (trivial) algebra of rank zero::
2470
2471         sage: QuaternionHermitianEJA(0)
2472         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2473
2474     """
2475     @classmethod
2476     def _denormalized_basis(cls, n, field):
2477         """
2478         Returns a basis for the space of quaternion Hermitian n-by-n matrices.
2479
2480         Why do we embed these? Basically, because all of numerical
2481         linear algebra assumes that you're working with vectors consisting
2482         of `n` entries from a field and scalars from the same field. There's
2483         no way to tell SageMath that (for example) the vectors contain
2484         complex numbers, while the scalar field is real.
2485
2486         SETUP::
2487
2488             sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2489
2490         TESTS::
2491
2492             sage: set_random_seed()
2493             sage: n = ZZ.random_element(1,5)
2494             sage: B = QuaternionHermitianEJA._denormalized_basis(n,ZZ)
2495             sage: all( M.is_symmetric() for M in B )
2496             True
2497
2498         """
2499         Q = QuaternionAlgebra(QQ,-1,-1)
2500         I,J,K = Q.gens()
2501
2502         # This is like the symmetric case, but we need to be careful:
2503         #
2504         #   * We want conjugate-symmetry, not just symmetry.
2505         #   * The diagonal will (as a result) be real.
2506         #
2507         S = []
2508         Eij = matrix.zero(Q,n)
2509         for i in range(n):
2510             for j in range(i+1):
2511                 # "build" E_ij
2512                 Eij[i,j] = 1
2513                 if i == j:
2514                     Sij = cls.real_embed(Eij)
2515                     S.append(Sij)
2516                 else:
2517                     # The second, third, and fourth ones have a minus
2518                     # because they're conjugated.
2519                     # Eij = Eij + Eij.transpose()
2520                     Eij[j,i] = 1
2521                     Sij_real = cls.real_embed(Eij)
2522                     S.append(Sij_real)
2523                     # Eij = I*(Eij - Eij.transpose())
2524                     Eij[i,j] = I
2525                     Eij[j,i] = -I
2526                     Sij_I = cls.real_embed(Eij)
2527                     S.append(Sij_I)
2528                     # Eij = J*(Eij - Eij.transpose())
2529                     Eij[i,j] = J
2530                     Eij[j,i] = -J
2531                     Sij_J = cls.real_embed(Eij)
2532                     S.append(Sij_J)
2533                     # Eij = K*(Eij - Eij.transpose())
2534                     Eij[i,j] = K
2535                     Eij[j,i] = -K
2536                     Sij_K = cls.real_embed(Eij)
2537                     S.append(Sij_K)
2538                     Eij[j,i] = 0
2539                 # "erase" E_ij
2540                 Eij[i,j] = 0
2541
2542         # Since we embedded the entries, we can drop back to the
2543         # desired real "field" instead of the quaternion algebra "Q".
2544         return tuple( s.change_ring(field) for s in S )
2545
2546
2547     def __init__(self, n, field=AA, **kwargs):
2548         # We know this is a valid EJA, but will double-check
2549         # if the user passes check_axioms=True.
2550         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2551
2552         associative = False
2553         if n <= 1:
2554             associative = True
2555
2556         super().__init__(self._denormalized_basis(n,field),
2557                          self.jordan_product,
2558                          self.trace_inner_product,
2559                          field=field,
2560                          associative=associative,
2561                          **kwargs)
2562
2563         # TODO: this could be factored out somehow, but is left here
2564         # because the MatrixEJA is not presently a subclass of the
2565         # FDEJA class that defines rank() and one().
2566         self.rank.set_cache(n)
2567         idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
2568         self.one.set_cache(self(idV))
2569
2570
2571     @staticmethod
2572     def _max_random_instance_size():
2573         r"""
2574         The maximum rank of a random QuaternionHermitianEJA.
2575         """
2576         return 2 # Dimension 6
2577
2578     @classmethod
2579     def random_instance(cls, **kwargs):
2580         """
2581         Return a random instance of this type of algebra.
2582         """
2583         n = ZZ.random_element(cls._max_random_instance_size() + 1)
2584         return cls(n, **kwargs)
2585
2586 class OctonionHermitianEJA(FiniteDimensionalEJA, MatrixEJA):
2587     r"""
2588     SETUP::
2589
2590         sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
2591         ....:                                  OctonionHermitianEJA)
2592
2593     EXAMPLES:
2594
2595     The 3-by-3 algebra satisfies the axioms of an EJA::
2596
2597         sage: OctonionHermitianEJA(3,                    # long time
2598         ....:                      field=QQ,             # long time
2599         ....:                      orthonormalize=False, # long time
2600         ....:                      check_axioms=True)    # long time
2601         Euclidean Jordan algebra of dimension 27 over Rational Field
2602
2603     After a change-of-basis, the 2-by-2 algebra has the same
2604     multiplication table as the ten-dimensional Jordan spin algebra::
2605
2606         sage: b = OctonionHermitianEJA._denormalized_basis(2,QQ)
2607         sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
2608         sage: jp = OctonionHermitianEJA.jordan_product
2609         sage: ip = OctonionHermitianEJA.trace_inner_product
2610         sage: J = FiniteDimensionalEJA(basis,
2611         ....:                          jp,
2612         ....:                          ip,
2613         ....:                          field=QQ,
2614         ....:                          orthonormalize=False)
2615         sage: J.multiplication_table()
2616         +----++----+----+----+----+----+----+----+----+----+----+
2617         | *  || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2618         +====++====+====+====+====+====+====+====+====+====+====+
2619         | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2620         +----++----+----+----+----+----+----+----+----+----+----+
2621         | b1 || b1 | b0 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | 0  |
2622         +----++----+----+----+----+----+----+----+----+----+----+
2623         | b2 || b2 | 0  | b0 | 0  | 0  | 0  | 0  | 0  | 0  | 0  |
2624         +----++----+----+----+----+----+----+----+----+----+----+
2625         | b3 || b3 | 0  | 0  | b0 | 0  | 0  | 0  | 0  | 0  | 0  |
2626         +----++----+----+----+----+----+----+----+----+----+----+
2627         | b4 || b4 | 0  | 0  | 0  | b0 | 0  | 0  | 0  | 0  | 0  |
2628         +----++----+----+----+----+----+----+----+----+----+----+
2629         | b5 || b5 | 0  | 0  | 0  | 0  | b0 | 0  | 0  | 0  | 0  |
2630         +----++----+----+----+----+----+----+----+----+----+----+
2631         | b6 || b6 | 0  | 0  | 0  | 0  | 0  | b0 | 0  | 0  | 0  |
2632         +----++----+----+----+----+----+----+----+----+----+----+
2633         | b7 || b7 | 0  | 0  | 0  | 0  | 0  | 0  | b0 | 0  | 0  |
2634         +----++----+----+----+----+----+----+----+----+----+----+
2635         | b8 || b8 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | b0 | 0  |
2636         +----++----+----+----+----+----+----+----+----+----+----+
2637         | b9 || b9 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | 0  | b0 |
2638         +----++----+----+----+----+----+----+----+----+----+----+
2639
2640     TESTS:
2641
2642     We can actually construct the 27-dimensional Albert algebra,
2643     and we get the right unit element if we recompute it::
2644
2645         sage: J = OctonionHermitianEJA(3,                    # long time
2646         ....:                          field=QQ,             # long time
2647         ....:                          orthonormalize=False) # long time
2648         sage: J.one.clear_cache()                            # long time
2649         sage: J.one()                                        # long time
2650         b0 + b9 + b26
2651         sage: J.one().to_matrix()                            # long time
2652         +----+----+----+
2653         | e0 | 0  | 0  |
2654         +----+----+----+
2655         | 0  | e0 | 0  |
2656         +----+----+----+
2657         | 0  | 0  | e0 |
2658         +----+----+----+
2659
2660     The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
2661     spin algebra, but just to be sure, we recompute its rank::
2662
2663         sage: J = OctonionHermitianEJA(2,                    # long time
2664         ....:                          field=QQ,             # long time
2665         ....:                          orthonormalize=False) # long time
2666         sage: J.rank.clear_cache()                           # long time
2667         sage: J.rank()                                       # long time
2668         2
2669     """
2670     def __init__(self, n, field=AA, **kwargs):
2671         if n > 3:
2672             # Otherwise we don't get an EJA.
2673             raise ValueError("n cannot exceed 3")
2674
2675         # We know this is a valid EJA, but will double-check
2676         # if the user passes check_axioms=True.
2677         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2678
2679         super().__init__(self._denormalized_basis(n,field),
2680                          self.jordan_product,
2681                          self.trace_inner_product,
2682                          field=field,
2683                          **kwargs)
2684
2685         # TODO: this could be factored out somehow, but is left here
2686         # because the MatrixEJA is not presently a subclass of the
2687         # FDEJA class that defines rank() and one().
2688         self.rank.set_cache(n)
2689         idV = self.matrix_space().one()
2690         self.one.set_cache(self(idV))
2691
2692
2693     @classmethod
2694     def _denormalized_basis(cls, n, field):
2695         """
2696         Returns a basis for the space of octonion Hermitian n-by-n
2697         matrices.
2698
2699         SETUP::
2700
2701             sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2702
2703         EXAMPLES::
2704
2705             sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ)
2706             sage: all( M.is_hermitian() for M in B )
2707             True
2708             sage: len(B)
2709             27
2710
2711         """
2712         from mjo.octonions import OctonionMatrixAlgebra
2713         MS = OctonionMatrixAlgebra(n, scalars=field)
2714         es = MS.entry_algebra().gens()
2715
2716         basis = []
2717         for i in range(n):
2718             for j in range(i+1):
2719                 if i == j:
2720                     E_ii = MS.monomial( (i,j,es[0]) )
2721                     basis.append(E_ii)
2722                 else:
2723                     for e in es:
2724                         E_ij  = MS.monomial( (i,j,e)             )
2725                         ec = e.conjugate()
2726                         # If the conjugate has a negative sign in front
2727                         # of it, (j,i,ec) won't be a monomial!
2728                         if (j,i,ec) in MS.indices():
2729                             E_ij += MS.monomial( (j,i,ec) )
2730                         else:
2731                             E_ij -= MS.monomial( (j,i,-ec) )
2732                         basis.append(E_ij)
2733
2734         return tuple( basis )
2735
2736     @staticmethod
2737     def trace_inner_product(X,Y):
2738         r"""
2739         The octonions don't know that the reals are embedded in them,
2740         so we have to take the e0 component ourselves.
2741
2742         SETUP::
2743
2744             sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
2745
2746         TESTS::
2747
2748             sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
2749             sage: I = J.one().to_matrix()
2750             sage: J.trace_inner_product(I, -I)
2751             -2
2752
2753         """
2754         return (X*Y).trace().real().coefficient(0)
2755
2756 class HadamardEJA(ConcreteEJA):
2757     """
2758     Return the Euclidean Jordan Algebra corresponding to the set
2759     `R^n` under the Hadamard product.
2760
2761     Note: this is nothing more than the Cartesian product of ``n``
2762     copies of the spin algebra. Once Cartesian product algebras
2763     are implemented, this can go.
2764
2765     SETUP::
2766
2767         sage: from mjo.eja.eja_algebra import HadamardEJA
2768
2769     EXAMPLES:
2770
2771     This multiplication table can be verified by hand::
2772
2773         sage: J = HadamardEJA(3)
2774         sage: b0,b1,b2 = J.gens()
2775         sage: b0*b0
2776         b0
2777         sage: b0*b1
2778         0
2779         sage: b0*b2
2780         0
2781         sage: b1*b1
2782         b1
2783         sage: b1*b2
2784         0
2785         sage: b2*b2
2786         b2
2787
2788     TESTS:
2789
2790     We can change the generator prefix::
2791
2792         sage: HadamardEJA(3, prefix='r').gens()
2793         (r0, r1, r2)
2794
2795     """
2796     def __init__(self, n, field=AA, **kwargs):
2797         if n == 0:
2798             jordan_product = lambda x,y: x
2799             inner_product = lambda x,y: x
2800         else:
2801             def jordan_product(x,y):
2802                 P = x.parent()
2803                 return P( xi*yi for (xi,yi) in zip(x,y) )
2804
2805             def inner_product(x,y):
2806                 return (x.T*y)[0,0]
2807
2808         # New defaults for keyword arguments. Don't orthonormalize
2809         # because our basis is already orthonormal with respect to our
2810         # inner-product. Don't check the axioms, because we know this
2811         # is a valid EJA... but do double-check if the user passes
2812         # check_axioms=True. Note: we DON'T override the "check_field"
2813         # default here, because the user can pass in a field!
2814         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
2815         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2816
2817         column_basis = tuple( b.column()
2818                               for b in FreeModule(field, n).basis() )
2819         super().__init__(column_basis,
2820                          jordan_product,
2821                          inner_product,
2822                          field=field,
2823                          associative=True,
2824                          **kwargs)
2825         self.rank.set_cache(n)
2826
2827         if n == 0:
2828             self.one.set_cache( self.zero() )
2829         else:
2830             self.one.set_cache( sum(self.gens()) )
2831
2832     @staticmethod
2833     def _max_random_instance_size():
2834         r"""
2835         The maximum dimension of a random HadamardEJA.
2836         """
2837         return 5
2838
2839     @classmethod
2840     def random_instance(cls, **kwargs):
2841         """
2842         Return a random instance of this type of algebra.
2843         """
2844         n = ZZ.random_element(cls._max_random_instance_size() + 1)
2845         return cls(n, **kwargs)
2846
2847
2848 class BilinearFormEJA(ConcreteEJA):
2849     r"""
2850     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2851     with the half-trace inner product and jordan product ``x*y =
2852     (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2853     a symmetric positive-definite "bilinear form" matrix. Its
2854     dimension is the size of `B`, and it has rank two in dimensions
2855     larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2856     the identity matrix of order ``n``.
2857
2858     We insist that the one-by-one upper-left identity block of `B` be
2859     passed in as well so that we can be passed a matrix of size zero
2860     to construct a trivial algebra.
2861
2862     SETUP::
2863
2864         sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2865         ....:                                  JordanSpinEJA)
2866
2867     EXAMPLES:
2868
2869     When no bilinear form is specified, the identity matrix is used,
2870     and the resulting algebra is the Jordan spin algebra::
2871
2872         sage: B = matrix.identity(AA,3)
2873         sage: J0 = BilinearFormEJA(B)
2874         sage: J1 = JordanSpinEJA(3)
2875         sage: J0.multiplication_table() == J0.multiplication_table()
2876         True
2877
2878     An error is raised if the matrix `B` does not correspond to a
2879     positive-definite bilinear form::
2880
2881         sage: B = matrix.random(QQ,2,3)
2882         sage: J = BilinearFormEJA(B)
2883         Traceback (most recent call last):
2884         ...
2885         ValueError: bilinear form is not positive-definite
2886         sage: B = matrix.zero(QQ,3)
2887         sage: J = BilinearFormEJA(B)
2888         Traceback (most recent call last):
2889         ...
2890         ValueError: bilinear form is not positive-definite
2891
2892     TESTS:
2893
2894     We can create a zero-dimensional algebra::
2895
2896         sage: B = matrix.identity(AA,0)
2897         sage: J = BilinearFormEJA(B)
2898         sage: J.basis()
2899         Finite family {}
2900
2901     We can check the multiplication condition given in the Jordan, von
2902     Neumann, and Wigner paper (and also discussed on my "On the
2903     symmetry..." paper). Note that this relies heavily on the standard
2904     choice of basis, as does anything utilizing the bilinear form
2905     matrix.  We opt not to orthonormalize the basis, because if we
2906     did, we would have to normalize the `s_{i}` in a similar manner::
2907
2908         sage: set_random_seed()
2909         sage: n = ZZ.random_element(5)
2910         sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2911         sage: B11 = matrix.identity(QQ,1)
2912         sage: B22 = M.transpose()*M
2913         sage: B = block_matrix(2,2,[ [B11,0  ],
2914         ....:                        [0, B22 ] ])
2915         sage: J = BilinearFormEJA(B, orthonormalize=False)
2916         sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2917         sage: V = J.vector_space()
2918         sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2919         ....:         for ei in eis ]
2920         sage: actual = [ sis[i]*sis[j]
2921         ....:            for i in range(n-1)
2922         ....:            for j in range(n-1) ]
2923         sage: expected = [ J.one() if i == j else J.zero()
2924         ....:              for i in range(n-1)
2925         ....:              for j in range(n-1) ]
2926         sage: actual == expected
2927         True
2928
2929     """
2930     def __init__(self, B, field=AA, **kwargs):
2931         # The matrix "B" is supplied by the user in most cases,
2932         # so it makes sense to check whether or not its positive-
2933         # definite unless we are specifically asked not to...
2934         if ("check_axioms" not in kwargs) or kwargs["check_axioms"]:
2935             if not B.is_positive_definite():
2936                 raise ValueError("bilinear form is not positive-definite")
2937
2938         # However, all of the other data for this EJA is computed
2939         # by us in manner that guarantees the axioms are
2940         # satisfied. So, again, unless we are specifically asked to
2941         # verify things, we'll skip the rest of the checks.
2942         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2943
2944         def inner_product(x,y):
2945             return (y.T*B*x)[0,0]
2946
2947         def jordan_product(x,y):
2948             P = x.parent()
2949             x0 = x[0,0]
2950             xbar = x[1:,0]
2951             y0 = y[0,0]
2952             ybar = y[1:,0]
2953             z0 = inner_product(y,x)
2954             zbar = y0*xbar + x0*ybar
2955             return P([z0] + zbar.list())
2956
2957         n = B.nrows()
2958         column_basis = tuple( b.column()
2959                               for b in FreeModule(field, n).basis() )
2960
2961         # TODO: I haven't actually checked this, but it seems legit.
2962         associative = False
2963         if n <= 2:
2964             associative = True
2965
2966         super().__init__(column_basis,
2967                          jordan_product,
2968                          inner_product,
2969                          field=field,
2970                          associative=associative,
2971                          **kwargs)
2972
2973         # The rank of this algebra is two, unless we're in a
2974         # one-dimensional ambient space (because the rank is bounded
2975         # by the ambient dimension).
2976         self.rank.set_cache(min(n,2))
2977
2978         if n == 0:
2979             self.one.set_cache( self.zero() )
2980         else:
2981             self.one.set_cache( self.monomial(0) )
2982
2983     @staticmethod
2984     def _max_random_instance_size():
2985         r"""
2986         The maximum dimension of a random BilinearFormEJA.
2987         """
2988         return 5
2989
2990     @classmethod
2991     def random_instance(cls, **kwargs):
2992         """
2993         Return a random instance of this algebra.
2994         """
2995         n = ZZ.random_element(cls._max_random_instance_size() + 1)
2996         if n.is_zero():
2997             B = matrix.identity(ZZ, n)
2998             return cls(B, **kwargs)
2999
3000         B11 = matrix.identity(ZZ, 1)
3001         M = matrix.random(ZZ, n-1)
3002         I = matrix.identity(ZZ, n-1)
3003         alpha = ZZ.zero()
3004         while alpha.is_zero():
3005             alpha = ZZ.random_element().abs()
3006         B22 = M.transpose()*M + alpha*I
3007
3008         from sage.matrix.special import block_matrix
3009         B = block_matrix(2,2, [ [B11,   ZZ(0) ],
3010                                 [ZZ(0), B22 ] ])
3011
3012         return cls(B, **kwargs)
3013
3014
3015 class JordanSpinEJA(BilinearFormEJA):
3016     """
3017     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
3018     with the usual inner product and jordan product ``x*y =
3019     (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
3020     the reals.
3021
3022     SETUP::
3023
3024         sage: from mjo.eja.eja_algebra import JordanSpinEJA
3025
3026     EXAMPLES:
3027
3028     This multiplication table can be verified by hand::
3029
3030         sage: J = JordanSpinEJA(4)
3031         sage: b0,b1,b2,b3 = J.gens()
3032         sage: b0*b0
3033         b0
3034         sage: b0*b1
3035         b1
3036         sage: b0*b2
3037         b2
3038         sage: b0*b3
3039         b3
3040         sage: b1*b2
3041         0
3042         sage: b1*b3
3043         0
3044         sage: b2*b3
3045         0
3046
3047     We can change the generator prefix::
3048
3049         sage: JordanSpinEJA(2, prefix='B').gens()
3050         (B0, B1)
3051
3052     TESTS:
3053
3054         Ensure that we have the usual inner product on `R^n`::
3055
3056             sage: set_random_seed()
3057             sage: J = JordanSpinEJA.random_instance()
3058             sage: x,y = J.random_elements(2)
3059             sage: actual = x.inner_product(y)
3060             sage: expected = x.to_vector().inner_product(y.to_vector())
3061             sage: actual == expected
3062             True
3063
3064     """
3065     def __init__(self, n, *args, **kwargs):
3066         # This is a special case of the BilinearFormEJA with the
3067         # identity matrix as its bilinear form.
3068         B = matrix.identity(ZZ, n)
3069
3070         # Don't orthonormalize because our basis is already
3071         # orthonormal with respect to our inner-product.
3072         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
3073
3074         # But also don't pass check_field=False here, because the user
3075         # can pass in a field!
3076         super().__init__(B, *args, **kwargs)
3077
3078     @staticmethod
3079     def _max_random_instance_size():
3080         r"""
3081         The maximum dimension of a random JordanSpinEJA.
3082         """
3083         return 5
3084
3085     @classmethod
3086     def random_instance(cls, **kwargs):
3087         """
3088         Return a random instance of this type of algebra.
3089
3090         Needed here to override the implementation for ``BilinearFormEJA``.
3091         """
3092         n = ZZ.random_element(cls._max_random_instance_size() + 1)
3093         return cls(n, **kwargs)
3094
3095
3096 class TrivialEJA(ConcreteEJA):
3097     """
3098     The trivial Euclidean Jordan algebra consisting of only a zero element.
3099
3100     SETUP::
3101
3102         sage: from mjo.eja.eja_algebra import TrivialEJA
3103
3104     EXAMPLES::
3105
3106         sage: J = TrivialEJA()
3107         sage: J.dimension()
3108         0
3109         sage: J.zero()
3110         0
3111         sage: J.one()
3112         0
3113         sage: 7*J.one()*12*J.one()
3114         0
3115         sage: J.one().inner_product(J.one())
3116         0
3117         sage: J.one().norm()
3118         0
3119         sage: J.one().subalgebra_generated_by()
3120         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
3121         sage: J.rank()
3122         0
3123
3124     """
3125     def __init__(self, **kwargs):
3126         jordan_product = lambda x,y: x
3127         inner_product = lambda x,y: 0
3128         basis = ()
3129
3130         # New defaults for keyword arguments
3131         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
3132         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
3133
3134         super().__init__(basis,
3135                          jordan_product,
3136                          inner_product,
3137                          associative=True,
3138                          **kwargs)
3139
3140         # The rank is zero using my definition, namely the dimension of the
3141         # largest subalgebra generated by any element.
3142         self.rank.set_cache(0)
3143         self.one.set_cache( self.zero() )
3144
3145     @classmethod
3146     def random_instance(cls, **kwargs):
3147         # We don't take a "size" argument so the superclass method is
3148         # inappropriate for us.
3149         return cls(**kwargs)
3150
3151
3152 class CartesianProductEJA(FiniteDimensionalEJA):
3153     r"""
3154     The external (orthogonal) direct sum of two or more Euclidean
3155     Jordan algebras. Every Euclidean Jordan algebra decomposes into an
3156     orthogonal direct sum of simple Euclidean Jordan algebras which is
3157     then isometric to a Cartesian product, so no generality is lost by
3158     providing only this construction.
3159
3160     SETUP::
3161
3162         sage: from mjo.eja.eja_algebra import (random_eja,
3163         ....:                                  CartesianProductEJA,
3164         ....:                                  HadamardEJA,
3165         ....:                                  JordanSpinEJA,
3166         ....:                                  RealSymmetricEJA)
3167
3168     EXAMPLES:
3169
3170     The Jordan product is inherited from our factors and implemented by
3171     our CombinatorialFreeModule Cartesian product superclass::
3172
3173         sage: set_random_seed()
3174         sage: J1 = HadamardEJA(2)
3175         sage: J2 = RealSymmetricEJA(2)
3176         sage: J = cartesian_product([J1,J2])
3177         sage: x,y = J.random_elements(2)
3178         sage: x*y in J
3179         True
3180
3181     The ability to retrieve the original factors is implemented by our
3182     CombinatorialFreeModule Cartesian product superclass::
3183
3184         sage: J1 = HadamardEJA(2, field=QQ)
3185         sage: J2 = JordanSpinEJA(3, field=QQ)
3186         sage: J = cartesian_product([J1,J2])
3187         sage: J.cartesian_factors()
3188         (Euclidean Jordan algebra of dimension 2 over Rational Field,
3189          Euclidean Jordan algebra of dimension 3 over Rational Field)
3190
3191     You can provide more than two factors::
3192
3193         sage: J1 = HadamardEJA(2)
3194         sage: J2 = JordanSpinEJA(3)
3195         sage: J3 = RealSymmetricEJA(3)
3196         sage: cartesian_product([J1,J2,J3])
3197         Euclidean Jordan algebra of dimension 2 over Algebraic Real
3198         Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
3199         Real Field (+) Euclidean Jordan algebra of dimension 6 over
3200         Algebraic Real Field
3201
3202     Rank is additive on a Cartesian product::
3203
3204         sage: J1 = HadamardEJA(1)
3205         sage: J2 = RealSymmetricEJA(2)
3206         sage: J = cartesian_product([J1,J2])
3207         sage: J1.rank.clear_cache()
3208         sage: J2.rank.clear_cache()
3209         sage: J.rank.clear_cache()
3210         sage: J.rank()
3211         3
3212         sage: J.rank() == J1.rank() + J2.rank()
3213         True
3214
3215     The same rank computation works over the rationals, with whatever
3216     basis you like::
3217
3218         sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
3219         sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
3220         sage: J = cartesian_product([J1,J2])
3221         sage: J1.rank.clear_cache()
3222         sage: J2.rank.clear_cache()
3223         sage: J.rank.clear_cache()
3224         sage: J.rank()
3225         3
3226         sage: J.rank() == J1.rank() + J2.rank()
3227         True
3228
3229     The product algebra will be associative if and only if all of its
3230     components are associative::
3231
3232         sage: J1 = HadamardEJA(2)
3233         sage: J1.is_associative()
3234         True
3235         sage: J2 = HadamardEJA(3)
3236         sage: J2.is_associative()
3237         True
3238         sage: J3 = RealSymmetricEJA(3)
3239         sage: J3.is_associative()
3240         False
3241         sage: CP1 = cartesian_product([J1,J2])
3242         sage: CP1.is_associative()
3243         True
3244         sage: CP2 = cartesian_product([J1,J3])
3245         sage: CP2.is_associative()
3246         False
3247
3248     Cartesian products of Cartesian products work::
3249
3250         sage: J1 = JordanSpinEJA(1)
3251         sage: J2 = JordanSpinEJA(1)
3252         sage: J3 = JordanSpinEJA(1)
3253         sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
3254         sage: J.multiplication_table()
3255         +----++----+----+----+
3256         | *  || b0 | b1 | b2 |
3257         +====++====+====+====+
3258         | b0 || b0 | 0  | 0  |
3259         +----++----+----+----+
3260         | b1 || 0  | b1 | 0  |
3261         +----++----+----+----+
3262         | b2 || 0  | 0  | b2 |
3263         +----++----+----+----+
3264         sage: HadamardEJA(3).multiplication_table()
3265         +----++----+----+----+
3266         | *  || b0 | b1 | b2 |
3267         +====++====+====+====+
3268         | b0 || b0 | 0  | 0  |
3269         +----++----+----+----+
3270         | b1 || 0  | b1 | 0  |
3271         +----++----+----+----+
3272         | b2 || 0  | 0  | b2 |
3273         +----++----+----+----+
3274
3275     TESTS:
3276
3277     All factors must share the same base field::
3278
3279         sage: J1 = HadamardEJA(2, field=QQ)
3280         sage: J2 = RealSymmetricEJA(2)
3281         sage: CartesianProductEJA((J1,J2))
3282         Traceback (most recent call last):
3283         ...
3284         ValueError: all factors must share the same base field
3285
3286     The cached unit element is the same one that would be computed::
3287
3288         sage: set_random_seed()              # long time
3289         sage: J1 = random_eja()              # long time
3290         sage: J2 = random_eja()              # long time
3291         sage: J = cartesian_product([J1,J2]) # long time
3292         sage: actual = J.one()               # long time
3293         sage: J.one.clear_cache()            # long time
3294         sage: expected = J.one()             # long time
3295         sage: actual == expected             # long time
3296         True
3297
3298     """
3299     Element = FiniteDimensionalEJAElement
3300
3301
3302     def __init__(self, factors, **kwargs):
3303         m = len(factors)
3304         if m == 0:
3305             return TrivialEJA()
3306
3307         self._sets = factors
3308
3309         field = factors[0].base_ring()
3310         if not all( J.base_ring() == field for J in factors ):
3311             raise ValueError("all factors must share the same base field")
3312
3313         associative = all( f.is_associative() for f in factors )
3314
3315         MS = self.matrix_space()
3316         basis = []
3317         zero = MS.zero()
3318         for i in range(m):
3319             for b in factors[i].matrix_basis():
3320                 z = list(zero)
3321                 z[i] = b
3322                 basis.append(z)
3323
3324         basis = tuple( MS(b) for b in basis )
3325
3326         # Define jordan/inner products that operate on that matrix_basis.
3327         def jordan_product(x,y):
3328             return MS(tuple(
3329                 (factors[i](x[i])*factors[i](y[i])).to_matrix()
3330                 for i in range(m)
3331             ))
3332
3333         def inner_product(x, y):
3334             return sum(
3335                 factors[i](x[i]).inner_product(factors[i](y[i]))
3336                 for i in range(m)
3337             )
3338
3339         # There's no need to check the field since it already came
3340         # from an EJA. Likewise the axioms are guaranteed to be
3341         # satisfied, unless the guy writing this class sucks.
3342         #
3343         # If you want the basis to be orthonormalized, orthonormalize
3344         # the factors.
3345         FiniteDimensionalEJA.__init__(self,
3346                                       basis,
3347                                       jordan_product,
3348                                       inner_product,
3349                                       field=field,
3350                                       orthonormalize=False,
3351                                       associative=associative,
3352                                       cartesian_product=True,
3353                                       check_field=False,
3354                                       check_axioms=False)
3355
3356         ones = tuple(J.one().to_matrix() for J in factors)
3357         self.one.set_cache(self(ones))
3358         self.rank.set_cache(sum(J.rank() for J in factors))
3359
3360     def cartesian_factors(self):
3361         # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3362         return self._sets
3363
3364     def cartesian_factor(self, i):
3365         r"""
3366         Return the ``i``th factor of this algebra.
3367         """
3368         return self._sets[i]
3369
3370     def _repr_(self):
3371         # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
3372         from sage.categories.cartesian_product import cartesian_product
3373         return cartesian_product.symbol.join("%s" % factor
3374                                              for factor in self._sets)
3375
3376     def matrix_space(self):
3377         r"""
3378         Return the space that our matrix basis lives in as a Cartesian
3379         product.
3380
3381         SETUP::
3382
3383             sage: from mjo.eja.eja_algebra import (HadamardEJA,
3384             ....:                                  RealSymmetricEJA)
3385
3386         EXAMPLES::
3387
3388             sage: J1 = HadamardEJA(1)
3389             sage: J2 = RealSymmetricEJA(2)
3390             sage: J = cartesian_product([J1,J2])
3391             sage: J.matrix_space()
3392             The Cartesian product of (Full MatrixSpace of 1 by 1 dense
3393             matrices over Algebraic Real Field, Full MatrixSpace of 2
3394             by 2 dense matrices over Algebraic Real Field)
3395
3396         """
3397         from sage.categories.cartesian_product import cartesian_product
3398         return cartesian_product( [J.matrix_space()
3399                                    for J in self.cartesian_factors()] )
3400
3401     @cached_method
3402     def cartesian_projection(self, i):
3403         r"""
3404         SETUP::
3405
3406             sage: from mjo.eja.eja_algebra import (random_eja,
3407             ....:                                  JordanSpinEJA,
3408             ....:                                  HadamardEJA,
3409             ....:                                  RealSymmetricEJA,
3410             ....:                                  ComplexHermitianEJA)
3411
3412         EXAMPLES:
3413
3414         The projection morphisms are Euclidean Jordan algebra
3415         operators::
3416
3417             sage: J1 = HadamardEJA(2)
3418             sage: J2 = RealSymmetricEJA(2)
3419             sage: J = cartesian_product([J1,J2])
3420             sage: J.cartesian_projection(0)
3421             Linear operator between finite-dimensional Euclidean Jordan
3422             algebras represented by the matrix:
3423             [1 0 0 0 0]
3424             [0 1 0 0 0]
3425             Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3426             Real Field (+) Euclidean Jordan algebra of dimension 3 over
3427             Algebraic Real Field
3428             Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3429             Real Field
3430             sage: J.cartesian_projection(1)
3431             Linear operator between finite-dimensional Euclidean Jordan
3432             algebras represented by the matrix:
3433             [0 0 1 0 0]
3434             [0 0 0 1 0]
3435             [0 0 0 0 1]
3436             Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3437             Real Field (+) Euclidean Jordan algebra of dimension 3 over
3438             Algebraic Real Field
3439             Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3440             Real Field
3441
3442         The projections work the way you'd expect on the vector
3443         representation of an element::
3444
3445             sage: J1 = JordanSpinEJA(2)
3446             sage: J2 = ComplexHermitianEJA(2)
3447             sage: J = cartesian_product([J1,J2])
3448             sage: pi_left = J.cartesian_projection(0)
3449             sage: pi_right = J.cartesian_projection(1)
3450             sage: pi_left(J.one()).to_vector()
3451             (1, 0)
3452             sage: pi_right(J.one()).to_vector()
3453             (1, 0, 0, 1)
3454             sage: J.one().to_vector()
3455             (1, 0, 1, 0, 0, 1)
3456
3457         TESTS:
3458
3459         The answer never changes::
3460
3461             sage: set_random_seed()
3462             sage: J1 = random_eja()
3463             sage: J2 = random_eja()
3464             sage: J = cartesian_product([J1,J2])
3465             sage: P0 = J.cartesian_projection(0)
3466             sage: P1 = J.cartesian_projection(0)
3467             sage: P0 == P1
3468             True
3469
3470         """
3471         offset = sum( self.cartesian_factor(k).dimension()
3472                       for k in range(i) )
3473         Ji = self.cartesian_factor(i)
3474         Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
3475                                    codomain=Ji)
3476
3477         return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
3478
3479     @cached_method
3480     def cartesian_embedding(self, i):
3481         r"""
3482         SETUP::
3483
3484             sage: from mjo.eja.eja_algebra import (random_eja,
3485             ....:                                  JordanSpinEJA,
3486             ....:                                  HadamardEJA,
3487             ....:                                  RealSymmetricEJA)
3488
3489         EXAMPLES:
3490
3491         The embedding morphisms are Euclidean Jordan algebra
3492         operators::
3493
3494             sage: J1 = HadamardEJA(2)
3495             sage: J2 = RealSymmetricEJA(2)
3496             sage: J = cartesian_product([J1,J2])
3497             sage: J.cartesian_embedding(0)
3498             Linear operator between finite-dimensional Euclidean Jordan
3499             algebras represented by the matrix:
3500             [1 0]
3501             [0 1]
3502             [0 0]
3503             [0 0]
3504             [0 0]
3505             Domain: Euclidean Jordan algebra of dimension 2 over
3506             Algebraic Real Field
3507             Codomain: Euclidean Jordan algebra of dimension 2 over
3508             Algebraic Real Field (+) Euclidean Jordan algebra of
3509             dimension 3 over Algebraic Real Field
3510             sage: J.cartesian_embedding(1)
3511             Linear operator between finite-dimensional Euclidean Jordan
3512             algebras represented by the matrix:
3513             [0 0 0]
3514             [0 0 0]
3515             [1 0 0]
3516             [0 1 0]
3517             [0 0 1]
3518             Domain: Euclidean Jordan algebra of dimension 3 over
3519             Algebraic Real Field
3520             Codomain: Euclidean Jordan algebra of dimension 2 over
3521             Algebraic Real Field (+) Euclidean Jordan algebra of
3522             dimension 3 over Algebraic Real Field
3523
3524         The embeddings work the way you'd expect on the vector
3525         representation of an element::
3526
3527             sage: J1 = JordanSpinEJA(3)
3528             sage: J2 = RealSymmetricEJA(2)
3529             sage: J = cartesian_product([J1,J2])
3530             sage: iota_left = J.cartesian_embedding(0)
3531             sage: iota_right = J.cartesian_embedding(1)
3532             sage: iota_left(J1.zero()) == J.zero()
3533             True
3534             sage: iota_right(J2.zero()) == J.zero()
3535             True
3536             sage: J1.one().to_vector()
3537             (1, 0, 0)
3538             sage: iota_left(J1.one()).to_vector()
3539             (1, 0, 0, 0, 0, 0)
3540             sage: J2.one().to_vector()
3541             (1, 0, 1)
3542             sage: iota_right(J2.one()).to_vector()
3543             (0, 0, 0, 1, 0, 1)
3544             sage: J.one().to_vector()
3545             (1, 0, 0, 1, 0, 1)
3546
3547         TESTS:
3548
3549         The answer never changes::
3550
3551             sage: set_random_seed()
3552             sage: J1 = random_eja()
3553             sage: J2 = random_eja()
3554             sage: J = cartesian_product([J1,J2])
3555             sage: E0 = J.cartesian_embedding(0)
3556             sage: E1 = J.cartesian_embedding(0)
3557             sage: E0 == E1
3558             True
3559
3560         Composing a projection with the corresponding inclusion should
3561         produce the identity map, and mismatching them should produce
3562         the zero map::
3563
3564             sage: set_random_seed()
3565             sage: J1 = random_eja()
3566             sage: J2 = random_eja()
3567             sage: J = cartesian_product([J1,J2])
3568             sage: iota_left = J.cartesian_embedding(0)
3569             sage: iota_right = J.cartesian_embedding(1)
3570             sage: pi_left = J.cartesian_projection(0)
3571             sage: pi_right = J.cartesian_projection(1)
3572             sage: pi_left*iota_left == J1.one().operator()
3573             True
3574             sage: pi_right*iota_right == J2.one().operator()
3575             True
3576             sage: (pi_left*iota_right).is_zero()
3577             True
3578             sage: (pi_right*iota_left).is_zero()
3579             True
3580
3581         """
3582         offset = sum( self.cartesian_factor(k).dimension()
3583                       for k in range(i) )
3584         Ji = self.cartesian_factor(i)
3585         Ei = Ji._module_morphism(lambda j: self.monomial(j + offset),
3586                                  codomain=self)
3587         return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
3588
3589
3590
3591 FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
3592
3593 class RationalBasisCartesianProductEJA(CartesianProductEJA,
3594                                        RationalBasisEJA):
3595     r"""
3596     A separate class for products of algebras for which we know a
3597     rational basis.
3598
3599     SETUP::
3600
3601         sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
3602         ....:                                  RealSymmetricEJA)
3603
3604     EXAMPLES:
3605
3606     This gives us fast characteristic polynomial computations in
3607     product algebras, too::
3608
3609
3610         sage: J1 = JordanSpinEJA(2)
3611         sage: J2 = RealSymmetricEJA(3)
3612         sage: J = cartesian_product([J1,J2])
3613         sage: J.characteristic_polynomial_of().degree()
3614         5
3615         sage: J.rank()
3616         5
3617
3618     """
3619     def __init__(self, algebras, **kwargs):
3620         CartesianProductEJA.__init__(self, algebras, **kwargs)
3621
3622         self._rational_algebra = None
3623         if self.vector_space().base_field() is not QQ:
3624             self._rational_algebra = cartesian_product([
3625                 r._rational_algebra for r in algebras
3626             ])
3627
3628
3629 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
3630
3631 def random_eja(*args, **kwargs):
3632     J1 = ConcreteEJA.random_instance(*args, **kwargs)
3633
3634     # This might make Cartesian products appear roughly as often as
3635     # any other ConcreteEJA.
3636     if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
3637         # Use random_eja() again so we can get more than two factors.
3638         J2 = random_eja(*args, **kwargs)
3639         J = cartesian_product([J1,J2])
3640         return J
3641     else:
3642         return J1