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