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