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