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_size():
1719         """
1720         Return an integer "size" that is an upper bound on the size of
1721         this algebra when it is used in a random test
1722         case. Unfortunately, the term "size" is ambiguous -- when
1723         dealing with `R^n` under either the Hadamard or Jordan spin
1724         product, the "size" refers to the dimension `n`. When dealing
1725         with a matrix algebra (real symmetric or complex/quaternion
1726         Hermitian), it refers to the size of the matrix, which is far
1727         less than the dimension of the underlying vector space.
1728
1729         This method must be implemented in each subclass.
1730         """
1731         raise NotImplementedError
1732
1733     @classmethod
1734     def random_instance(cls, *args, **kwargs):
1735         """
1736         Return a random instance of this type of algebra.
1737
1738         This method should be implemented in each subclass.
1739         """
1740         from sage.misc.prandom import choice
1741         eja_class = choice(cls.__subclasses__())
1742
1743         # These all bubble up to the RationalBasisEJA superclass
1744         # constructor, so any (kw)args valid there are also valid
1745         # here.
1746         return eja_class.random_instance(*args, **kwargs)
1747
1748
1749 class MatrixEJA:
1750     @staticmethod
1751     def _denormalized_basis(A):
1752         """
1753         Returns a basis for the space of complex Hermitian n-by-n matrices.
1754
1755         Why do we embed these? Basically, because all of numerical linear
1756         algebra assumes that you're working with vectors consisting of `n`
1757         entries from a field and scalars from the same field. There's no way
1758         to tell SageMath that (for example) the vectors contain complex
1759         numbers, while the scalar field is real.
1760
1761         SETUP::
1762
1763             sage: from mjo.hurwitz import (ComplexMatrixAlgebra,
1764             ....:                          QuaternionMatrixAlgebra,
1765             ....:                          OctonionMatrixAlgebra)
1766             sage: from mjo.eja.eja_algebra import MatrixEJA
1767
1768         TESTS::
1769
1770             sage: set_random_seed()
1771             sage: n = ZZ.random_element(1,5)
1772             sage: A = MatrixSpace(QQ, n)
1773             sage: B = MatrixEJA._denormalized_basis(A)
1774             sage: all( M.is_hermitian() for M in  B)
1775             True
1776
1777         ::
1778
1779             sage: set_random_seed()
1780             sage: n = ZZ.random_element(1,5)
1781             sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
1782             sage: B = MatrixEJA._denormalized_basis(A)
1783             sage: all( M.is_hermitian() for M in  B)
1784             True
1785
1786         ::
1787
1788             sage: set_random_seed()
1789             sage: n = ZZ.random_element(1,5)
1790             sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
1791             sage: B = MatrixEJA._denormalized_basis(A)
1792             sage: all( M.is_hermitian() for M in B )
1793             True
1794
1795         ::
1796
1797             sage: set_random_seed()
1798             sage: n = ZZ.random_element(1,5)
1799             sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
1800             sage: B = MatrixEJA._denormalized_basis(A)
1801             sage: all( M.is_hermitian() for M in B )
1802             True
1803
1804         """
1805         # These work for real MatrixSpace, whose monomials only have
1806         # two coordinates (because the last one would always be "1").
1807         es = A.base_ring().gens()
1808         gen = lambda A,m: A.monomial(m[:2])
1809
1810         if hasattr(A, 'entry_algebra_gens'):
1811             # We've got a MatrixAlgebra, and its monomials will have
1812             # three coordinates.
1813             es = A.entry_algebra_gens()
1814             gen = lambda A,m: A.monomial(m)
1815
1816         basis = []
1817         for i in range(A.nrows()):
1818             for j in range(i+1):
1819                 if i == j:
1820                     E_ii = gen(A, (i,j,es[0]))
1821                     basis.append(E_ii)
1822                 else:
1823                     for e in es:
1824                         E_ij  = gen(A, (i,j,e))
1825                         E_ij += E_ij.conjugate_transpose()
1826                         basis.append(E_ij)
1827
1828         return tuple( basis )
1829
1830     @staticmethod
1831     def jordan_product(X,Y):
1832         return (X*Y + Y*X)/2
1833
1834     @staticmethod
1835     def trace_inner_product(X,Y):
1836         r"""
1837         A trace inner-product for matrices that aren't embedded in the
1838         reals. It takes MATRICES as arguments, not EJA elements.
1839
1840         SETUP::
1841
1842             sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
1843             ....:                                  ComplexHermitianEJA,
1844             ....:                                  QuaternionHermitianEJA,
1845             ....:                                  OctonionHermitianEJA)
1846
1847         EXAMPLES::
1848
1849             sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False)
1850             sage: I = J.one().to_matrix()
1851             sage: J.trace_inner_product(I, -I)
1852             -2
1853
1854         ::
1855
1856             sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
1857             sage: I = J.one().to_matrix()
1858             sage: J.trace_inner_product(I, -I)
1859             -2
1860
1861         ::
1862
1863             sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
1864             sage: I = J.one().to_matrix()
1865             sage: J.trace_inner_product(I, -I)
1866             -2
1867
1868         ::
1869
1870             sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
1871             sage: I = J.one().to_matrix()
1872             sage: J.trace_inner_product(I, -I)
1873             -2
1874
1875         """
1876         tr = (X*Y).trace()
1877         if hasattr(tr, 'coefficient'):
1878             # Works for octonions, and has to come first because they
1879             # also have a "real()" method that doesn't return an
1880             # element of the scalar ring.
1881             return tr.coefficient(0)
1882         elif hasattr(tr, 'coefficient_tuple'):
1883             # Works for quaternions.
1884             return tr.coefficient_tuple()[0]
1885
1886         # Works for real and complex numbers.
1887         return tr.real()
1888
1889
1890
1891 class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
1892     """
1893     The rank-n simple EJA consisting of real symmetric n-by-n
1894     matrices, the usual symmetric Jordan product, and the trace inner
1895     product. It has dimension `(n^2 + n)/2` over the reals.
1896
1897     SETUP::
1898
1899         sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1900
1901     EXAMPLES::
1902
1903         sage: J = RealSymmetricEJA(2)
1904         sage: b0, b1, b2 = J.gens()
1905         sage: b0*b0
1906         b0
1907         sage: b1*b1
1908         1/2*b0 + 1/2*b2
1909         sage: b2*b2
1910         b2
1911
1912     In theory, our "field" can be any subfield of the reals::
1913
1914         sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
1915         Euclidean Jordan algebra of dimension 3 over Real Double Field
1916         sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
1917         Euclidean Jordan algebra of dimension 3 over Real Field with
1918         53 bits of precision
1919
1920     TESTS:
1921
1922     The dimension of this algebra is `(n^2 + n) / 2`::
1923
1924         sage: set_random_seed()
1925         sage: n_max = RealSymmetricEJA._max_random_instance_size()
1926         sage: n = ZZ.random_element(1, n_max)
1927         sage: J = RealSymmetricEJA(n)
1928         sage: J.dimension() == (n^2 + n)/2
1929         True
1930
1931     The Jordan multiplication is what we think it is::
1932
1933         sage: set_random_seed()
1934         sage: J = RealSymmetricEJA.random_instance()
1935         sage: x,y = J.random_elements(2)
1936         sage: actual = (x*y).to_matrix()
1937         sage: X = x.to_matrix()
1938         sage: Y = y.to_matrix()
1939         sage: expected = (X*Y + Y*X)/2
1940         sage: actual == expected
1941         True
1942         sage: J(expected) == x*y
1943         True
1944
1945     We can change the generator prefix::
1946
1947         sage: RealSymmetricEJA(3, prefix='q').gens()
1948         (q0, q1, q2, q3, q4, q5)
1949
1950     We can construct the (trivial) algebra of rank zero::
1951
1952         sage: RealSymmetricEJA(0)
1953         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
1954
1955     """
1956     @staticmethod
1957     def _max_random_instance_size():
1958         return 4 # Dimension 10
1959
1960     @classmethod
1961     def random_instance(cls, **kwargs):
1962         """
1963         Return a random instance of this type of algebra.
1964         """
1965         n = ZZ.random_element(cls._max_random_instance_size() + 1)
1966         return cls(n, **kwargs)
1967
1968     def __init__(self, n, field=AA, **kwargs):
1969         # We know this is a valid EJA, but will double-check
1970         # if the user passes check_axioms=True.
1971         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
1972
1973         A = MatrixSpace(field, n)
1974         super().__init__(self._denormalized_basis(A),
1975                          self.jordan_product,
1976                          self.trace_inner_product,
1977                          field=field,
1978                          matrix_space=A,
1979                          **kwargs)
1980
1981         # TODO: this could be factored out somehow, but is left here
1982         # because the MatrixEJA is not presently a subclass of the
1983         # FDEJA class that defines rank() and one().
1984         self.rank.set_cache(n)
1985         self.one.set_cache( self(A.one()) )
1986
1987
1988
1989 class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
1990     """
1991     The rank-n simple EJA consisting of complex Hermitian n-by-n
1992     matrices over the real numbers, the usual symmetric Jordan product,
1993     and the real-part-of-trace inner product. It has dimension `n^2` over
1994     the reals.
1995
1996     SETUP::
1997
1998         sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1999
2000     EXAMPLES:
2001
2002     In theory, our "field" can be any subfield of the reals, but we
2003     can't use inexact real fields at the moment because SageMath
2004     doesn't know how to convert their elements into complex numbers,
2005     or even into algebraic reals::
2006
2007         sage: QQbar(RDF(1))
2008         Traceback (most recent call last):
2009         ...
2010         TypeError: Illegal initializer for algebraic number
2011         sage: AA(RR(1))
2012         Traceback (most recent call last):
2013         ...
2014         TypeError: Illegal initializer for algebraic number
2015
2016     This causes the following error when we try to scale a matrix of
2017     complex numbers by an inexact real number::
2018
2019         sage: ComplexHermitianEJA(2,field=RR)
2020         Traceback (most recent call last):
2021         ...
2022         TypeError: Unable to coerce entries (=(1.00000000000000,
2023         -0.000000000000000)) to coefficients in Algebraic Real Field
2024
2025     TESTS:
2026
2027     The dimension of this algebra is `n^2`::
2028
2029         sage: set_random_seed()
2030         sage: n_max = ComplexHermitianEJA._max_random_instance_size()
2031         sage: n = ZZ.random_element(1, n_max)
2032         sage: J = ComplexHermitianEJA(n)
2033         sage: J.dimension() == n^2
2034         True
2035
2036     The Jordan multiplication is what we think it is::
2037
2038         sage: set_random_seed()
2039         sage: J = ComplexHermitianEJA.random_instance()
2040         sage: x,y = J.random_elements(2)
2041         sage: actual = (x*y).to_matrix()
2042         sage: X = x.to_matrix()
2043         sage: Y = y.to_matrix()
2044         sage: expected = (X*Y + Y*X)/2
2045         sage: actual == expected
2046         True
2047         sage: J(expected) == x*y
2048         True
2049
2050     We can change the generator prefix::
2051
2052         sage: ComplexHermitianEJA(2, prefix='z').gens()
2053         (z0, z1, z2, z3)
2054
2055     We can construct the (trivial) algebra of rank zero::
2056
2057         sage: ComplexHermitianEJA(0)
2058         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2059     """
2060     def __init__(self, n, field=AA, **kwargs):
2061         # We know this is a valid EJA, but will double-check
2062         # if the user passes check_axioms=True.
2063         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2064
2065         from mjo.hurwitz import ComplexMatrixAlgebra
2066         A = ComplexMatrixAlgebra(n, scalars=field)
2067         super().__init__(self._denormalized_basis(A),
2068                          self.jordan_product,
2069                          self.trace_inner_product,
2070                          field=field,
2071                          matrix_space=A,
2072                          **kwargs)
2073
2074         # TODO: this could be factored out somehow, but is left here
2075         # because the MatrixEJA is not presently a subclass of the
2076         # FDEJA class that defines rank() and one().
2077         self.rank.set_cache(n)
2078         self.one.set_cache( self(A.one()) )
2079
2080
2081     @staticmethod
2082     def _max_random_instance_size():
2083         return 3 # Dimension 9
2084
2085     @classmethod
2086     def random_instance(cls, **kwargs):
2087         """
2088         Return a random instance of this type of algebra.
2089         """
2090         n = ZZ.random_element(cls._max_random_instance_size() + 1)
2091         return cls(n, **kwargs)
2092
2093
2094 class QuaternionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
2095     r"""
2096     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
2097     matrices, the usual symmetric Jordan product, and the
2098     real-part-of-trace inner product. It has dimension `2n^2 - n` over
2099     the reals.
2100
2101     SETUP::
2102
2103         sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
2104
2105     EXAMPLES:
2106
2107     In theory, our "field" can be any subfield of the reals::
2108
2109         sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
2110         Euclidean Jordan algebra of dimension 6 over Real Double Field
2111         sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
2112         Euclidean Jordan algebra of dimension 6 over Real Field with
2113         53 bits of precision
2114
2115     TESTS:
2116
2117     The dimension of this algebra is `2*n^2 - n`::
2118
2119         sage: set_random_seed()
2120         sage: n_max = QuaternionHermitianEJA._max_random_instance_size()
2121         sage: n = ZZ.random_element(1, n_max)
2122         sage: J = QuaternionHermitianEJA(n)
2123         sage: J.dimension() == 2*(n^2) - n
2124         True
2125
2126     The Jordan multiplication is what we think it is::
2127
2128         sage: set_random_seed()
2129         sage: J = QuaternionHermitianEJA.random_instance()
2130         sage: x,y = J.random_elements(2)
2131         sage: actual = (x*y).to_matrix()
2132         sage: X = x.to_matrix()
2133         sage: Y = y.to_matrix()
2134         sage: expected = (X*Y + Y*X)/2
2135         sage: actual == expected
2136         True
2137         sage: J(expected) == x*y
2138         True
2139
2140     We can change the generator prefix::
2141
2142         sage: QuaternionHermitianEJA(2, prefix='a').gens()
2143         (a0, a1, a2, a3, a4, a5)
2144
2145     We can construct the (trivial) algebra of rank zero::
2146
2147         sage: QuaternionHermitianEJA(0)
2148         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2149
2150     """
2151     def __init__(self, n, field=AA, **kwargs):
2152         # We know this is a valid EJA, but will double-check
2153         # if the user passes check_axioms=True.
2154         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2155
2156         from mjo.hurwitz import QuaternionMatrixAlgebra
2157         A = QuaternionMatrixAlgebra(n, scalars=field)
2158         super().__init__(self._denormalized_basis(A),
2159                          self.jordan_product,
2160                          self.trace_inner_product,
2161                          field=field,
2162                          matrix_space=A,
2163                          **kwargs)
2164
2165         # TODO: this could be factored out somehow, but is left here
2166         # because the MatrixEJA is not presently a subclass of the
2167         # FDEJA class that defines rank() and one().
2168         self.rank.set_cache(n)
2169         self.one.set_cache( self(A.one()) )
2170
2171
2172     @staticmethod
2173     def _max_random_instance_size():
2174         r"""
2175         The maximum rank of a random QuaternionHermitianEJA.
2176         """
2177         return 2 # Dimension 6
2178
2179     @classmethod
2180     def random_instance(cls, **kwargs):
2181         """
2182         Return a random instance of this type of algebra.
2183         """
2184         n = ZZ.random_element(cls._max_random_instance_size() + 1)
2185         return cls(n, **kwargs)
2186
2187 class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
2188     r"""
2189     SETUP::
2190
2191         sage: from mjo.eja.eja_algebra import (FiniteDimensionalEJA,
2192         ....:                                  OctonionHermitianEJA)
2193         sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra
2194
2195     EXAMPLES:
2196
2197     The 3-by-3 algebra satisfies the axioms of an EJA::
2198
2199         sage: OctonionHermitianEJA(3,                    # long time
2200         ....:                      field=QQ,             # long time
2201         ....:                      orthonormalize=False, # long time
2202         ....:                      check_axioms=True)    # long time
2203         Euclidean Jordan algebra of dimension 27 over Rational Field
2204
2205     After a change-of-basis, the 2-by-2 algebra has the same
2206     multiplication table as the ten-dimensional Jordan spin algebra::
2207
2208         sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ)
2209         sage: b = OctonionHermitianEJA._denormalized_basis(A)
2210         sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],)
2211         sage: jp = OctonionHermitianEJA.jordan_product
2212         sage: ip = OctonionHermitianEJA.trace_inner_product
2213         sage: J = FiniteDimensionalEJA(basis,
2214         ....:                          jp,
2215         ....:                          ip,
2216         ....:                          field=QQ,
2217         ....:                          orthonormalize=False)
2218         sage: J.multiplication_table()
2219         +----++----+----+----+----+----+----+----+----+----+----+
2220         | *  || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2221         +====++====+====+====+====+====+====+====+====+====+====+
2222         | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 |
2223         +----++----+----+----+----+----+----+----+----+----+----+
2224         | b1 || b1 | b0 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | 0  |
2225         +----++----+----+----+----+----+----+----+----+----+----+
2226         | b2 || b2 | 0  | b0 | 0  | 0  | 0  | 0  | 0  | 0  | 0  |
2227         +----++----+----+----+----+----+----+----+----+----+----+
2228         | b3 || b3 | 0  | 0  | b0 | 0  | 0  | 0  | 0  | 0  | 0  |
2229         +----++----+----+----+----+----+----+----+----+----+----+
2230         | b4 || b4 | 0  | 0  | 0  | b0 | 0  | 0  | 0  | 0  | 0  |
2231         +----++----+----+----+----+----+----+----+----+----+----+
2232         | b5 || b5 | 0  | 0  | 0  | 0  | b0 | 0  | 0  | 0  | 0  |
2233         +----++----+----+----+----+----+----+----+----+----+----+
2234         | b6 || b6 | 0  | 0  | 0  | 0  | 0  | b0 | 0  | 0  | 0  |
2235         +----++----+----+----+----+----+----+----+----+----+----+
2236         | b7 || b7 | 0  | 0  | 0  | 0  | 0  | 0  | b0 | 0  | 0  |
2237         +----++----+----+----+----+----+----+----+----+----+----+
2238         | b8 || b8 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | b0 | 0  |
2239         +----++----+----+----+----+----+----+----+----+----+----+
2240         | b9 || b9 | 0  | 0  | 0  | 0  | 0  | 0  | 0  | 0  | b0 |
2241         +----++----+----+----+----+----+----+----+----+----+----+
2242
2243     TESTS:
2244
2245     We can actually construct the 27-dimensional Albert algebra,
2246     and we get the right unit element if we recompute it::
2247
2248         sage: J = OctonionHermitianEJA(3,                    # long time
2249         ....:                          field=QQ,             # long time
2250         ....:                          orthonormalize=False) # long time
2251         sage: J.one.clear_cache()                            # long time
2252         sage: J.one()                                        # long time
2253         b0 + b9 + b26
2254         sage: J.one().to_matrix()                            # long time
2255         +----+----+----+
2256         | e0 | 0  | 0  |
2257         +----+----+----+
2258         | 0  | e0 | 0  |
2259         +----+----+----+
2260         | 0  | 0  | e0 |
2261         +----+----+----+
2262
2263     The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan
2264     spin algebra, but just to be sure, we recompute its rank::
2265
2266         sage: J = OctonionHermitianEJA(2,                    # long time
2267         ....:                          field=QQ,             # long time
2268         ....:                          orthonormalize=False) # long time
2269         sage: J.rank.clear_cache()                           # long time
2270         sage: J.rank()                                       # long time
2271         2
2272
2273     """
2274     @staticmethod
2275     def _max_random_instance_size():
2276         r"""
2277         The maximum rank of a random QuaternionHermitianEJA.
2278         """
2279         return 1 # Dimension 1
2280
2281     @classmethod
2282     def random_instance(cls, **kwargs):
2283         """
2284         Return a random instance of this type of algebra.
2285         """
2286         n = ZZ.random_element(cls._max_random_instance_size() + 1)
2287         return cls(n, **kwargs)
2288
2289     def __init__(self, n, field=AA, **kwargs):
2290         if n > 3:
2291             # Otherwise we don't get an EJA.
2292             raise ValueError("n cannot exceed 3")
2293
2294         # We know this is a valid EJA, but will double-check
2295         # if the user passes check_axioms=True.
2296         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2297
2298         from mjo.hurwitz import OctonionMatrixAlgebra
2299         A = OctonionMatrixAlgebra(n, scalars=field)
2300         super().__init__(self._denormalized_basis(A),
2301                          self.jordan_product,
2302                          self.trace_inner_product,
2303                          field=field,
2304                          matrix_space=A,
2305                          **kwargs)
2306
2307         # TODO: this could be factored out somehow, but is left here
2308         # because the MatrixEJA is not presently a subclass of the
2309         # FDEJA class that defines rank() and one().
2310         self.rank.set_cache(n)
2311         self.one.set_cache( self(A.one()) )
2312
2313
2314 class AlbertEJA(OctonionHermitianEJA):
2315     r"""
2316     The Albert algebra is the algebra of three-by-three Hermitian
2317     matrices whose entries are octonions.
2318
2319     SETUP::
2320
2321         sage: from mjo.eja.eja_algebra import AlbertEJA
2322
2323     EXAMPLES::
2324
2325         sage: AlbertEJA(field=QQ, orthonormalize=False)
2326         Euclidean Jordan algebra of dimension 27 over Rational Field
2327         sage: AlbertEJA() # long time
2328         Euclidean Jordan algebra of dimension 27 over Algebraic Real Field
2329
2330     """
2331     def __init__(self, *args, **kwargs):
2332         super().__init__(3, *args, **kwargs)
2333
2334
2335 class HadamardEJA(RationalBasisEJA, ConcreteEJA):
2336     """
2337     Return the Euclidean Jordan algebra on `R^n` with the Hadamard
2338     (pointwise real-number multiplication) Jordan product and the
2339     usual inner-product.
2340
2341     This is nothing more than the Cartesian product of ``n`` copies of
2342     the one-dimensional Jordan spin algebra, and is the most common
2343     example of a non-simple Euclidean Jordan algebra.
2344
2345     SETUP::
2346
2347         sage: from mjo.eja.eja_algebra import HadamardEJA
2348
2349     EXAMPLES:
2350
2351     This multiplication table can be verified by hand::
2352
2353         sage: J = HadamardEJA(3)
2354         sage: b0,b1,b2 = J.gens()
2355         sage: b0*b0
2356         b0
2357         sage: b0*b1
2358         0
2359         sage: b0*b2
2360         0
2361         sage: b1*b1
2362         b1
2363         sage: b1*b2
2364         0
2365         sage: b2*b2
2366         b2
2367
2368     TESTS:
2369
2370     We can change the generator prefix::
2371
2372         sage: HadamardEJA(3, prefix='r').gens()
2373         (r0, r1, r2)
2374     """
2375     def __init__(self, n, field=AA, **kwargs):
2376         MS = MatrixSpace(field, n, 1)
2377
2378         if n == 0:
2379             jordan_product = lambda x,y: x
2380             inner_product = lambda x,y: x
2381         else:
2382             def jordan_product(x,y):
2383                 return MS( xi*yi for (xi,yi) in zip(x,y) )
2384
2385             def inner_product(x,y):
2386                 return (x.T*y)[0,0]
2387
2388         # New defaults for keyword arguments. Don't orthonormalize
2389         # because our basis is already orthonormal with respect to our
2390         # inner-product. Don't check the axioms, because we know this
2391         # is a valid EJA... but do double-check if the user passes
2392         # check_axioms=True. Note: we DON'T override the "check_field"
2393         # default here, because the user can pass in a field!
2394         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
2395         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2396
2397         column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )
2398         super().__init__(column_basis,
2399                          jordan_product,
2400                          inner_product,
2401                          field=field,
2402                          matrix_space=MS,
2403                          associative=True,
2404                          **kwargs)
2405         self.rank.set_cache(n)
2406
2407         self.one.set_cache( self.sum(self.gens()) )
2408
2409     @staticmethod
2410     def _max_random_instance_size():
2411         r"""
2412         The maximum dimension of a random HadamardEJA.
2413         """
2414         return 5
2415
2416     @classmethod
2417     def random_instance(cls, **kwargs):
2418         """
2419         Return a random instance of this type of algebra.
2420         """
2421         n = ZZ.random_element(cls._max_random_instance_size() + 1)
2422         return cls(n, **kwargs)
2423
2424
2425 class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
2426     r"""
2427     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2428     with the half-trace inner product and jordan product ``x*y =
2429     (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
2430     a symmetric positive-definite "bilinear form" matrix. Its
2431     dimension is the size of `B`, and it has rank two in dimensions
2432     larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
2433     the identity matrix of order ``n``.
2434
2435     We insist that the one-by-one upper-left identity block of `B` be
2436     passed in as well so that we can be passed a matrix of size zero
2437     to construct a trivial algebra.
2438
2439     SETUP::
2440
2441         sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
2442         ....:                                  JordanSpinEJA)
2443
2444     EXAMPLES:
2445
2446     When no bilinear form is specified, the identity matrix is used,
2447     and the resulting algebra is the Jordan spin algebra::
2448
2449         sage: B = matrix.identity(AA,3)
2450         sage: J0 = BilinearFormEJA(B)
2451         sage: J1 = JordanSpinEJA(3)
2452         sage: J0.multiplication_table() == J0.multiplication_table()
2453         True
2454
2455     An error is raised if the matrix `B` does not correspond to a
2456     positive-definite bilinear form::
2457
2458         sage: B = matrix.random(QQ,2,3)
2459         sage: J = BilinearFormEJA(B)
2460         Traceback (most recent call last):
2461         ...
2462         ValueError: bilinear form is not positive-definite
2463         sage: B = matrix.zero(QQ,3)
2464         sage: J = BilinearFormEJA(B)
2465         Traceback (most recent call last):
2466         ...
2467         ValueError: bilinear form is not positive-definite
2468
2469     TESTS:
2470
2471     We can create a zero-dimensional algebra::
2472
2473         sage: B = matrix.identity(AA,0)
2474         sage: J = BilinearFormEJA(B)
2475         sage: J.basis()
2476         Finite family {}
2477
2478     We can check the multiplication condition given in the Jordan, von
2479     Neumann, and Wigner paper (and also discussed on my "On the
2480     symmetry..." paper). Note that this relies heavily on the standard
2481     choice of basis, as does anything utilizing the bilinear form
2482     matrix.  We opt not to orthonormalize the basis, because if we
2483     did, we would have to normalize the `s_{i}` in a similar manner::
2484
2485         sage: set_random_seed()
2486         sage: n = ZZ.random_element(5)
2487         sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
2488         sage: B11 = matrix.identity(QQ,1)
2489         sage: B22 = M.transpose()*M
2490         sage: B = block_matrix(2,2,[ [B11,0  ],
2491         ....:                        [0, B22 ] ])
2492         sage: J = BilinearFormEJA(B, orthonormalize=False)
2493         sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
2494         sage: V = J.vector_space()
2495         sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() )
2496         ....:         for ei in eis ]
2497         sage: actual = [ sis[i]*sis[j]
2498         ....:            for i in range(n-1)
2499         ....:            for j in range(n-1) ]
2500         sage: expected = [ J.one() if i == j else J.zero()
2501         ....:              for i in range(n-1)
2502         ....:              for j in range(n-1) ]
2503         sage: actual == expected
2504         True
2505
2506     """
2507     def __init__(self, B, field=AA, **kwargs):
2508         # The matrix "B" is supplied by the user in most cases,
2509         # so it makes sense to check whether or not its positive-
2510         # definite unless we are specifically asked not to...
2511         if ("check_axioms" not in kwargs) or kwargs["check_axioms"]:
2512             if not B.is_positive_definite():
2513                 raise ValueError("bilinear form is not positive-definite")
2514
2515         # However, all of the other data for this EJA is computed
2516         # by us in manner that guarantees the axioms are
2517         # satisfied. So, again, unless we are specifically asked to
2518         # verify things, we'll skip the rest of the checks.
2519         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2520
2521         n = B.nrows()
2522         MS = MatrixSpace(field, n, 1)
2523
2524         def inner_product(x,y):
2525             return (y.T*B*x)[0,0]
2526
2527         def jordan_product(x,y):
2528             x0 = x[0,0]
2529             xbar = x[1:,0]
2530             y0 = y[0,0]
2531             ybar = y[1:,0]
2532             z0 = inner_product(y,x)
2533             zbar = y0*xbar + x0*ybar
2534             return MS([z0] + zbar.list())
2535
2536         column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() )
2537
2538         # TODO: I haven't actually checked this, but it seems legit.
2539         associative = False
2540         if n <= 2:
2541             associative = True
2542
2543         super().__init__(column_basis,
2544                          jordan_product,
2545                          inner_product,
2546                          field=field,
2547                          matrix_space=MS,
2548                          associative=associative,
2549                          **kwargs)
2550
2551         # The rank of this algebra is two, unless we're in a
2552         # one-dimensional ambient space (because the rank is bounded
2553         # by the ambient dimension).
2554         self.rank.set_cache(min(n,2))
2555         if n == 0:
2556             self.one.set_cache( self.zero() )
2557         else:
2558             self.one.set_cache( self.monomial(0) )
2559
2560     @staticmethod
2561     def _max_random_instance_size():
2562         r"""
2563         The maximum dimension of a random BilinearFormEJA.
2564         """
2565         return 5
2566
2567     @classmethod
2568     def random_instance(cls, **kwargs):
2569         """
2570         Return a random instance of this algebra.
2571         """
2572         n = ZZ.random_element(cls._max_random_instance_size() + 1)
2573         if n.is_zero():
2574             B = matrix.identity(ZZ, n)
2575             return cls(B, **kwargs)
2576
2577         B11 = matrix.identity(ZZ, 1)
2578         M = matrix.random(ZZ, n-1)
2579         I = matrix.identity(ZZ, n-1)
2580         alpha = ZZ.zero()
2581         while alpha.is_zero():
2582             alpha = ZZ.random_element().abs()
2583         B22 = M.transpose()*M + alpha*I
2584
2585         from sage.matrix.special import block_matrix
2586         B = block_matrix(2,2, [ [B11,   ZZ(0) ],
2587                                 [ZZ(0), B22 ] ])
2588
2589         return cls(B, **kwargs)
2590
2591
2592 class JordanSpinEJA(BilinearFormEJA):
2593     """
2594     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
2595     with the usual inner product and jordan product ``x*y =
2596     (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
2597     the reals.
2598
2599     SETUP::
2600
2601         sage: from mjo.eja.eja_algebra import JordanSpinEJA
2602
2603     EXAMPLES:
2604
2605     This multiplication table can be verified by hand::
2606
2607         sage: J = JordanSpinEJA(4)
2608         sage: b0,b1,b2,b3 = J.gens()
2609         sage: b0*b0
2610         b0
2611         sage: b0*b1
2612         b1
2613         sage: b0*b2
2614         b2
2615         sage: b0*b3
2616         b3
2617         sage: b1*b2
2618         0
2619         sage: b1*b3
2620         0
2621         sage: b2*b3
2622         0
2623
2624     We can change the generator prefix::
2625
2626         sage: JordanSpinEJA(2, prefix='B').gens()
2627         (B0, B1)
2628
2629     TESTS:
2630
2631         Ensure that we have the usual inner product on `R^n`::
2632
2633             sage: set_random_seed()
2634             sage: J = JordanSpinEJA.random_instance()
2635             sage: x,y = J.random_elements(2)
2636             sage: actual = x.inner_product(y)
2637             sage: expected = x.to_vector().inner_product(y.to_vector())
2638             sage: actual == expected
2639             True
2640
2641     """
2642     def __init__(self, n, *args, **kwargs):
2643         # This is a special case of the BilinearFormEJA with the
2644         # identity matrix as its bilinear form.
2645         B = matrix.identity(ZZ, n)
2646
2647         # Don't orthonormalize because our basis is already
2648         # orthonormal with respect to our inner-product.
2649         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
2650
2651         # But also don't pass check_field=False here, because the user
2652         # can pass in a field!
2653         super().__init__(B, *args, **kwargs)
2654
2655     @staticmethod
2656     def _max_random_instance_size():
2657         r"""
2658         The maximum dimension of a random JordanSpinEJA.
2659         """
2660         return 5
2661
2662     @classmethod
2663     def random_instance(cls, **kwargs):
2664         """
2665         Return a random instance of this type of algebra.
2666
2667         Needed here to override the implementation for ``BilinearFormEJA``.
2668         """
2669         n = ZZ.random_element(cls._max_random_instance_size() + 1)
2670         return cls(n, **kwargs)
2671
2672
2673 class TrivialEJA(RationalBasisEJA, ConcreteEJA):
2674     """
2675     The trivial Euclidean Jordan algebra consisting of only a zero element.
2676
2677     SETUP::
2678
2679         sage: from mjo.eja.eja_algebra import TrivialEJA
2680
2681     EXAMPLES::
2682
2683         sage: J = TrivialEJA()
2684         sage: J.dimension()
2685         0
2686         sage: J.zero()
2687         0
2688         sage: J.one()
2689         0
2690         sage: 7*J.one()*12*J.one()
2691         0
2692         sage: J.one().inner_product(J.one())
2693         0
2694         sage: J.one().norm()
2695         0
2696         sage: J.one().subalgebra_generated_by()
2697         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
2698         sage: J.rank()
2699         0
2700
2701     """
2702     def __init__(self, field=AA, **kwargs):
2703         jordan_product = lambda x,y: x
2704         inner_product = lambda x,y: field.zero()
2705         basis = ()
2706         MS = MatrixSpace(field,0)
2707
2708         # New defaults for keyword arguments
2709         if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False
2710         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
2711
2712         super().__init__(basis,
2713                          jordan_product,
2714                          inner_product,
2715                          associative=True,
2716                          field=field,
2717                          matrix_space=MS,
2718                          **kwargs)
2719
2720         # The rank is zero using my definition, namely the dimension of the
2721         # largest subalgebra generated by any element.
2722         self.rank.set_cache(0)
2723         self.one.set_cache( self.zero() )
2724
2725     @classmethod
2726     def random_instance(cls, **kwargs):
2727         # We don't take a "size" argument so the superclass method is
2728         # inappropriate for us.
2729         return cls(**kwargs)
2730
2731
2732 class CartesianProductEJA(FiniteDimensionalEJA):
2733     r"""
2734     The external (orthogonal) direct sum of two or more Euclidean
2735     Jordan algebras. Every Euclidean Jordan algebra decomposes into an
2736     orthogonal direct sum of simple Euclidean Jordan algebras which is
2737     then isometric to a Cartesian product, so no generality is lost by
2738     providing only this construction.
2739
2740     SETUP::
2741
2742         sage: from mjo.eja.eja_algebra import (random_eja,
2743         ....:                                  CartesianProductEJA,
2744         ....:                                  HadamardEJA,
2745         ....:                                  JordanSpinEJA,
2746         ....:                                  RealSymmetricEJA)
2747
2748     EXAMPLES:
2749
2750     The Jordan product is inherited from our factors and implemented by
2751     our CombinatorialFreeModule Cartesian product superclass::
2752
2753         sage: set_random_seed()
2754         sage: J1 = HadamardEJA(2)
2755         sage: J2 = RealSymmetricEJA(2)
2756         sage: J = cartesian_product([J1,J2])
2757         sage: x,y = J.random_elements(2)
2758         sage: x*y in J
2759         True
2760
2761     The ability to retrieve the original factors is implemented by our
2762     CombinatorialFreeModule Cartesian product superclass::
2763
2764         sage: J1 = HadamardEJA(2, field=QQ)
2765         sage: J2 = JordanSpinEJA(3, field=QQ)
2766         sage: J = cartesian_product([J1,J2])
2767         sage: J.cartesian_factors()
2768         (Euclidean Jordan algebra of dimension 2 over Rational Field,
2769          Euclidean Jordan algebra of dimension 3 over Rational Field)
2770
2771     You can provide more than two factors::
2772
2773         sage: J1 = HadamardEJA(2)
2774         sage: J2 = JordanSpinEJA(3)
2775         sage: J3 = RealSymmetricEJA(3)
2776         sage: cartesian_product([J1,J2,J3])
2777         Euclidean Jordan algebra of dimension 2 over Algebraic Real
2778         Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic
2779         Real Field (+) Euclidean Jordan algebra of dimension 6 over
2780         Algebraic Real Field
2781
2782     Rank is additive on a Cartesian product::
2783
2784         sage: J1 = HadamardEJA(1)
2785         sage: J2 = RealSymmetricEJA(2)
2786         sage: J = cartesian_product([J1,J2])
2787         sage: J1.rank.clear_cache()
2788         sage: J2.rank.clear_cache()
2789         sage: J.rank.clear_cache()
2790         sage: J.rank()
2791         3
2792         sage: J.rank() == J1.rank() + J2.rank()
2793         True
2794
2795     The same rank computation works over the rationals, with whatever
2796     basis you like::
2797
2798         sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
2799         sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
2800         sage: J = cartesian_product([J1,J2])
2801         sage: J1.rank.clear_cache()
2802         sage: J2.rank.clear_cache()
2803         sage: J.rank.clear_cache()
2804         sage: J.rank()
2805         3
2806         sage: J.rank() == J1.rank() + J2.rank()
2807         True
2808
2809     The product algebra will be associative if and only if all of its
2810     components are associative::
2811
2812         sage: J1 = HadamardEJA(2)
2813         sage: J1.is_associative()
2814         True
2815         sage: J2 = HadamardEJA(3)
2816         sage: J2.is_associative()
2817         True
2818         sage: J3 = RealSymmetricEJA(3)
2819         sage: J3.is_associative()
2820         False
2821         sage: CP1 = cartesian_product([J1,J2])
2822         sage: CP1.is_associative()
2823         True
2824         sage: CP2 = cartesian_product([J1,J3])
2825         sage: CP2.is_associative()
2826         False
2827
2828     Cartesian products of Cartesian products work::
2829
2830         sage: J1 = JordanSpinEJA(1)
2831         sage: J2 = JordanSpinEJA(1)
2832         sage: J3 = JordanSpinEJA(1)
2833         sage: J = cartesian_product([J1,cartesian_product([J2,J3])])
2834         sage: J.multiplication_table()
2835         +----++----+----+----+
2836         | *  || b0 | b1 | b2 |
2837         +====++====+====+====+
2838         | b0 || b0 | 0  | 0  |
2839         +----++----+----+----+
2840         | b1 || 0  | b1 | 0  |
2841         +----++----+----+----+
2842         | b2 || 0  | 0  | b2 |
2843         +----++----+----+----+
2844         sage: HadamardEJA(3).multiplication_table()
2845         +----++----+----+----+
2846         | *  || b0 | b1 | b2 |
2847         +====++====+====+====+
2848         | b0 || b0 | 0  | 0  |
2849         +----++----+----+----+
2850         | b1 || 0  | b1 | 0  |
2851         +----++----+----+----+
2852         | b2 || 0  | 0  | b2 |
2853         +----++----+----+----+
2854
2855     TESTS:
2856
2857     All factors must share the same base field::
2858
2859         sage: J1 = HadamardEJA(2, field=QQ)
2860         sage: J2 = RealSymmetricEJA(2)
2861         sage: CartesianProductEJA((J1,J2))
2862         Traceback (most recent call last):
2863         ...
2864         ValueError: all factors must share the same base field
2865
2866     The cached unit element is the same one that would be computed::
2867
2868         sage: set_random_seed()              # long time
2869         sage: J1 = random_eja()              # long time
2870         sage: J2 = random_eja()              # long time
2871         sage: J = cartesian_product([J1,J2]) # long time
2872         sage: actual = J.one()               # long time
2873         sage: J.one.clear_cache()            # long time
2874         sage: expected = J.one()             # long time
2875         sage: actual == expected             # long time
2876         True
2877
2878     """
2879     Element = FiniteDimensionalEJAElement
2880
2881
2882     def __init__(self, factors, **kwargs):
2883         m = len(factors)
2884         if m == 0:
2885             return TrivialEJA()
2886
2887         self._sets = factors
2888
2889         field = factors[0].base_ring()
2890         if not all( J.base_ring() == field for J in factors ):
2891             raise ValueError("all factors must share the same base field")
2892
2893         associative = all( f.is_associative() for f in factors )
2894
2895         # Compute my matrix space. This category isn't perfect, but
2896         # is good enough for what we need to do.
2897         MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
2898         MS_cat = MS_cat.Unital().CartesianProducts()
2899         MS_factors = tuple( J.matrix_space() for J in factors )
2900         from sage.sets.cartesian_product import CartesianProduct
2901         MS = CartesianProduct(MS_factors, MS_cat)
2902
2903         basis = []
2904         zero = MS.zero()
2905         for i in range(m):
2906             for b in factors[i].matrix_basis():
2907                 z = list(zero)
2908                 z[i] = b
2909                 basis.append(z)
2910
2911         basis = tuple( MS(b) for b in basis )
2912
2913         # Define jordan/inner products that operate on that matrix_basis.
2914         def jordan_product(x,y):
2915             return MS(tuple(
2916                 (factors[i](x[i])*factors[i](y[i])).to_matrix()
2917                 for i in range(m)
2918             ))
2919
2920         def inner_product(x, y):
2921             return sum(
2922                 factors[i](x[i]).inner_product(factors[i](y[i]))
2923                 for i in range(m)
2924             )
2925
2926         # There's no need to check the field since it already came
2927         # from an EJA. Likewise the axioms are guaranteed to be
2928         # satisfied, unless the guy writing this class sucks.
2929         #
2930         # If you want the basis to be orthonormalized, orthonormalize
2931         # the factors.
2932         FiniteDimensionalEJA.__init__(self,
2933                                       basis,
2934                                       jordan_product,
2935                                       inner_product,
2936                                       field=field,
2937                                       matrix_space=MS,
2938                                       orthonormalize=False,
2939                                       associative=associative,
2940                                       cartesian_product=True,
2941                                       check_field=False,
2942                                       check_axioms=False)
2943
2944         self.rank.set_cache(sum(J.rank() for J in factors))
2945         ones = tuple(J.one().to_matrix() for J in factors)
2946         self.one.set_cache(self(ones))
2947
2948     def cartesian_factors(self):
2949         # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
2950         return self._sets
2951
2952     def cartesian_factor(self, i):
2953         r"""
2954         Return the ``i``th factor of this algebra.
2955         """
2956         return self._sets[i]
2957
2958     def _repr_(self):
2959         # Copy/pasted from CombinatorialFreeModule_CartesianProduct.
2960         from sage.categories.cartesian_product import cartesian_product
2961         return cartesian_product.symbol.join("%s" % factor
2962                                              for factor in self._sets)
2963
2964     def matrix_space(self):
2965         r"""
2966         Return the space that our matrix basis lives in as a Cartesian
2967         product.
2968
2969         We don't simply use the ``cartesian_product()`` functor here
2970         because it acts differently on SageMath MatrixSpaces and our
2971         custom MatrixAlgebras, which are CombinatorialFreeModules. We
2972         always want the result to be represented (and indexed) as
2973         an ordered tuple.
2974
2975         SETUP::
2976
2977             sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
2978             ....:                                  HadamardEJA,
2979             ....:                                  OctonionHermitianEJA,
2980             ....:                                  RealSymmetricEJA)
2981
2982         EXAMPLES::
2983
2984             sage: J1 = HadamardEJA(1)
2985             sage: J2 = RealSymmetricEJA(2)
2986             sage: J = cartesian_product([J1,J2])
2987             sage: J.matrix_space()
2988             The Cartesian product of (Full MatrixSpace of 1 by 1 dense
2989             matrices over Algebraic Real Field, Full MatrixSpace of 2
2990             by 2 dense matrices over Algebraic Real Field)
2991
2992         ::
2993
2994             sage: J1 = ComplexHermitianEJA(1)
2995             sage: J2 = ComplexHermitianEJA(1)
2996             sage: J = cartesian_product([J1,J2])
2997             sage: J.one().to_matrix()[0]
2998             +---+
2999             | 1 |
3000             +---+
3001             sage: J.one().to_matrix()[1]
3002             +---+
3003             | 1 |
3004             +---+
3005
3006         ::
3007
3008             sage: J1 = OctonionHermitianEJA(1)
3009             sage: J2 = OctonionHermitianEJA(1)
3010             sage: J = cartesian_product([J1,J2])
3011             sage: J.one().to_matrix()[0]
3012             +----+
3013             | e0 |
3014             +----+
3015             sage: J.one().to_matrix()[1]
3016             +----+
3017             | e0 |
3018             +----+
3019
3020         """
3021         return super().matrix_space()
3022
3023
3024     @cached_method
3025     def cartesian_projection(self, i):
3026         r"""
3027         SETUP::
3028
3029             sage: from mjo.eja.eja_algebra import (random_eja,
3030             ....:                                  JordanSpinEJA,
3031             ....:                                  HadamardEJA,
3032             ....:                                  RealSymmetricEJA,
3033             ....:                                  ComplexHermitianEJA)
3034
3035         EXAMPLES:
3036
3037         The projection morphisms are Euclidean Jordan algebra
3038         operators::
3039
3040             sage: J1 = HadamardEJA(2)
3041             sage: J2 = RealSymmetricEJA(2)
3042             sage: J = cartesian_product([J1,J2])
3043             sage: J.cartesian_projection(0)
3044             Linear operator between finite-dimensional Euclidean Jordan
3045             algebras represented by the matrix:
3046             [1 0 0 0 0]
3047             [0 1 0 0 0]
3048             Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3049             Real Field (+) Euclidean Jordan algebra of dimension 3 over
3050             Algebraic Real Field
3051             Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic
3052             Real Field
3053             sage: J.cartesian_projection(1)
3054             Linear operator between finite-dimensional Euclidean Jordan
3055             algebras represented by the matrix:
3056             [0 0 1 0 0]
3057             [0 0 0 1 0]
3058             [0 0 0 0 1]
3059             Domain: Euclidean Jordan algebra of dimension 2 over Algebraic
3060             Real Field (+) Euclidean Jordan algebra of dimension 3 over
3061             Algebraic Real Field
3062             Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic
3063             Real Field
3064
3065         The projections work the way you'd expect on the vector
3066         representation of an element::
3067
3068             sage: J1 = JordanSpinEJA(2)
3069             sage: J2 = ComplexHermitianEJA(2)
3070             sage: J = cartesian_product([J1,J2])
3071             sage: pi_left = J.cartesian_projection(0)
3072             sage: pi_right = J.cartesian_projection(1)
3073             sage: pi_left(J.one()).to_vector()
3074             (1, 0)
3075             sage: pi_right(J.one()).to_vector()
3076             (1, 0, 0, 1)
3077             sage: J.one().to_vector()
3078             (1, 0, 1, 0, 0, 1)
3079
3080         TESTS:
3081
3082         The answer never changes::
3083
3084             sage: set_random_seed()
3085             sage: J1 = random_eja()
3086             sage: J2 = random_eja()
3087             sage: J = cartesian_product([J1,J2])
3088             sage: P0 = J.cartesian_projection(0)
3089             sage: P1 = J.cartesian_projection(0)
3090             sage: P0 == P1
3091             True
3092
3093         """
3094         offset = sum( self.cartesian_factor(k).dimension()
3095                       for k in range(i) )
3096         Ji = self.cartesian_factor(i)
3097         Pi = self._module_morphism(lambda j: Ji.monomial(j - offset),
3098                                    codomain=Ji)
3099
3100         return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix())
3101
3102     @cached_method
3103     def cartesian_embedding(self, i):
3104         r"""
3105         SETUP::
3106
3107             sage: from mjo.eja.eja_algebra import (random_eja,
3108             ....:                                  JordanSpinEJA,
3109             ....:                                  HadamardEJA,
3110             ....:                                  RealSymmetricEJA)
3111
3112         EXAMPLES:
3113
3114         The embedding morphisms are Euclidean Jordan algebra
3115         operators::
3116
3117             sage: J1 = HadamardEJA(2)
3118             sage: J2 = RealSymmetricEJA(2)
3119             sage: J = cartesian_product([J1,J2])
3120             sage: J.cartesian_embedding(0)
3121             Linear operator between finite-dimensional Euclidean Jordan
3122             algebras represented by the matrix:
3123             [1 0]
3124             [0 1]
3125             [0 0]
3126             [0 0]
3127             [0 0]
3128             Domain: Euclidean Jordan algebra of dimension 2 over
3129             Algebraic Real Field
3130             Codomain: Euclidean Jordan algebra of dimension 2 over
3131             Algebraic Real Field (+) Euclidean Jordan algebra of
3132             dimension 3 over Algebraic Real Field
3133             sage: J.cartesian_embedding(1)
3134             Linear operator between finite-dimensional Euclidean Jordan
3135             algebras represented by the matrix:
3136             [0 0 0]
3137             [0 0 0]
3138             [1 0 0]
3139             [0 1 0]
3140             [0 0 1]
3141             Domain: Euclidean Jordan algebra of dimension 3 over
3142             Algebraic Real Field
3143             Codomain: Euclidean Jordan algebra of dimension 2 over
3144             Algebraic Real Field (+) Euclidean Jordan algebra of
3145             dimension 3 over Algebraic Real Field
3146
3147         The embeddings work the way you'd expect on the vector
3148         representation of an element::
3149
3150             sage: J1 = JordanSpinEJA(3)
3151             sage: J2 = RealSymmetricEJA(2)
3152             sage: J = cartesian_product([J1,J2])
3153             sage: iota_left = J.cartesian_embedding(0)
3154             sage: iota_right = J.cartesian_embedding(1)
3155             sage: iota_left(J1.zero()) == J.zero()
3156             True
3157             sage: iota_right(J2.zero()) == J.zero()
3158             True
3159             sage: J1.one().to_vector()
3160             (1, 0, 0)
3161             sage: iota_left(J1.one()).to_vector()
3162             (1, 0, 0, 0, 0, 0)
3163             sage: J2.one().to_vector()
3164             (1, 0, 1)
3165             sage: iota_right(J2.one()).to_vector()
3166             (0, 0, 0, 1, 0, 1)
3167             sage: J.one().to_vector()
3168             (1, 0, 0, 1, 0, 1)
3169
3170         TESTS:
3171
3172         The answer never changes::
3173
3174             sage: set_random_seed()
3175             sage: J1 = random_eja()
3176             sage: J2 = random_eja()
3177             sage: J = cartesian_product([J1,J2])
3178             sage: E0 = J.cartesian_embedding(0)
3179             sage: E1 = J.cartesian_embedding(0)
3180             sage: E0 == E1
3181             True
3182
3183         Composing a projection with the corresponding inclusion should
3184         produce the identity map, and mismatching them should produce
3185         the zero map::
3186
3187             sage: set_random_seed()
3188             sage: J1 = random_eja()
3189             sage: J2 = random_eja()
3190             sage: J = cartesian_product([J1,J2])
3191             sage: iota_left = J.cartesian_embedding(0)
3192             sage: iota_right = J.cartesian_embedding(1)
3193             sage: pi_left = J.cartesian_projection(0)
3194             sage: pi_right = J.cartesian_projection(1)
3195             sage: pi_left*iota_left == J1.one().operator()
3196             True
3197             sage: pi_right*iota_right == J2.one().operator()
3198             True
3199             sage: (pi_left*iota_right).is_zero()
3200             True
3201             sage: (pi_right*iota_left).is_zero()
3202             True
3203
3204         """
3205         offset = sum( self.cartesian_factor(k).dimension()
3206                       for k in range(i) )
3207         Ji = self.cartesian_factor(i)
3208         Ei = Ji._module_morphism(lambda j: self.monomial(j + offset),
3209                                  codomain=self)
3210         return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
3211
3212
3213
3214 FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
3215
3216 class RationalBasisCartesianProductEJA(CartesianProductEJA,
3217                                        RationalBasisEJA):
3218     r"""
3219     A separate class for products of algebras for which we know a
3220     rational basis.
3221
3222     SETUP::
3223
3224         sage: from mjo.eja.eja_algebra import (HadamardEJA,
3225         ....:                                  JordanSpinEJA,
3226         ....:                                  OctonionHermitianEJA,
3227         ....:                                  RealSymmetricEJA)
3228
3229     EXAMPLES:
3230
3231     This gives us fast characteristic polynomial computations in
3232     product algebras, too::
3233
3234
3235         sage: J1 = JordanSpinEJA(2)
3236         sage: J2 = RealSymmetricEJA(3)
3237         sage: J = cartesian_product([J1,J2])
3238         sage: J.characteristic_polynomial_of().degree()
3239         5
3240         sage: J.rank()
3241         5
3242
3243     TESTS:
3244
3245     The ``cartesian_product()`` function only uses the first factor to
3246     decide where the result will live; thus we have to be careful to
3247     check that all factors do indeed have a `_rational_algebra` member
3248     before we try to access it::
3249
3250         sage: J1 = OctonionHermitianEJA(1) # no rational basis
3251         sage: J2 = HadamardEJA(2)
3252         sage: cartesian_product([J1,J2])
3253         Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3254         (+) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3255         sage: cartesian_product([J2,J1])
3256         Euclidean Jordan algebra of dimension 2 over Algebraic Real Field
3257         (+) Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
3258
3259     """
3260     def __init__(self, algebras, **kwargs):
3261         CartesianProductEJA.__init__(self, algebras, **kwargs)
3262
3263         self._rational_algebra = None
3264         if self.vector_space().base_field() is not QQ:
3265             if all( hasattr(r, "_rational_algebra") for r in algebras ):
3266                 self._rational_algebra = cartesian_product([
3267                     r._rational_algebra for r in algebras
3268                 ])
3269
3270
3271 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
3272
3273 def random_eja(*args, **kwargs):
3274     J1 = ConcreteEJA.random_instance(*args, **kwargs)
3275
3276     # This might make Cartesian products appear roughly as often as
3277     # any other ConcreteEJA.
3278     if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
3279         # Use random_eja() again so we can get more than two factors.
3280         J2 = random_eja(*args, **kwargs)
3281         J = cartesian_product([J1,J2])
3282         return J
3283     else:
3284         return J1