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