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