X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=d3eac4f6d3bfbad50f6bbc4b371aaa9d39f8859b;hb=HEAD;hp=70be6ed43bb8d091d6f2b4565400f368bd75b741;hpb=5afa0afc19e6a4695687ea9f0c26fc531e18e5b8;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 70be6ed..adcc343 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -1,9 +1,149 @@ -""" -Euclidean Jordan Algebras. These are formally-real Jordan Algebras; -specifically those where u^2 + v^2 = 0 implies that u = v = 0. They -are used in optimization, and have some additional nice methods beyond -what can be supported in a general Jordan Algebra. - +r""" +Representations and constructions for Euclidean Jordan algebras. + +A Euclidean Jordan algebra is a Jordan algebra that has some +additional properties: + + 1. It is finite-dimensional. + 2. Its scalar field is the real numbers. + 3a. An inner product is defined on it, and... + 3b. That inner product is compatible with the Jordan product + in the sense that ` = ` for all elements + `x,y,z` in the algebra. + +Every Euclidean Jordan algebra is formally-real: for any two elements +`x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y = +0`. Conversely, every finite-dimensional formally-real Jordan algebra +can be made into a Euclidean Jordan algebra with an appropriate choice +of inner-product. + +Formally-real Jordan algebras were originally studied as a framework +for quantum mechanics. Today, Euclidean Jordan algebras are crucial in +symmetric cone optimization, since every symmetric cone arises as the +cone of squares in some Euclidean Jordan algebra. + +It is known that every Euclidean Jordan algebra decomposes into an +orthogonal direct sum (essentially, a Cartesian product) of simple +algebras, and that moreover, up to Jordan-algebra isomorphism, there +are only five families of simple algebras. We provide constructions +for these simple algebras: + + * :class:`BilinearFormEJA` + * :class:`RealSymmetricEJA` + * :class:`ComplexHermitianEJA` + * :class:`QuaternionHermitianEJA` + * :class:`OctonionHermitianEJA` + +In addition to these, we provide a few other example constructions, + + * :class:`JordanSpinEJA` + * :class:`HadamardEJA` + * :class:`AlbertEJA` + * :class:`TrivialEJA` + * :class:`ComplexSkewSymmetricEJA` + +The Jordan spin algebra is a bilinear form algebra where the bilinear +form is the identity. The Hadamard EJA is simply a Cartesian product +of one-dimensional spin algebras. The Albert EJA is simply a special +case of the :class:`OctonionHermitianEJA` where the matrices are +three-by-three and the resulting space has dimension 27. And +last/least, the trivial EJA is exactly what you think it is; it could +also be obtained by constructing a dimension-zero instance of any of +the other algebras. Cartesian products of these are also supported +using the usual ``cartesian_product()`` function; as a result, we +support (up to isomorphism) all Euclidean Jordan algebras. + +At a minimum, the following are required to construct a Euclidean +Jordan algebra: + + * A basis of matrices, column vectors, or MatrixAlgebra elements + * A Jordan product defined on the basis + * Its inner product defined on the basis + +The real numbers form a Euclidean Jordan algebra when both the Jordan +and inner products are the usual multiplication. We use this as our +example, and demonstrate a few ways to construct an EJA. + +First, we can use one-by-one SageMath matrices with algebraic real +entries to represent real numbers. We define the Jordan and inner +products to be essentially real-number multiplication, with the only +difference being that the Jordan product again returns a one-by-one +matrix, whereas the inner product must return a scalar. Our basis for +the one-by-one matrices is of course the set consisting of a single +matrix with its sole entry non-zero:: + + sage: from mjo.eja.eja_algebra import EJA + sage: jp = lambda X,Y: X*Y + sage: ip = lambda X,Y: X[0,0]*Y[0,0] + sage: b1 = matrix(AA, [[1]]) + sage: J1 = EJA((b1,), jp, ip) + sage: J1 + Euclidean Jordan algebra of dimension 1 over Algebraic Real Field + +In fact, any positive scalar multiple of that inner-product would work:: + + sage: ip2 = lambda X,Y: 16*ip(X,Y) + sage: J2 = EJA((b1,), jp, ip2) + sage: J2 + Euclidean Jordan algebra of dimension 1 over Algebraic Real Field + +But beware that your basis will be orthonormalized _with respect to the +given inner-product_ unless you pass ``orthonormalize=False`` to the +constructor. For example:: + + sage: J3 = EJA((b1,), jp, ip2, orthonormalize=False) + sage: J3 + Euclidean Jordan algebra of dimension 1 over Algebraic Real Field + +To see the difference, you can take the first and only basis element +of the resulting algebra, and ask for it to be converted back into +matrix form:: + + sage: J1.basis()[0].to_matrix() + [1] + sage: J2.basis()[0].to_matrix() + [1/4] + sage: J3.basis()[0].to_matrix() + [1] + +Since square roots are used in that process, the default scalar field +that we use is the field of algebraic real numbers, ``AA``. You can +also Use rational numbers, but only if you either pass +``orthonormalize=False`` or know that orthonormalizing your basis +won't stray beyond the rational numbers. The example above would +have worked only because ``sqrt(16) == 4`` is rational. + +Another option for your basis is to use elemebts of a +:class:`MatrixAlgebra`:: + + sage: from mjo.matrix_algebra import MatrixAlgebra + sage: A = MatrixAlgebra(1,AA,AA) + sage: J4 = EJA(A.gens(), jp, ip) + sage: J4 + Euclidean Jordan algebra of dimension 1 over Algebraic Real Field + sage: J4.basis()[0].to_matrix() + +---+ + | 1 | + +---+ + +An easier way to view the entire EJA basis in its original (but +perhaps orthonormalized) matrix form is to use the ``matrix_basis`` +method:: + + sage: J4.matrix_basis() + (+---+ + | 1 | + +---+,) + +In particular, a :class:`MatrixAlgebra` is needed to work around the +fact that matrices in SageMath must have entries in the same +(commutative and associative) ring as its scalars. There are many +Euclidean Jordan algebras whose elements are matrices that violate +those assumptions. The complex, quaternion, and octonion Hermitian +matrices all have entries in a ring (the complex numbers, quaternions, +or octonions...) that differs from the algebra's scalar ring (the real +numbers). Quaternions are also non-commutative; the octonions are +neither commutative nor associative. SETUP:: @@ -13,13 +153,11 @@ EXAMPLES:: sage: random_eja() Euclidean Jordan algebra of dimension... - """ -from itertools import repeat - from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra from sage.categories.magmatic_algebras import MagmaticAlgebras +from sage.categories.sets_cat import cartesian_product from sage.combinat.free_module import CombinatorialFreeModule from sage.matrix.constructor import matrix from sage.matrix.matrix_space import MatrixSpace @@ -29,14 +167,257 @@ from sage.modules.free_module import FreeModule, VectorSpace from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF, PolynomialRing, QuadraticField) -from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement -from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator -from mjo.eja.eja_utils import _mat2vec +from mjo.eja.eja_element import (CartesianProductEJAElement, + EJAElement) +from mjo.eja.eja_operator import EJAOperator +from mjo.eja.eja_utils import _all2list + +def EuclideanJordanAlgebras(field): + r""" + The category of Euclidean Jordan algebras over ``field``, which + must be a subfield of the real numbers. For now this is just a + convenient wrapper around all of the other category axioms that + apply to all EJAs. + """ + category = MagmaticAlgebras(field).FiniteDimensional() + category = category.WithBasis().Unital().Commutative() + return category -class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): +class EJA(CombinatorialFreeModule): r""" - The lowest-level class for representing a Euclidean Jordan algebra. + A finite-dimensional Euclidean Jordan algebra. + + INPUT: + + - ``basis`` -- a tuple; a tuple of basis elements in "matrix + form," which must be the same form as the arguments to + ``jordan_product`` and ``inner_product``. In reality, "matrix + form" can be either vectors, matrices, or a Cartesian product + (ordered tuple) of vectors or matrices. All of these would + ideally be vector spaces in sage with no special-casing + needed; but in reality we turn vectors into column-matrices + and Cartesian products `(a,b)` into column matrices + `(a,b)^{T}` after converting `a` and `b` themselves. + + - ``jordan_product`` -- a function; afunction of two ``basis`` + elements (in matrix form) that returns their jordan product, + also in matrix form; this will be applied to ``basis`` to + compute a multiplication table for the algebra. + + - ``inner_product`` -- a function; a function of two ``basis`` + elements (in matrix form) that returns their inner + product. This will be applied to ``basis`` to compute an + inner-product table (basically a matrix) for this algebra. + + - ``matrix_space`` -- the space that your matrix basis lives in, + or ``None`` (the default). So long as your basis does not have + length zero you can omit this. But in trivial algebras, it is + required. + + - ``field`` -- a subfield of the reals (default: ``AA``); the scalar + field for the algebra. + + - ``orthonormalize`` -- boolean (default: ``True``); whether or + not to orthonormalize the basis. Doing so is expensive and + generally rules out using the rationals as your ``field``, but + is required for spectral decompositions. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS: + + We should compute that an element subalgebra is associative even + if we circumvent the element method:: + + sage: J = random_eja(field=QQ,orthonormalize=False) + sage: x = J.random_element() + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: basis = tuple(b.superalgebra_element() for b in A.basis()) + sage: J.subalgebra(basis, orthonormalize=False).is_associative() + True """ + Element = EJAElement + + @staticmethod + def _check_input_field(field): + if not field.is_subring(RR): + # Note: this does return true for the real algebraic + # field, the rationals, and any quadratic field where + # we've specified a real embedding. + raise ValueError("scalar field is not real") + + @staticmethod + def _check_input_axioms(basis, jordan_product, inner_product): + if not all( jordan_product(bi,bj) == jordan_product(bj,bi) + for bi in basis + for bj in basis ): + raise ValueError("Jordan product is not commutative") + + if not all( inner_product(bi,bj) == inner_product(bj,bi) + for bi in basis + for bj in basis ): + raise ValueError("inner-product is not commutative") + + def __init__(self, + basis, + jordan_product, + inner_product, + field=AA, + matrix_space=None, + orthonormalize=True, + associative=None, + check_field=True, + check_axioms=True, + prefix="b"): + + n = len(basis) + + if check_field: + self._check_input_field(field) + + if check_axioms: + # Check commutativity of the Jordan and inner-products. + # This has to be done before we build the multiplication + # and inner-product tables/matrices, because we take + # advantage of symmetry in the process. + self._check_input_axioms(basis, jordan_product, inner_product) + + if n <= 1: + # All zero- and one-dimensional algebras are just the real + # numbers with (some positive multiples of) the usual + # multiplication as its Jordan and inner-product. + associative = True + if associative is None: + # We should figure it out. As with check_axioms, we have to do + # this without the help of the _jordan_product_is_associative() + # method because we need to know the category before we + # initialize the algebra. + associative = all( jordan_product(jordan_product(bi,bj),bk) + == + jordan_product(bi,jordan_product(bj,bk)) + for bi in basis + for bj in basis + for bk in basis) + + category = EuclideanJordanAlgebras(field) + + if associative: + # Element subalgebras can take advantage of this. + category = category.Associative() + + # Call the superclass constructor so that we can use its from_vector() + # method to build our multiplication table. + CombinatorialFreeModule.__init__(self, + field, + range(n), + prefix=prefix, + category=category, + bracket=False) + + # Now comes all of the hard work. We'll be constructing an + # ambient vector space V that our (vectorized) basis lives in, + # as well as a subspace W of V spanned by those (vectorized) + # basis elements. The W-coordinates are the coefficients that + # we see in things like x = 1*b1 + 2*b2. + + degree = 0 + if n > 0: + degree = len(_all2list(basis[0])) + + # Build an ambient space that fits our matrix basis when + # written out as "long vectors." + V = VectorSpace(field, degree) + + # The matrix that will hold the orthonormal -> unorthonormal + # coordinate transformation. Default to an identity matrix of + # the appropriate size to avoid special cases for None + # everywhere. + self._deortho_matrix = matrix.identity(field,n) + + if orthonormalize: + # Save a copy of the un-orthonormalized basis for later. + # Convert it to ambient V (vector) coordinates while we're + # at it, because we'd have to do it later anyway. + deortho_vector_basis = tuple( V(_all2list(b)) for b in basis ) + + from mjo.eja.eja_utils import gram_schmidt + basis = tuple(gram_schmidt(basis, inner_product)) + + # Save the (possibly orthonormalized) matrix basis for + # later, as well as the space that its elements live in. + # In most cases we can deduce the matrix space, but when + # n == 0 (that is, there are no basis elements) we cannot. + self._matrix_basis = basis + if matrix_space is None: + self._matrix_space = self._matrix_basis[0].parent() + else: + self._matrix_space = matrix_space + + # Now create the vector space for the algebra, which will have + # its own set of non-ambient coordinates (in terms of the + # supplied basis). + vector_basis = tuple( V(_all2list(b)) for b in basis ) + + # Save the span of our matrix basis (when written out as long + # vectors) because otherwise we'll have to reconstruct it + # every time we want to coerce a matrix into the algebra. + self._matrix_span = V.span_of_basis( vector_basis, check=check_axioms) + + if orthonormalize: + # Now "self._matrix_span" is the vector space of our + # algebra coordinates. The variables "X0", "X1",... refer + # to the entries of vectors in self._matrix_span. Thus to + # convert back and forth between the orthonormal + # coordinates and the given ones, we need to stick the + # original basis in self._matrix_span. + U = V.span_of_basis( deortho_vector_basis, check=check_axioms) + self._deortho_matrix = matrix.column( U.coordinate_vector(q) + for q in vector_basis ) + + + # Now we actually compute the multiplication and inner-product + # tables/matrices using the possibly-orthonormalized basis. + self._inner_product_matrix = matrix.identity(field, n) + zed = self.zero() + self._multiplication_table = [ [zed for j in range(i+1)] + for i in range(n) ] + + # Note: the Jordan and inner-products are defined in terms + # of the ambient basis. It's important that their arguments + # are in ambient coordinates as well. + for i in range(n): + for j in range(i+1): + # ortho basis w.r.t. ambient coords + q_i = basis[i] + q_j = basis[j] + + # The jordan product returns a matrixy answer, so we + # have to convert it to the algebra coordinates. + elt = jordan_product(q_i, q_j) + elt = self._matrix_span.coordinate_vector(V(_all2list(elt))) + self._multiplication_table[i][j] = self.from_vector(elt) + + if not orthonormalize: + # If we're orthonormalizing the basis with respect + # to an inner-product, then the inner-product + # matrix with respect to the resulting basis is + # just going to be the identity. + ip = inner_product(q_i, q_j) + self._inner_product_matrix[i,j] = ip + self._inner_product_matrix[j,i] = ip + + self._inner_product_matrix._cache = {'hermitian': True} + self._inner_product_matrix.set_immutable() + + if check_axioms: + if not self._is_jordanian(): + raise ValueError("Jordan identity does not hold") + if not self._inner_product_is_associative(): + raise ValueError("inner product is not associative") + + def _coerce_map_from_base_ring(self): """ Disable the map from the base ring into the algebra. @@ -51,7 +432,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): TESTS:: - sage: set_random_seed() sage: J = random_eja() sage: J(1) Traceback (most recent call last): @@ -61,200 +441,272 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): """ return None - def __init__(self, - field, - multiplication_table, - inner_product_table, - prefix='e', - category=None, - matrix_basis=None, - check_field=True, - check_axioms=True): + + def product_on_basis(self, i, j): + r""" + Returns the Jordan product of the `i` and `j`th basis elements. + + This completely defines the Jordan product on the algebra, and + is used direclty by our superclass machinery to implement + :meth:`product`. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS:: + + sage: J = random_eja() + sage: n = J.dimension() + sage: bi = J.zero() + sage: bj = J.zero() + sage: bi_bj = J.zero()*J.zero() + sage: if n > 0: + ....: i = ZZ.random_element(n) + ....: j = ZZ.random_element(n) + ....: bi = J.monomial(i) + ....: bj = J.monomial(j) + ....: bi_bj = J.product_on_basis(i,j) + sage: bi*bj == bi_bj + True + """ - INPUT: + # We only stored the lower-triangular portion of the + # multiplication table. + if j <= i: + return self._multiplication_table[i][j] + else: + return self._multiplication_table[j][i] - * field -- the scalar field for this algebra (must be real) + def inner_product(self, x, y): + """ + The inner product associated with this Euclidean Jordan algebra. - * multiplication_table -- the multiplication table for this - algebra's implicit basis. Only the lower-triangular portion - of the table is used, since the multiplication is assumed - to be commutative. + Defaults to the trace inner product, but can be overridden by + subclasses if they are sure that the necessary properties are + satisfied. SETUP:: - sage: from mjo.eja.eja_algebra import ( - ....: FiniteDimensionalEuclideanJordanAlgebra, - ....: JordanSpinEJA, - ....: random_eja) + sage: from mjo.eja.eja_algebra import (random_eja, + ....: HadamardEJA, + ....: BilinearFormEJA) EXAMPLES: - By definition, Jordan multiplication commutes:: + Our inner product is "associative," which means the following for + a symmetric bilinear form:: - sage: set_random_seed() sage: J = random_eja() + sage: x,y,z = J.random_elements(3) + sage: (x*y).inner_product(z) == y.inner_product(x*z) + True + + TESTS: + + Ensure that this is the usual inner product for the algebras + over `R^n`:: + + sage: J = HadamardEJA.random_instance() sage: x,y = J.random_elements(2) - sage: x*y == y*x + sage: actual = x.inner_product(y) + sage: expected = x.to_vector().inner_product(y.to_vector()) + sage: actual == expected True - An error is raised if the Jordan product is not commutative:: + Ensure that this is one-half of the trace inner-product in a + BilinearFormEJA that isn't just the reals (when ``n`` isn't + one). This is in Faraut and Koranyi, and also my "On the + symmetry..." paper:: - sage: JP = ((1,2),(0,0)) - sage: IP = ((1,0),(0,1)) - sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP) - Traceback (most recent call last): - ... - ValueError: Jordan product is not commutative + sage: J = BilinearFormEJA.random_instance() + sage: n = J.dimension() + sage: x = J.random_element() + sage: y = J.random_element() + sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) + True + + """ + B = self._inner_product_matrix + return (B*x.to_vector()).inner_product(y.to_vector()) - An error is raised if the inner-product is not commutative:: - sage: JP = ((1,0),(0,1)) - sage: IP = ((1,2),(0,0)) - sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,JP,IP) - Traceback (most recent call last): - ... - ValueError: inner-product is not commutative + def is_associative(self): + r""" + Return whether or not this algebra's Jordan product is associative. + + SETUP:: + + sage: from mjo.eja.eja_algebra import ComplexHermitianEJA + + EXAMPLES:: + + sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False) + sage: J.is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: A.is_associative() + True + + """ + return "Associative" in self.category().axioms() + + def _is_commutative(self): + r""" + Whether or not this algebra's multiplication table is commutative. + + This method should of course always return ``True``, unless + this algebra was constructed with ``check_axioms=False`` and + passed an invalid multiplication table. + """ + return all( x*y == y*x for x in self.gens() for y in self.gens() ) + + def _is_jordanian(self): + r""" + Whether or not this algebra's multiplication table respects the + Jordan identity `(x^{2})(xy) = x(x^{2}y)`. + + We only check one arrangement of `x` and `y`, so for a + ``True`` result to be truly true, you should also check + :meth:`_is_commutative`. This method should of course always + return ``True``, unless this algebra was constructed with + ``check_axioms=False`` and passed an invalid multiplication table. + """ + return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j)) + == + (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j)) + for i in range(self.dimension()) + for j in range(self.dimension()) ) + + def _jordan_product_is_associative(self): + r""" + Return whether or not this algebra's Jordan product is + associative; that is, whether or not `x*(y*z) = (x*y)*z` + for all `x,y,x`. + + This method should agree with :meth:`is_associative` unless + you lied about the value of the ``associative`` parameter + when you constructed the algebra. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (random_eja, + ....: RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA) + + EXAMPLES:: + + sage: J = RealSymmetricEJA(4, orthonormalize=False) + sage: J._jordan_product_is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by() + sage: A._jordan_product_is_associative() + True + + :: + + sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: J._jordan_product_is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: A._jordan_product_is_associative() + True + + :: + + sage: J = QuaternionHermitianEJA(2) + sage: J._jordan_product_is_associative() + False + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by() + sage: A._jordan_product_is_associative() + True TESTS: - The ``field`` we're given must be real with ``check_field=True``:: + The values we've presupplied to the constructors agree with + the computation:: - sage: JordanSpinEJA(2,QQbar) - Traceback (most recent call last): - ... - ValueError: scalar field is not real + sage: J = random_eja() + sage: J.is_associative() == J._jordan_product_is_associative() + True - The multiplication table must be square with ``check_axioms=True``:: + """ + R = self.base_ring() - sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((),()),((1,),)) - Traceback (most recent call last): - ... - ValueError: multiplication table is not square + # Used to check whether or not something is zero. + epsilon = R.zero() + if not R.is_exact(): + # I don't know of any examples that make this magnitude + # necessary because I don't know how to make an + # associative algebra when the element subalgebra + # construction is unreliable (as it is over RDF; we can't + # find the degree of an element because we can't compute + # the rank of a matrix). But even multiplication of floats + # is non-associative, so *some* epsilon is needed... let's + # just take the one from _inner_product_is_associative? + epsilon = 1e-15 + + for i in range(self.dimension()): + for j in range(self.dimension()): + for k in range(self.dimension()): + x = self.monomial(i) + y = self.monomial(j) + z = self.monomial(k) + diff = (x*y)*z - x*(y*z) - The multiplication and inner-product tables must be the same - size (and in particular, the inner-product table must also be - square) with ``check_axioms=True``:: + if diff.norm() > epsilon: + return False - sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),(())) - Traceback (most recent call last): - ... - ValueError: multiplication and inner-product tables are - different sizes - sage: FiniteDimensionalEuclideanJordanAlgebra(QQ,((1,),),((1,2),)) - Traceback (most recent call last): - ... - ValueError: multiplication and inner-product tables are - different sizes + return True + def _inner_product_is_associative(self): + r""" + Return whether or not this algebra's inner product `B` is + associative; that is, whether or not `B(xy,z) = B(x,yz)`. + + This method should of course always return ``True``, unless + this algebra was constructed with ``check_axioms=False`` and + passed an invalid Jordan or inner-product. """ - if check_field: - if not field.is_subring(RR): - # Note: this does return true for the real algebraic - # field, the rationals, and any quadratic field where - # we've specified a real embedding. - raise ValueError("scalar field is not real") + R = self.base_ring() + # Used to check whether or not something is zero. + epsilon = R.zero() + if not R.is_exact(): + # This choice is sufficient to allow the construction of + # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True. + epsilon = 1e-15 - # The multiplication and inner-product tables should be square - # if the user wants us to verify them. And we verify them as - # soon as possible, because we want to exploit their symmetry. - n = len(multiplication_table) - if check_axioms: - if not all( len(l) == n for l in multiplication_table ): - raise ValueError("multiplication table is not square") - - # If the multiplication table is square, we can check if - # the inner-product table is square by comparing it to the - # multiplication table's dimensions. - msg = "multiplication and inner-product tables are different sizes" - if not len(inner_product_table) == n: - raise ValueError(msg) - - if not all( len(l) == n for l in inner_product_table ): - raise ValueError(msg) - - # Check commutativity of the Jordan product (symmetry of - # the multiplication table) and the commutativity of the - # inner-product (symmetry of the inner-product table) - # first if we're going to check them at all.. This has to - # be done before we define product_on_basis(), because - # that method assumes that self._multiplication_table is - # symmetric. And it has to be done before we build - # self._inner_product_matrix, because the process used to - # construct it assumes symmetry as well. - if not all( multiplication_table[j][i] - == multiplication_table[i][j] - for i in range(n) - for j in range(i+1) ): - raise ValueError("Jordan product is not commutative") - - if not all( inner_product_table[j][i] - == inner_product_table[i][j] - for i in range(n) - for j in range(i+1) ): - raise ValueError("inner-product is not commutative") - - self._matrix_basis = matrix_basis - - if category is None: - category = MagmaticAlgebras(field).FiniteDimensional() - category = category.WithBasis().Unital() - - fda = super(FiniteDimensionalEuclideanJordanAlgebra, self) - fda.__init__(field, - range(n), - prefix=prefix, - category=category) - self.print_options(bracket='') - - # The multiplication table we're given is necessarily in terms - # of vectors, because we don't have an algebra yet for - # anything to be an element of. However, it's faster in the - # long run to have the multiplication table be in terms of - # algebra elements. We do this after calling the superclass - # constructor so that from_vector() knows what to do. - # - # Note: we take advantage of symmetry here, and only store - # the lower-triangular portion of the table. - self._multiplication_table = [ [ self.vector_space().zero() - for j in range(i+1) ] - for i in range(n) ] + for i in range(self.dimension()): + for j in range(self.dimension()): + for k in range(self.dimension()): + x = self.monomial(i) + y = self.monomial(j) + z = self.monomial(k) + diff = (x*y).inner_product(z) - x.inner_product(y*z) - for i in range(n): - for j in range(i+1): - elt = self.from_vector(multiplication_table[i][j]) - self._multiplication_table[i][j] = elt - - self._multiplication_table = tuple(map(tuple, self._multiplication_table)) - - # Save our inner product as a matrix, since the efficiency of - # matrix multiplication will usually outweigh the fact that we - # have to store a redundant upper- or lower-triangular part. - # Pre-cache the fact that these are Hermitian (real symmetric, - # in fact) in case some e.g. matrix multiplication routine can - # take advantage of it. - ip_matrix_constructor = lambda i,j: inner_product_table[i][j] if j <= i else inner_product_table[j][i] - self._inner_product_matrix = matrix(field, n, ip_matrix_constructor) - self._inner_product_matrix._cache = {'hermitian': True} - self._inner_product_matrix.set_immutable() + if diff.abs() > epsilon: + return False - if check_axioms: - if not self._is_jordanian(): - raise ValueError("Jordan identity does not hold") - if not self._inner_product_is_associative(): - raise ValueError("inner product is not associative") + return True def _element_constructor_(self, elt): """ - Construct an element of this algebra from its vector or matrix - representation. + Construct an element of this algebra or a subalgebra from its + EJA element, vector, or matrix representation. This gets called only after the parent element _call_ method fails to find a coercion for the argument. SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + sage: from mjo.eja.eja_algebra import (random_eja, + ....: JordanSpinEJA, ....: HadamardEJA, ....: RealSymmetricEJA) @@ -276,34 +728,69 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): ... ValueError: not an element of this algebra + Tuples work as well, provided that the matrix basis for the + algebra consists of them:: + + sage: J1 = HadamardEJA(3) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) ) + b1 + b5 + + Subalgebra elements are embedded into the superalgebra:: + + sage: J = JordanSpinEJA(3) + sage: J.one() + b0 + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by() + sage: J(A.one()) + b0 + TESTS: - Ensure that we can convert any element of the two non-matrix - simple algebras (whose matrix representations are columns) - back and forth faithfully:: + Ensure that we can convert any element back and forth + faithfully between its matrix and algebra representations:: - sage: set_random_seed() - sage: J = HadamardEJA.random_instance() - sage: x = J.random_element() - sage: J(x.to_vector().column()) == x - True - sage: J = JordanSpinEJA.random_instance() + sage: J = random_eja() sage: x = J.random_element() - sage: J(x.to_vector().column()) == x + sage: J(x.to_matrix()) == x True + + We cannot coerce elements between algebras just because their + matrix representations are compatible:: + + sage: J1 = HadamardEJA(3) + sage: J2 = JordanSpinEJA(3) + sage: J2(J1.one()) + Traceback (most recent call last): + ... + ValueError: not an element of this algebra + sage: J1(J2.zero()) + Traceback (most recent call last): + ... + ValueError: not an element of this algebra + """ msg = "not an element of this algebra" - if elt == 0: - # The superclass implementation of random_element() - # needs to be able to coerce "0" into the algebra. - return self.zero() - elif elt in self.base_ring(): + if elt in self.base_ring(): # Ensure that no base ring -> algebra coercion is performed # by this method. There's some stupidity in sage that would # otherwise propagate to this method; for example, sage thinks # that the integer 3 belongs to the space of 2-by-2 matrices. raise ValueError(msg) + if hasattr(elt, 'superalgebra_element'): + # Handle subalgebra elements + if elt.parent().superalgebra() == self: + return elt.superalgebra_element() + + if hasattr(elt, 'sparse_vector'): + # Convert a vector into a column-matrix. We check for + # "sparse_vector" and not "column" because matrices also + # have a "column" method. + elt = elt.column() + if elt not in self.matrix_space(): raise ValueError(msg) @@ -313,105 +800,42 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # element's ring because the basis space might be an algebraic # closure whereas the base ring of the 3-by-3 identity matrix # could be QQ instead of QQbar. - V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols()) - W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() ) + # + # And, we also have to handle Cartesian product bases (when + # the matrix basis consists of tuples) here. The "good news" + # is that we're already converting everything to long vectors, + # and that strategy works for tuples as well. + # + elt = self._matrix_span.ambient_vector_space()(_all2list(elt)) try: - coords = W.coordinate_vector(_mat2vec(elt)) + coords = self._matrix_span.coordinate_vector(elt) except ArithmeticError: # vector is not in free module raise ValueError(msg) - return self.from_vector(coords) - - def _repr_(self): - """ - Return a string representation of ``self``. - - SETUP:: - - sage: from mjo.eja.eja_algebra import JordanSpinEJA - - TESTS: - - Ensure that it says what we think it says:: - - sage: JordanSpinEJA(2, field=AA) - Euclidean Jordan algebra of dimension 2 over Algebraic Real Field - sage: JordanSpinEJA(3, field=RDF) - Euclidean Jordan algebra of dimension 3 over Real Double Field - - """ - fmt = "Euclidean Jordan algebra of dimension {} over {}" - return fmt.format(self.dimension(), self.base_ring()) - - def product_on_basis(self, i, j): - # We only stored the lower-triangular portion of the - # multiplication table. - if j <= i: - return self._multiplication_table[i][j] - else: - return self._multiplication_table[j][i] - - def _is_commutative(self): - r""" - Whether or not this algebra's multiplication table is commutative. - - This method should of course always return ``True``, unless - this algebra was constructed with ``check_axioms=False`` and - passed an invalid multiplication table. - """ - return all( self.product_on_basis(i,j) == self.product_on_basis(i,j) - for i in range(self.dimension()) - for j in range(self.dimension()) ) - - def _is_jordanian(self): - r""" - Whether or not this algebra's multiplication table respects the - Jordan identity `(x^{2})(xy) = x(x^{2}y)`. + return self.from_vector(coords) - We only check one arrangement of `x` and `y`, so for a - ``True`` result to be truly true, you should also check - :meth:`_is_commutative`. This method should of course always - return ``True``, unless this algebra was constructed with - ``check_axioms=False`` and passed an invalid multiplication table. + def _repr_(self): """ - return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j)) - == - (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j)) - for i in range(self.dimension()) - for j in range(self.dimension()) ) + Return a string representation of ``self``. - def _inner_product_is_associative(self): - r""" - Return whether or not this algebra's inner product `B` is - associative; that is, whether or not `B(xy,z) = B(x,yz)`. + SETUP:: - This method should of course always return ``True``, unless - this algebra was constructed with ``check_axioms=False`` and - passed an invalid multiplication table. - """ + sage: from mjo.eja.eja_algebra import JordanSpinEJA - # Used to check whether or not something is zero in an inexact - # ring. This number is sufficient to allow the construction of - # QuaternionHermitianEJA(2, RDF) with check_axioms=True. - epsilon = 1e-16 + TESTS: - for i in range(self.dimension()): - for j in range(self.dimension()): - for k in range(self.dimension()): - x = self.monomial(i) - y = self.monomial(j) - z = self.monomial(k) - diff = (x*y).inner_product(z) - x.inner_product(y*z) + Ensure that it says what we think it says:: + + sage: JordanSpinEJA(2, field=AA) + Euclidean Jordan algebra of dimension 2 over Algebraic Real Field + sage: JordanSpinEJA(3, field=RDF) + Euclidean Jordan algebra of dimension 3 over Real Double Field - if self.base_ring().is_exact(): - if diff != 0: - return False - else: - if diff.abs() > epsilon: - return False + """ + fmt = "Euclidean Jordan algebra of dimension {} over {}" + return fmt.format(self.dimension(), self.base_ring()) - return True @cached_method def characteristic_polynomial_of(self): @@ -440,7 +864,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = JordanSpinEJA(3) sage: p = J.characteristic_polynomial_of(); p - X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2 + X0^2 - X1^2 - X2^2 + (-2*t)*X0 + t^2 sage: xvec = J.one().to_vector() sage: p(*xvec) t^2 - 2*t + 1 @@ -489,13 +913,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = HadamardEJA(2) sage: J.coordinate_polynomial_ring() - Multivariate Polynomial Ring in X1, X2... - sage: J = RealSymmetricEJA(3,QQ,orthonormalize=False) + Multivariate Polynomial Ring in X0, X1... + sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False) sage: J.coordinate_polynomial_ring() - Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6... + Multivariate Polynomial Ring in X0, X1, X2, X3, X4, X5... """ - var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) ) + var_names = tuple( "X%d" % z for z in range(self.dimension()) ) return PolynomialRing(self.base_ring(), var_names) def inner_product(self, x, y): @@ -517,7 +941,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Our inner product is "associative," which means the following for a symmetric bilinear form:: - sage: set_random_seed() sage: J = random_eja() sage: x,y,z = J.random_elements(3) sage: (x*y).inner_product(z) == y.inner_product(x*z) @@ -528,7 +951,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Ensure that this is the usual inner product for the algebras over `R^n`:: - sage: set_random_seed() sage: J = HadamardEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = x.inner_product(y) @@ -541,7 +963,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): one). This is in Faraut and Koranyi, and also my "On the symmetry..." paper:: - sage: set_random_seed() sage: J = BilinearFormEJA.random_instance() sage: n = J.dimension() sage: x = J.random_element() @@ -594,33 +1015,28 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = JordanSpinEJA(4) sage: J.multiplication_table() +----++----+----+----+----+ - | * || e0 | e1 | e2 | e3 | + | * || b0 | b1 | b2 | b3 | +====++====+====+====+====+ - | e0 || e0 | e1 | e2 | e3 | + | b0 || b0 | b1 | b2 | b3 | +----++----+----+----+----+ - | e1 || e1 | e0 | 0 | 0 | + | b1 || b1 | b0 | 0 | 0 | +----++----+----+----+----+ - | e2 || e2 | 0 | e0 | 0 | + | b2 || b2 | 0 | b0 | 0 | +----++----+----+----+----+ - | e3 || e3 | 0 | 0 | e0 | + | b3 || b3 | 0 | 0 | b0 | +----++----+----+----+----+ """ n = self.dimension() - M = [ [ self.zero() for j in range(n) ] - for i in range(n) ] - for i in range(n): - for j in range(i+1): - M[i][j] = self._multiplication_table[i][j] - M[j][i] = M[i][j] + # Prepend the header row. + M = [["*"] + list(self.gens())] - for i in range(n): - # Prepend the left "header" column entry Can't do this in - # the loop because it messes up the symmetry. - M[i] = [self.monomial(i)] + M[i] + # And to each subsequent row, prepend an entry that belongs to + # the left-side "header column." + M += [ [self.monomial(i)] + [ self.monomial(i)*self.monomial(j) + for j in range(n) ] + for i in range(n) ] - # Prepend the header row. - M = [["*"] + list(self.gens())] + M return table(M, header_row=True, header_column=True, frame=True) @@ -648,7 +1064,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Why implement this for non-matrix algebras? Avoiding special cases for the :class:`BilinearFormEJA` pays with simplicity in its own right. But mainly, we would like to be able to assume - that elements of a :class:`DirectSumEJA` can be displayed + that elements of a :class:`CartesianProductEJA` can be displayed nicely, without having to have special classes for direct sums one of whose components was a matrix algebra. @@ -661,7 +1077,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = RealSymmetricEJA(2) sage: J.basis() - Finite family {0: e0, 1: e1, 2: e2} + Finite family {0: b0, 1: b1, 2: b2} sage: J.matrix_basis() ( [1 0] [ 0 0.7071067811865475?] [0 0] @@ -672,18 +1088,14 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: J = JordanSpinEJA(2) sage: J.basis() - Finite family {0: e0, 1: e1} + Finite family {0: b0, 1: b1} sage: J.matrix_basis() ( [1] [0] [0], [1] ) """ - if self._matrix_basis is None: - M = self.matrix_space() - return tuple( M(b.to_vector()) for b in self.basis() ) - else: - return self._matrix_basis + return self._matrix_basis def matrix_space(self): @@ -692,19 +1104,54 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): we think of them as matrices (including column vectors of the appropriate size). - Generally this will be an `n`-by-`1` column-vector space, + "By default" this will be an `n`-by-`1` column-matrix space, except when the algebra is trivial. There it's `n`-by-`n` (where `n` is zero), to ensure that two elements of the matrix - space (empty matrices) can be multiplied. + space (empty matrices) can be multiplied. For algebras of + matrices, this returns the space in which their + real embeddings live. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, + ....: JordanSpinEJA, + ....: QuaternionHermitianEJA, + ....: TrivialEJA) + + EXAMPLES: + + By default, the matrix representation is just a column-matrix + equivalent to the vector representation:: + + sage: J = JordanSpinEJA(3) + sage: J.matrix_space() + Full MatrixSpace of 3 by 1 dense matrices over Algebraic + Real Field + + The matrix representation in the trivial algebra is + zero-by-zero instead of the usual `n`-by-one:: + + sage: J = TrivialEJA() + sage: J.matrix_space() + Full MatrixSpace of 0 by 0 dense matrices over Algebraic + Real Field + + The matrix space for complex/quaternion Hermitian matrix EJA + is the space in which their real-embeddings live, not the + original complex/quaternion matrix space:: + + sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: J.matrix_space() + Module of 2 by 2 matrices with entries in Algebraic Field over + the scalar ring Rational Field + sage: J = QuaternionHermitianEJA(1,field=QQ,orthonormalize=False) + sage: J.matrix_space() + Module of 1 by 1 matrices with entries in Quaternion + Algebra (-1, -1) with base ring Rational Field over + the scalar ring Rational Field - Matrix algebras override this with something more useful. """ - if self.is_trivial(): - return MatrixSpace(self.base_ring(), 0) - elif self._matrix_basis is None or len(self._matrix_basis) == 0: - return MatrixSpace(self.base_ring(), self.dimension(), 1) - else: - return self._matrix_basis[0].matrix_space() + return self._matrix_space @cached_method @@ -717,41 +1164,99 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: from mjo.eja.eja_algebra import (HadamardEJA, ....: random_eja) - EXAMPLES:: + EXAMPLES: + + We can compute unit element in the Hadamard EJA:: sage: J = HadamardEJA(5) sage: J.one() - e0 + e1 + e2 + e3 + e4 + b0 + b1 + b2 + b3 + b4 + + The unit element in the Hadamard EJA is inherited in the + subalgebras generated by its elements:: + + sage: J = HadamardEJA(5) + sage: J.one() + b0 + b1 + b2 + b3 + b4 + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: A.one() + c0 + sage: A.one().superalgebra_element() + b0 + b1 + b2 + b3 + b4 TESTS: - The identity element acts like the identity:: + The identity element acts like the identity, regardless of + whether or not we orthonormalize:: - sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: J.one()*x == x and x*J.one() == x True + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: y = A.random_element() + sage: A.one()*y == y and y*A.one() == y + True + + :: + + sage: J = random_eja(field=QQ, orthonormalize=False) + sage: x = J.random_element() + sage: J.one()*x == x and x*J.one() == x + True + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: y = A.random_element() + sage: A.one()*y == y and y*A.one() == y + True - The matrix of the unit element's operator is the identity:: + The matrix of the unit element's operator is the identity, + regardless of the base field and whether or not we + orthonormalize:: - sage: set_random_seed() sage: J = random_eja() sage: actual = J.one().operator().matrix() sage: expected = matrix.identity(J.base_ring(), J.dimension()) sage: actual == expected True + sage: x = J.random_element() + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: actual = A.one().operator().matrix() + sage: expected = matrix.identity(A.base_ring(), A.dimension()) + sage: actual == expected + True + + :: + + sage: J = random_eja(field=QQ, orthonormalize=False) + sage: actual = J.one().operator().matrix() + sage: expected = matrix.identity(J.base_ring(), J.dimension()) + sage: actual == expected + True + sage: x = J.random_element() + sage: A = x.subalgebra_generated_by(orthonormalize=False) + sage: actual = A.one().operator().matrix() + sage: expected = matrix.identity(A.base_ring(), A.dimension()) + sage: actual == expected + True Ensure that the cached unit element (often precomputed by hand) agrees with the computed one:: - sage: set_random_seed() sage: J = random_eja() sage: cached = J.one() sage: J.one.clear_cache() sage: J.one() == cached True + :: + + sage: J = random_eja(field=QQ, orthonormalize=False) + sage: cached = J.one() + sage: J.one.clear_cache() + sage: J.one() == cached + True + """ # We can brute-force compute the matrices of the operators # that correspond to the basis elements of this algebra. @@ -761,7 +1266,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # # Of course, matrices aren't vectors in sage, so we have to # appeal to the "long vectors" isometry. - oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ] + + V = VectorSpace(self.base_ring(), self.dimension()**2) + oper_vecs = [ V(g.operator().matrix().list()) for g in self.gens() ] # Now we use basic linear algebra to find the coefficients, # of the matrices-as-vectors-linear-combination, which should @@ -771,7 +1278,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # We used the isometry on the left-hand side already, but we # still need to do it for the right-hand side. Recall that we # wanted something that summed to the identity matrix. - b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) ) + b = V( matrix.identity(self.base_ring(), self.dimension()).list() ) # Now if there's an identity element in the algebra, this # should work. We solve on the left to avoid having to @@ -856,7 +1363,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Every algebra decomposes trivially with respect to its identity element:: - sage: set_random_seed() sage: J = random_eja() sage: J0,J5,J1 = J.peirce_decomposition(J.one()) sage: J0.dimension() == 0 and J5.dimension() == 0 @@ -869,7 +1375,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): elements in the two subalgebras are the projections onto their respective subspaces of the superalgebra's identity element:: - sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: if not J.is_trivial(): @@ -896,14 +1401,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): if not c.is_idempotent(): raise ValueError("element is not idempotent: %s" % c) - from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra - # Default these to what they should be if they turn out to be # trivial, because eigenspaces_left() won't return eigenvalues # corresponding to trivial spaces (e.g. it returns only the # eigenspace corresponding to lambda=1 if you take the # decomposition relative to the identity element). - trivial = FiniteDimensionalEuclideanJordanSubalgebra(self, ()) + trivial = self.subalgebra((), check_axioms=False) J0 = trivial # eigenvalue zero J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half J1 = trivial # eigenvalue one @@ -913,9 +1416,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): J5 = eigspace else: gens = tuple( self.from_vector(b) for b in eigspace.basis() ) - subalg = FiniteDimensionalEuclideanJordanSubalgebra(self, - gens, - check_axioms=False) + subalg = self.subalgebra(gens, check_axioms=False) if eigval == 0: J0 = subalg elif eigval == 1: @@ -947,26 +1448,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # For a general base ring... maybe we can trust this to do the # right thing? Unlikely, but. V = self.vector_space() - v = V.random_element() - - if self.base_ring() is AA: - # The "random element" method of the algebraic reals is - # stupid at the moment, and only returns integers between - # -2 and 2, inclusive: - # - # https://trac.sagemath.org/ticket/30875 - # - # Instead, we implement our own "random vector" method, - # and then coerce that into the algebra. We use the vector - # space degree here instead of the dimension because a - # subalgebra could (for example) be spanned by only two - # vectors, each with five coordinates. We need to - # generate all five coordinates. - if thorough: - v *= QQbar.random_element().real() - else: - v *= QQ.random_element() + if self.base_ring() is AA and not thorough: + # Now that AA generates actually random random elements + # (post Trac 30875), we only need to de-thorough the + # randomness when asked to. + V = V.change_ring(QQ) + v = V.random_element() return self.from_vector(V.coordinate_vector(v)) def random_elements(self, count, thorough=False): @@ -999,24 +1487,89 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): for idx in range(count) ) + def operator_polynomial_matrix(self): + r""" + Return the matrix of polynomials (over this algebra's + :meth:`coordinate_polynomial_ring`) that, when evaluated at + the basis coordinates of an element `x`, produces the basis + representation of `L_{x}`. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: JordanSpinEJA) + + EXAMPLES:: + + sage: J = HadamardEJA(4) + sage: L_x = J.operator_polynomial_matrix() + sage: L_x + [X0 0 0 0] + [ 0 X1 0 0] + [ 0 0 X2 0] + [ 0 0 0 X3] + sage: x = J.one() + sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector()) + sage: L_x.subs(dict(d)) + [1 0 0 0] + [0 1 0 0] + [0 0 1 0] + [0 0 0 1] + + :: + + sage: J = JordanSpinEJA(4) + sage: L_x = J.operator_polynomial_matrix() + sage: L_x + [X0 X1 X2 X3] + [X1 X0 0 0] + [X2 0 X0 0] + [X3 0 0 X0] + sage: x = J.one() + sage: d = zip(J.coordinate_polynomial_ring().gens(), x.to_vector()) + sage: L_x.subs(dict(d)) + [1 0 0 0] + [0 1 0 0] + [0 0 1 0] + [0 0 0 1] + + """ + R = self.coordinate_polynomial_ring() + + def L_x_i_j(i,j): + # From a result in my book, these are the entries of the + # basis representation of L_x. + return sum( v*self.monomial(k).operator().matrix()[i,j] + for (k,v) in enumerate(R.gens()) ) + + n = self.dimension() + return matrix(R, n, n, L_x_i_j) + @cached_method def _charpoly_coefficients(self): r""" The `r` polynomial coefficients of the "characteristic polynomial of" function. + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS: + + The theory shows that these are all homogeneous polynomials of + a known degree:: + + sage: J = random_eja() + sage: all(p.is_homogeneous() for p in J._charpoly_coefficients()) + True + """ n = self.dimension() R = self.coordinate_polynomial_ring() - vars = R.gens() F = R.fraction_field() - def L_x_i_j(i,j): - # From a result in my book, these are the entries of the - # basis representation of L_x. - return sum( vars[k]*self.monomial(k).operator().matrix()[i,j] - for k in range(n) ) - - L_x = matrix(F, n, n, L_x_i_j) + L_x = self.operator_polynomial_matrix() r = None if self.rank.is_in_cache(): @@ -1039,10 +1592,17 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # The theory says that only the first "r" coefficients are # nonzero, and they actually live in the original polynomial - # ring and not the fraction field. We negate them because - # in the actual characteristic polynomial, they get moved - # to the other side where x^r lives. - return -A_rref.solve_right(E*b).change_ring(R)[:r] + # ring and not the fraction field. We negate them because in + # the actual characteristic polynomial, they get moved to the + # other side where x^r lives. We don't bother to trim A_rref + # down to a square matrix and solve the resulting system, + # because the upper-left r-by-r portion of A_rref is + # guaranteed to be the identity matrix, so e.g. + # + # A_rref.solve_right(Y) + # + # would just be returning Y. + return (-E*b)[:r].change_ring(R) @cached_method def rank(self): @@ -1090,7 +1650,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): positive integer rank, unless the algebra is trivial in which case its rank will be zero:: - sage: set_random_seed() sage: J = random_eja() sage: r = J.rank() sage: r in ZZ @@ -1101,9 +1660,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Ensure that computing the rank actually works, since the ranks of all simple algebras are known and will be cached by default:: - sage: set_random_seed() # long time sage: J = random_eja() # long time - sage: caches = J.rank() # long time + sage: cached = J.rank() # long time sage: J.rank.clear_cache() # long time sage: J.rank() == cached # long time True @@ -1112,6 +1670,14 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return len(self._charpoly_coefficients()) + def subalgebra(self, basis, **kwargs): + r""" + Create a subalgebra of this algebra from the given basis. + """ + from mjo.eja.eja_subalgebra import EJASubalgebra + return EJASubalgebra(self, basis, **kwargs) + + def vector_space(self): """ Return the vector space that underlies this algebra. @@ -1130,11 +1696,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return self.zero().to_vector().parent().ambient_vector_space() - Element = FiniteDimensionalEuclideanJordanAlgebraElement -class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): +class RationalBasisEJA(EJA): r""" - New class for algebras whose supplied basis elements have all rational entries. + Algebras whose supplied basis elements have all rational entries. SETUP:: @@ -1155,181 +1720,57 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr """ def __init__(self, - field, basis, jordan_product, inner_product, - orthonormalize=True, - prefix='e', - category=None, + field=AA, check_field=True, - check_axioms=True): + **kwargs): if check_field: # Abuse the check_field parameter to check that the entries of # out basis (in ambient coordinates) are in the field QQ. - if not all( all(b_i in QQ for b_i in b.list()) for b in basis ): + # Use _all2list to get the vector coordinates of octonion + # entries and not the octonions themselves (which are not + # rational). + if not all( all(b_i in QQ for b_i in _all2list(b)) + for b in basis ): raise TypeError("basis not rational") - n = len(basis) - vector_basis = basis - - from sage.structure.element import is_Matrix - basis_is_matrices = False - - degree = 0 - if n > 0: - if is_Matrix(basis[0]): - basis_is_matrices = True - from mjo.eja.eja_utils import _vec2mat - vector_basis = tuple( map(_mat2vec,basis) ) - degree = basis[0].nrows()**2 - else: - degree = basis[0].degree() - - V = VectorSpace(field, degree) + super().__init__(basis, + jordan_product, + inner_product, + field=field, + check_field=check_field, + **kwargs) - # If we were asked to orthonormalize, and if the orthonormal - # basis is different from the given one, then we also want to - # compute multiplication and inner-product tables for the - # deorthonormalized basis. These can be used later to - # construct a deorthonormalized copy of this algebra over QQ - # in which several operations are much faster. self._rational_algebra = None - - if orthonormalize: - if self.base_ring() is not QQ: - # There's no point in constructing the extra algebra if this - # one is already rational. If the original basis is rational - # but normalization would make it irrational, then this whole - # constructor will just fail anyway as it tries to stick an - # irrational number into a rational algebra. - # - # Note: the same Jordan and inner-products work here, - # because they are necessarily defined with respect to - # ambient coordinates and not any particular basis. - self._rational_algebra = RationalBasisEuclideanJordanAlgebra( - QQ, - basis, - jordan_product, - inner_product, - orthonormalize=False, - prefix=prefix, - category=category, - check_field=False, - check_axioms=False) - - # Compute the deorthonormalized tables before we orthonormalize - # the given basis. - W = V.span_of_basis( vector_basis ) - - # Note: the Jordan and inner-products are defined in terms - # of the ambient basis. It's important that their arguments - # are in ambient coordinates as well. - for i in range(n): - for j in range(i+1): - # given basis w.r.t. ambient coords - q_i = vector_basis[i] - q_j = vector_basis[j] - - if basis_is_matrices: - q_i = _vec2mat(q_i) - q_j = _vec2mat(q_j) - - elt = jordan_product(q_i, q_j) - ip = inner_product(q_i, q_j) - - if basis_is_matrices: - # do another mat2vec because the multiplication - # table is in terms of vectors - elt = _mat2vec(elt) - - # We overwrite the name "vector_basis" in a second, but never modify it - # in place, to this effectively makes a copy of it. - deortho_vector_basis = vector_basis - self._deortho_matrix = None - - if orthonormalize: - from mjo.eja.eja_utils import gram_schmidt - if basis_is_matrices: - vector_ip = lambda x,y: inner_product(_vec2mat(x), _vec2mat(y)) - vector_basis = gram_schmidt(vector_basis, vector_ip) - else: - vector_basis = gram_schmidt(vector_basis, inner_product) - - W = V.span_of_basis( vector_basis ) - - # Normalize the "matrix" basis, too! - basis = vector_basis - - if basis_is_matrices: - basis = tuple( map(_vec2mat,basis) ) - - W = V.span_of_basis( vector_basis ) - - # Now "W" is the vector space of our algebra coordinates. The - # variables "X1", "X2",... refer to the entries of vectors in - # W. Thus to convert back and forth between the orthonormal - # coordinates and the given ones, we need to stick the original - # basis in W. - U = V.span_of_basis( deortho_vector_basis ) - self._deortho_matrix = matrix( U.coordinate_vector(q) - for q in vector_basis ) - - # If the superclass constructor is going to verify the - # symmetry of this table, it has better at least be - # square... - if check_axioms: - mult_table = [ [0 for j in range(n)] for i in range(n) ] - ip_table = [ [0 for j in range(n)] for i in range(n) ] - else: - mult_table = [ [0 for j in range(i+1)] for i in range(n) ] - ip_table = [ [0 for j in range(i+1)] for i in range(n) ] - - # Note: the Jordan and inner-products are defined in terms - # of the ambient basis. It's important that their arguments - # are in ambient coordinates as well. - for i in range(n): - for j in range(i+1): - # ortho basis w.r.t. ambient coords - q_i = vector_basis[i] - q_j = vector_basis[j] - - if basis_is_matrices: - q_i = _vec2mat(q_i) - q_j = _vec2mat(q_j) - - elt = jordan_product(q_i, q_j) - ip = inner_product(q_i, q_j) - - if basis_is_matrices: - # do another mat2vec because the multiplication - # table is in terms of vectors - elt = _mat2vec(elt) - - elt = W.coordinate_vector(elt) - mult_table[i][j] = elt - ip_table[i][j] = ip - if check_axioms: - # The tables are square if we're verifying that they - # are commutative. - mult_table[j][i] = elt - ip_table[j][i] = ip - - if basis_is_matrices: - for m in basis: - m.set_immutable() + if field is not QQ: + # There's no point in constructing the extra algebra if this + # one is already rational. + # + # Note: the same Jordan and inner-products work here, + # because they are necessarily defined with respect to + # ambient coordinates and not any particular basis. + self._rational_algebra = EJA( + basis, + jordan_product, + inner_product, + field=QQ, + matrix_space=self.matrix_space(), + associative=self.is_associative(), + orthonormalize=False, + check_field=False, + check_axioms=False) + + def rational_algebra(self): + # Using None as a flag here (rather than just assigning "self" + # to self._rational_algebra by default) feels a little bit + # more sane to me in a garbage-collected environment. + if self._rational_algebra is None: + return self else: - basis = tuple( x.column() for x in basis ) - - super().__init__(field, - mult_table, - ip_table, - prefix, - category, - basis, # matrix basis - check_field, - check_axioms) + return self._rational_algebra @cached_method def _charpoly_coefficients(self): @@ -1347,7 +1788,7 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr sage: J = JordanSpinEJA(3) sage: J._charpoly_coefficients() - (X1^2 - X2^2 - X3^2, -2*X1) + (X0^2 - X1^2 - X2^2, -2*X0) sage: a0 = J._charpoly_coefficients()[0] sage: J.base_ring() Algebraic Real Field @@ -1355,20 +1796,17 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr Algebraic Real Field """ - if self.base_ring() is QQ: - # There's no need to construct *another* algebra over the - # rationals if this one is already over the rationals. - superclass = super(RationalBasisEuclideanJordanAlgebra, self) - return superclass._charpoly_coefficients() - - # Do the computation over the rationals. The answer will be - # the same, because all we've done is a change of basis. - # Then, change back from QQ to our real base ring + if self.rational_algebra() is self: + # Bypass the hijinks if they won't benefit us. + return super()._charpoly_coefficients() + + # Do the computation over the rationals. a = ( a_i.change_ring(self.base_ring()) - for a_i in self._rational_algebra._charpoly_coefficients() ) + for a_i in self.rational_algebra()._charpoly_coefficients() ) - # Now convert the coordinate variables back to the - # deorthonormalized ones. + # Convert our coordinate variables into deorthonormalized ones + # and substitute them into the deorthonormalized charpoly + # coefficients. R = self.coordinate_polynomial_ring() from sage.modules.free_module_element import vector X = vector(R, R.gens()) @@ -1377,7 +1815,7 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr subs_dict = { X[i]: BX[i] for i in range(len(X)) } return tuple( a_i.subs(subs_dict) for a_i in a ) -class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra): +class ConcreteEJA(EJA): r""" A class for the Euclidean Jordan algebras that we know by name. @@ -1388,15 +1826,14 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra): SETUP:: - sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra + sage: from mjo.eja.eja_algebra import ConcreteEJA TESTS: Our basis is normalized with respect to the algebra's inner product, unless we specify otherwise:: - sage: set_random_seed() - sage: J = ConcreteEuclideanJordanAlgebra.random_instance() + sage: J = ConcreteEJA.random_instance() sage: all( b.norm() == 1 for b in J.gens() ) True @@ -1406,110 +1843,232 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra): natural->EJA basis representation is an isometry and within the EJA the operator is self-adjoint by the Jordan axiom:: - sage: set_random_seed() - sage: J = ConcreteEuclideanJordanAlgebra.random_instance() + sage: J = ConcreteEJA.random_instance() sage: x = J.random_element() sage: x.operator().is_self_adjoint() True """ @staticmethod - def _max_random_instance_size(): + def _max_random_instance_dimension(): + r""" + The maximum dimension of any random instance. Ten dimensions seems + to be about the point where everything takes a turn for the + worse. And dimension ten (but not nine) allows the 4-by-4 real + Hermitian matrices, the 2-by-2 quaternion Hermitian matrices, + and the 2-by-2 octonion Hermitian matrices. + """ + return 10 + + @staticmethod + def _max_random_instance_size(max_dimension): """ Return an integer "size" that is an upper bound on the size of - this algebra when it is used in a random test - case. Unfortunately, the term "size" is ambiguous -- when - dealing with `R^n` under either the Hadamard or Jordan spin - product, the "size" refers to the dimension `n`. When dealing - with a matrix algebra (real symmetric or complex/quaternion - Hermitian), it refers to the size of the matrix, which is far - less than the dimension of the underlying vector space. + this algebra when it is used in a random test case. This size + (which can be passed to the algebra's constructor) is itself + based on the ``max_dimension`` parameter. + + This method must be implemented in each subclass. + """ + raise NotImplementedError + + @classmethod + def random_instance(cls, max_dimension=None, *args, **kwargs): + """ + Return a random instance of this type of algebra whose dimension + is less than or equal to the lesser of ``max_dimension`` and + the value returned by ``_max_random_instance_dimension()``. If + the dimension bound is omitted, then only the + ``_max_random_instance_dimension()`` is used as a bound. + + This method should be implemented in each subclass. + + SETUP:: + + sage: from mjo.eja.eja_algebra import ConcreteEJA + + TESTS: + + Both the class bound and the ``max_dimension`` argument are upper + bounds on the dimension of the algebra returned:: + + sage: from sage.misc.prandom import choice + sage: eja_class = choice(ConcreteEJA.__subclasses__()) + sage: class_max_d = eja_class._max_random_instance_dimension() + sage: J = eja_class.random_instance(max_dimension=20, + ....: field=QQ, + ....: orthonormalize=False) + sage: J.dimension() <= class_max_d + True + sage: J = eja_class.random_instance(max_dimension=2, + ....: field=QQ, + ....: orthonormalize=False) + sage: J.dimension() <= 2 + True + + """ + from sage.misc.prandom import choice + eja_class = choice(cls.__subclasses__()) + + # These all bubble up to the RationalBasisEJA superclass + # constructor, so any (kw)args valid there are also valid + # here. + return eja_class.random_instance(max_dimension, *args, **kwargs) + + +class HermitianMatrixEJA(EJA): + @staticmethod + def _denormalized_basis(A): + """ + Returns a basis for the given Hermitian matrix space. + + Why do we embed these? Basically, because all of numerical linear + algebra assumes that you're working with vectors consisting of `n` + entries from a field and scalars from the same field. There's no way + to tell SageMath that (for example) the vectors contain complex + numbers, while the scalar field is real. + + SETUP:: + + sage: from mjo.hurwitz import (ComplexMatrixAlgebra, + ....: QuaternionMatrixAlgebra, + ....: OctonionMatrixAlgebra) + sage: from mjo.eja.eja_algebra import HermitianMatrixEJA + + TESTS:: + + sage: n = ZZ.random_element(1,5) + sage: A = MatrixSpace(QQ, n) + sage: B = HermitianMatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B) + True + + :: + + sage: n = ZZ.random_element(1,5) + sage: A = ComplexMatrixAlgebra(n, scalars=QQ) + sage: B = HermitianMatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B) + True + + :: + + sage: n = ZZ.random_element(1,5) + sage: A = QuaternionMatrixAlgebra(n, scalars=QQ) + sage: B = HermitianMatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B ) + True + + :: - This method must be implemented in each subclass. - """ - raise NotImplementedError + sage: n = ZZ.random_element(1,5) + sage: A = OctonionMatrixAlgebra(n, scalars=QQ) + sage: B = HermitianMatrixEJA._denormalized_basis(A) + sage: all( M.is_hermitian() for M in B ) + True - @classmethod - def random_instance(cls, *args, **kwargs): """ - Return a random instance of this type of algebra. + # These work for real MatrixSpace, whose monomials only have + # two coordinates (because the last one would always be "1"). + es = A.base_ring().gens() + gen = lambda A,m: A.monomial(m[:2]) - This method should be implemented in each subclass. - """ - from sage.misc.prandom import choice - eja_class = choice(cls.__subclasses__()) + if hasattr(A, 'entry_algebra_gens'): + # We've got a MatrixAlgebra, and its monomials will have + # three coordinates. + es = A.entry_algebra_gens() + gen = lambda A,m: A.monomial(m) - # These all bubble up to the RationalBasisEuclideanJordanAlgebra - # superclass constructor, so any (kw)args valid there are also - # valid here. - return eja_class.random_instance(*args, **kwargs) + basis = [] + for i in range(A.nrows()): + for j in range(i+1): + if i == j: + E_ii = gen(A, (i,j,es[0])) + basis.append(E_ii) + else: + for e in es: + E_ij = gen(A, (i,j,e)) + E_ij += E_ij.conjugate_transpose() + basis.append(E_ij) + return tuple( basis ) -class MatrixEuclideanJordanAlgebra: @staticmethod - def real_embed(M): - """ - Embed the matrix ``M`` into a space of real matrices. + def jordan_product(X,Y): + return (X*Y + Y*X)/2 - The matrix ``M`` can have entries in any field at the moment: - the real numbers, complex numbers, or quaternions. And although - they are not a field, we can probably support octonions at some - point, too. This function returns a real matrix that "acts like" - the original with respect to matrix multiplication; i.e. + @staticmethod + def trace_inner_product(X,Y): + r""" + A trace inner-product for matrices that aren't embedded in the + reals. It takes MATRICES as arguments, not EJA elements. - real_embed(M*N) = real_embed(M)*real_embed(N) + SETUP:: - """ - raise NotImplementedError + sage: from mjo.eja.eja_algebra import (RealSymmetricEJA, + ....: ComplexHermitianEJA, + ....: QuaternionHermitianEJA, + ....: OctonionHermitianEJA) + EXAMPLES:: - @staticmethod - def real_unembed(M): - """ - The inverse of :meth:`real_embed`. - """ - raise NotImplementedError + sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 - @staticmethod - def jordan_product(X,Y): - return (X*Y + Y*X)/2 + :: - @classmethod - def trace_inner_product(cls,X,Y): - Xu = cls.real_unembed(X) - Yu = cls.real_unembed(Y) - tr = (Xu*Yu).trace() + sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 - try: - # Works in QQ, AA, RDF, et cetera. - return tr.real() - except AttributeError: - # A quaternion doesn't have a real() method, but does - # have coefficient_tuple() method that returns the - # coefficients of 1, i, j, and k -- in that order. - return tr.coefficient_tuple()[0] + :: + sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 -class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): - @staticmethod - def real_embed(M): - """ - The identity function, for embedding real matrices into real - matrices. - """ - return M + :: + + sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False) + sage: I = J.one().to_matrix() + sage: J.trace_inner_product(I, -I) + -2 - @staticmethod - def real_unembed(M): - """ - The identity function, for unembedding real matrices from real - matrices. """ - return M + tr = (X*Y).trace() + if hasattr(tr, 'coefficient'): + # Works for octonions, and has to come first because they + # also have a "real()" method that doesn't return an + # element of the scalar ring. + return tr.coefficient(0) + elif hasattr(tr, 'coefficient_tuple'): + # Works for quaternions. + return tr.coefficient_tuple()[0] + + # Works for real and complex numbers. + return tr.real() + + + def __init__(self, matrix_space, **kwargs): + # We know this is a valid EJA, but will double-check + # if the user passes check_axioms=True. + if "check_axioms" not in kwargs: kwargs["check_axioms"] = False + super().__init__(self._denormalized_basis(matrix_space), + self.jordan_product, + self.trace_inner_product, + field=matrix_space.base_ring(), + matrix_space=matrix_space, + **kwargs) -class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra, - RealMatrixEuclideanJordanAlgebra): + self.rank.set_cache(matrix_space.nrows()) + self.one.set_cache( self(matrix_space.one()) ) + +class RealSymmetricEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1522,19 +2081,19 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra, EXAMPLES:: sage: J = RealSymmetricEJA(2) - sage: e0, e1, e2 = J.gens() - sage: e0*e0 - e0 - sage: e1*e1 - 1/2*e0 + 1/2*e2 - sage: e2*e2 - e2 + sage: b0, b1, b2 = J.gens() + sage: b0*b0 + b0 + sage: b1*b1 + 1/2*b0 + 1/2*b2 + sage: b2*b2 + b2 In theory, our "field" can be any subfield of the reals:: - sage: RealSymmetricEJA(2, RDF) + sage: RealSymmetricEJA(2, field=RDF, check_axioms=True) Euclidean Jordan algebra of dimension 3 over Real Double Field - sage: RealSymmetricEJA(2, RR) + sage: RealSymmetricEJA(2, field=RR, check_axioms=True) Euclidean Jordan algebra of dimension 3 over Real Field with 53 bits of precision @@ -1542,16 +2101,14 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra, The dimension of this algebra is `(n^2 + n) / 2`:: - sage: set_random_seed() - sage: n_max = RealSymmetricEJA._max_random_instance_size() - sage: n = ZZ.random_element(1, n_max) + sage: d = RealSymmetricEJA._max_random_instance_dimension() + sage: n = RealSymmetricEJA._max_random_instance_size(d) sage: J = RealSymmetricEJA(n) sage: J.dimension() == (n^2 + n)/2 True The Jordan multiplication is what we think it is:: - sage: set_random_seed() sage: J = RealSymmetricEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).to_matrix() @@ -1574,221 +2131,36 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra, Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ - @classmethod - def _denormalized_basis(cls, n, field): - """ - Return a basis for the space of real symmetric n-by-n matrices. - - SETUP:: - - sage: from mjo.eja.eja_algebra import RealSymmetricEJA - - TESTS:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: B = RealSymmetricEJA._denormalized_basis(n,QQ) - sage: all( M.is_symmetric() for M in B) - True - - """ - # The basis of symmetric matrices, as matrices, in their R^(n-by-n) - # coordinates. - S = [] - for i in range(n): - for j in range(i+1): - Eij = matrix(field, n, lambda k,l: k==i and l==j) - if i == j: - Sij = Eij - else: - Sij = Eij + Eij.transpose() - S.append(Sij) - return tuple(S) - - @staticmethod - def _max_random_instance_size(): - return 4 # Dimension 10 + def _max_random_instance_size(max_dimension): + # Obtained by solving d = (n^2 + n)/2. + # The ZZ-int-ZZ thing is just "floor." + return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/2 - 1/2)) @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this type of algebra. """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) - return cls(n, field, **kwargs) + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) + return cls(n, **kwargs) def __init__(self, n, field=AA, **kwargs): - basis = self._denormalized_basis(n, field) - super(RealSymmetricEJA, self).__init__(field, - basis, - self.jordan_product, - self.trace_inner_product, - **kwargs) - self.rank.set_cache(n) - self.one.set_cache(self(matrix.identity(field,n))) - - -class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): - @staticmethod - def real_embed(M): - """ - Embed the n-by-n complex matrix ``M`` into the space of real - matrices of size 2n-by-2n via the map the sends each entry `z = a + - bi` to the block matrix ``[[a,b],[-b,a]]``. - - SETUP:: - - sage: from mjo.eja.eja_algebra import \ - ....: ComplexMatrixEuclideanJordanAlgebra - - EXAMPLES:: - - sage: F = QuadraticField(-1, 'I') - sage: x1 = F(4 - 2*i) - sage: x2 = F(1 + 2*i) - sage: x3 = F(-i) - sage: x4 = F(6) - sage: M = matrix(F,2,[[x1,x2],[x3,x4]]) - sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M) - [ 4 -2| 1 2] - [ 2 4|-2 1] - [-----+-----] - [ 0 -1| 6 0] - [ 1 0| 0 6] - - TESTS: - - Embedding is a homomorphism (isomorphism, in fact):: - - sage: set_random_seed() - sage: n = ZZ.random_element(3) - sage: F = QuadraticField(-1, 'I') - sage: X = random_matrix(F, n) - sage: Y = random_matrix(F, n) - sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X) - sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y) - sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y) - sage: Xe*Ye == XYe - True - - """ - n = M.nrows() - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - - # We don't need any adjoined elements... - field = M.base_ring().base_ring() - - blocks = [] - for z in M.list(): - a = z.list()[0] # real part, I guess - b = z.list()[1] # imag part, I guess - blocks.append(matrix(field, 2, [[a,b],[-b,a]])) - - return matrix.block(field, n, blocks) - - - @staticmethod - def real_unembed(M): - """ - The inverse of _embed_complex_matrix(). - - SETUP:: - - sage: from mjo.eja.eja_algebra import \ - ....: ComplexMatrixEuclideanJordanAlgebra - - EXAMPLES:: - - sage: A = matrix(QQ,[ [ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [ 9, 10, 11, 12], - ....: [-10, 9, -12, 11] ]) - sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A) - [ 2*I + 1 4*I + 3] - [ 10*I + 9 12*I + 11] - - TESTS: - - Unembedding is the inverse of embedding:: - - sage: set_random_seed() - sage: F = QuadraticField(-1, 'I') - sage: M = random_matrix(F, 3) - sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M) - sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M - True - - """ - n = ZZ(M.nrows()) - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - if not n.mod(2).is_zero(): - raise ValueError("the matrix 'M' must be a complex embedding") - - # If "M" was normalized, its base ring might have roots - # adjoined and they can stick around after unembedding. - field = M.base_ring() - R = PolynomialRing(field, 'z') - z = R.gen() - if field is AA: - # Sage doesn't know how to embed AA into QQbar, i.e. how - # to adjoin sqrt(-1) to AA. - F = QQbar - else: - F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) - i = F.gen() - - # Go top-left to bottom-right (reading order), converting every - # 2-by-2 block we see to a single complex element. - elements = [] - for k in range(n/2): - for j in range(n/2): - submat = M[2*k:2*k+2,2*j:2*j+2] - if submat[0,0] != submat[1,1]: - raise ValueError('bad on-diagonal submatrix') - if submat[0,1] != -submat[1,0]: - raise ValueError('bad off-diagonal submatrix') - z = submat[0,0] + submat[0,1]*i - elements.append(z) - - return matrix(F, n/2, elements) - - - @classmethod - def trace_inner_product(cls,X,Y): - """ - Compute a matrix inner product in this algebra directly from - its real embedding. + A = MatrixSpace(field, n) + super().__init__(A, **kwargs) - SETUP:: - - sage: from mjo.eja.eja_algebra import ComplexHermitianEJA - - TESTS: - - This gives the same answer as the slow, default method implemented - in :class:`MatrixEuclideanJordanAlgebra`:: - - sage: set_random_seed() - sage: J = ComplexHermitianEJA.random_instance() - sage: x,y = J.random_elements(2) - sage: Xe = x.to_matrix() - sage: Ye = y.to_matrix() - sage: X = ComplexHermitianEJA.real_unembed(Xe) - sage: Y = ComplexHermitianEJA.real_unembed(Ye) - sage: expected = (X*Y).trace().real() - sage: actual = ComplexHermitianEJA.trace_inner_product(Xe,Ye) - sage: actual == expected - True + from mjo.eja.eja_cache import real_symmetric_eja_coeffs + a = real_symmetric_eja_coeffs(self) + if a is not None: + self.rational_algebra()._charpoly_coefficients.set_cache(a) - """ - return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2 -class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, - ComplexMatrixEuclideanJordanAlgebra): +class ComplexHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -1801,28 +2173,32 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, EXAMPLES: - In theory, our "field" can be any subfield of the reals:: + In theory, our "field" can be any subfield of the reals, but we + can't use inexact real fields at the moment because SageMath + doesn't know how to convert their elements into complex numbers, + or even into algebraic reals:: - sage: ComplexHermitianEJA(2, RDF) - Euclidean Jordan algebra of dimension 4 over Real Double Field - sage: ComplexHermitianEJA(2, RR) - Euclidean Jordan algebra of dimension 4 over Real Field with - 53 bits of precision + sage: QQbar(RDF(1)) + Traceback (most recent call last): + ... + TypeError: Illegal initializer for algebraic number + sage: AA(RR(1)) + Traceback (most recent call last): + ... + TypeError: Illegal initializer for algebraic number TESTS: The dimension of this algebra is `n^2`:: - sage: set_random_seed() - sage: n_max = ComplexHermitianEJA._max_random_instance_size() - sage: n = ZZ.random_element(1, n_max) + sage: d = ComplexHermitianEJA._max_random_instance_dimension() + sage: n = ComplexHermitianEJA._max_random_instance_size(d) sage: J = ComplexHermitianEJA(n) sage: J.dimension() == n^2 True The Jordan multiplication is what we think it is:: - sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).to_matrix() @@ -1845,247 +2221,36 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra, Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ - - @classmethod - def _denormalized_basis(cls, n, field): - """ - Returns a basis for the space of complex Hermitian n-by-n matrices. - - Why do we embed these? Basically, because all of numerical linear - algebra assumes that you're working with vectors consisting of `n` - entries from a field and scalars from the same field. There's no way - to tell SageMath that (for example) the vectors contain complex - numbers, while the scalar field is real. - - SETUP:: - - sage: from mjo.eja.eja_algebra import ComplexHermitianEJA - - TESTS:: - - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: field = QuadraticField(2, 'sqrt2') - sage: B = ComplexHermitianEJA._denormalized_basis(n, field) - sage: all( M.is_symmetric() for M in B) - True - - """ - R = PolynomialRing(field, 'z') - z = R.gen() - F = field.extension(z**2 + 1, 'I') - I = F.gen() - - # This is like the symmetric case, but we need to be careful: - # - # * We want conjugate-symmetry, not just symmetry. - # * The diagonal will (as a result) be real. - # - S = [] - for i in range(n): - for j in range(i+1): - Eij = matrix(F, n, lambda k,l: k==i and l==j) - if i == j: - Sij = cls.real_embed(Eij) - S.append(Sij) - else: - # The second one has a minus because it's conjugated. - Sij_real = cls.real_embed(Eij + Eij.transpose()) - S.append(Sij_real) - Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose()) - S.append(Sij_imag) - - # Since we embedded these, we can drop back to the "field" that we - # started with instead of the complex extension "F". - return tuple( s.change_ring(field) for s in S ) - - def __init__(self, n, field=AA, **kwargs): - basis = self._denormalized_basis(n,field) - super(ComplexHermitianEJA, self).__init__(field, - basis, - self.jordan_product, - self.trace_inner_product, - **kwargs) - self.rank.set_cache(n) - # TODO: pre-cache the identity! - - @staticmethod - def _max_random_instance_size(): - return 3 # Dimension 9 - - @classmethod - def random_instance(cls, field=AA, **kwargs): - """ - Return a random instance of this type of algebra. - """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) - return cls(n, field, **kwargs) - -class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): - @staticmethod - def real_embed(M): - """ - Embed the n-by-n quaternion matrix ``M`` into the space of real - matrices of size 4n-by-4n by first sending each quaternion entry `z - = a + bi + cj + dk` to the block-complex matrix ``[[a + bi, - c+di],[-c + di, a-bi]]`, and then embedding those into a real - matrix. - - SETUP:: - - sage: from mjo.eja.eja_algebra import \ - ....: QuaternionMatrixEuclideanJordanAlgebra - - EXAMPLES:: - - sage: Q = QuaternionAlgebra(QQ,-1,-1) - sage: i,j,k = Q.gens() - sage: x = 1 + 2*i + 3*j + 4*k - sage: M = matrix(Q, 1, [[x]]) - sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M) - [ 1 2 3 4] - [-2 1 -4 3] - [-3 4 1 -2] - [-4 -3 2 1] - - Embedding is a homomorphism (isomorphism, in fact):: - - sage: set_random_seed() - sage: n = ZZ.random_element(2) - sage: Q = QuaternionAlgebra(QQ,-1,-1) - sage: X = random_matrix(Q, n) - sage: Y = random_matrix(Q, n) - sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X) - sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y) - sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y) - sage: Xe*Ye == XYe - True - - """ - quaternions = M.base_ring() - n = M.nrows() - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - - F = QuadraticField(-1, 'I') - i = F.gen() - - blocks = [] - for z in M.list(): - t = z.coefficient_tuple() - a = t[0] - b = t[1] - c = t[2] - d = t[3] - cplxM = matrix(F, 2, [[ a + b*i, c + d*i], - [-c + d*i, a - b*i]]) - realM = ComplexMatrixEuclideanJordanAlgebra.real_embed(cplxM) - blocks.append(realM) - - # We should have real entries by now, so use the realest field - # we've got for the return value. - return matrix.block(quaternions.base_ring(), n, blocks) - + from mjo.hurwitz import ComplexMatrixAlgebra + A = ComplexMatrixAlgebra(n, scalars=field) + super().__init__(A, **kwargs) + from mjo.eja.eja_cache import complex_hermitian_eja_coeffs + a = complex_hermitian_eja_coeffs(self) + if a is not None: + self.rational_algebra()._charpoly_coefficients.set_cache(a) @staticmethod - def real_unembed(M): - """ - The inverse of _embed_quaternion_matrix(). - - SETUP:: - - sage: from mjo.eja.eja_algebra import \ - ....: QuaternionMatrixEuclideanJordanAlgebra - - EXAMPLES:: - - sage: M = matrix(QQ, [[ 1, 2, 3, 4], - ....: [-2, 1, -4, 3], - ....: [-3, 4, 1, -2], - ....: [-4, -3, 2, 1]]) - sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M) - [1 + 2*i + 3*j + 4*k] - - TESTS: - - Unembedding is the inverse of embedding:: - - sage: set_random_seed() - sage: Q = QuaternionAlgebra(QQ, -1, -1) - sage: M = random_matrix(Q, 3) - sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M) - sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M - True - - """ - n = ZZ(M.nrows()) - if M.ncols() != n: - raise ValueError("the matrix 'M' must be square") - if not n.mod(4).is_zero(): - raise ValueError("the matrix 'M' must be a quaternion embedding") - - # Use the base ring of the matrix to ensure that its entries can be - # multiplied by elements of the quaternion algebra. - field = M.base_ring() - Q = QuaternionAlgebra(field,-1,-1) - i,j,k = Q.gens() - - # Go top-left to bottom-right (reading order), converting every - # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1 - # quaternion block. - elements = [] - for l in range(n/4): - for m in range(n/4): - submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed( - M[4*l:4*l+4,4*m:4*m+4] ) - if submat[0,0] != submat[1,1].conjugate(): - raise ValueError('bad on-diagonal submatrix') - if submat[0,1] != -submat[1,0].conjugate(): - raise ValueError('bad off-diagonal submatrix') - z = submat[0,0].real() - z += submat[0,0].imag()*i - z += submat[0,1].real()*j - z += submat[0,1].imag()*k - elements.append(z) - - return matrix(Q, n/4, elements) - + def _max_random_instance_size(max_dimension): + # Obtained by solving d = n^2. + # The ZZ-int-ZZ thing is just "floor." + return ZZ(int(ZZ(max_dimension).sqrt())) @classmethod - def trace_inner_product(cls,X,Y): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ - Compute a matrix inner product in this algebra directly from - its real embedding. - - SETUP:: - - sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA - - TESTS: - - This gives the same answer as the slow, default method implemented - in :class:`MatrixEuclideanJordanAlgebra`:: - - sage: set_random_seed() - sage: J = QuaternionHermitianEJA.random_instance() - sage: x,y = J.random_elements(2) - sage: Xe = x.to_matrix() - sage: Ye = y.to_matrix() - sage: X = QuaternionHermitianEJA.real_unembed(Xe) - sage: Y = QuaternionHermitianEJA.real_unembed(Ye) - sage: expected = (X*Y).trace().coefficient_tuple()[0] - sage: actual = QuaternionHermitianEJA.trace_inner_product(Xe,Ye) - sage: actual == expected - True - + Return a random instance of this type of algebra. """ - return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4 + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) + return cls(n, **kwargs) -class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, - QuaternionMatrixEuclideanJordanAlgebra): +class QuaternionHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA): r""" The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the @@ -2100,9 +2265,9 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, In theory, our "field" can be any subfield of the reals:: - sage: QuaternionHermitianEJA(2, RDF) + sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True) Euclidean Jordan algebra of dimension 6 over Real Double Field - sage: QuaternionHermitianEJA(2, RR) + sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True) Euclidean Jordan algebra of dimension 6 over Real Field with 53 bits of precision @@ -2110,16 +2275,14 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, The dimension of this algebra is `2*n^2 - n`:: - sage: set_random_seed() - sage: n_max = QuaternionHermitianEJA._max_random_instance_size() - sage: n = ZZ.random_element(1, n_max) + sage: d = QuaternionHermitianEJA._max_random_instance_dimension() + sage: n = QuaternionHermitianEJA._max_random_instance_size(d) sage: J = QuaternionHermitianEJA(n) sage: J.dimension() == 2*(n^2) - n True The Jordan multiplication is what we think it is:: - sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).to_matrix() @@ -2142,96 +2305,203 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra, Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ + def __init__(self, n, field=AA, **kwargs): + from mjo.hurwitz import QuaternionMatrixAlgebra + A = QuaternionMatrixAlgebra(n, scalars=field) + super().__init__(A, **kwargs) + + from mjo.eja.eja_cache import quaternion_hermitian_eja_coeffs + a = quaternion_hermitian_eja_coeffs(self) + if a is not None: + self.rational_algebra()._charpoly_coefficients.set_cache(a) + + + + @staticmethod + def _max_random_instance_size(max_dimension): + r""" + The maximum rank of a random QuaternionHermitianEJA. + """ + # Obtained by solving d = 2n^2 - n. + # The ZZ-int-ZZ thing is just "floor." + return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4)) + @classmethod - def _denormalized_basis(cls, n, field): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ - Returns a basis for the space of quaternion Hermitian n-by-n matrices. + Return a random instance of this type of algebra. + """ + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) + return cls(n, **kwargs) + +class OctonionHermitianEJA(HermitianMatrixEJA, RationalBasisEJA, ConcreteEJA): + r""" + SETUP:: - Why do we embed these? Basically, because all of numerical - linear algebra assumes that you're working with vectors consisting - of `n` entries from a field and scalars from the same field. There's - no way to tell SageMath that (for example) the vectors contain - complex numbers, while the scalar field is real. + sage: from mjo.eja.eja_algebra import (EJA, + ....: OctonionHermitianEJA) + sage: from mjo.hurwitz import Octonions, OctonionMatrixAlgebra - SETUP:: + EXAMPLES: - sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA + The 3-by-3 algebra satisfies the axioms of an EJA:: + + sage: OctonionHermitianEJA(3, # long time + ....: field=QQ, # long time + ....: orthonormalize=False, # long time + ....: check_axioms=True) # long time + Euclidean Jordan algebra of dimension 27 over Rational Field + + After a change-of-basis, the 2-by-2 algebra has the same + multiplication table as the ten-dimensional Jordan spin algebra:: + + sage: A = OctonionMatrixAlgebra(2,Octonions(QQ),QQ) + sage: b = OctonionHermitianEJA._denormalized_basis(A) + sage: basis = (b[0] + b[9],) + b[1:9] + (b[0] - b[9],) + sage: jp = OctonionHermitianEJA.jordan_product + sage: ip = OctonionHermitianEJA.trace_inner_product + sage: J = EJA(basis, + ....: jp, + ....: ip, + ....: field=QQ, + ....: orthonormalize=False) + sage: J.multiplication_table() + +----++----+----+----+----+----+----+----+----+----+----+ + | * || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 | + +====++====+====+====+====+====+====+====+====+====+====+ + | b0 || b0 | b1 | b2 | b3 | b4 | b5 | b6 | b7 | b8 | b9 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b1 || b1 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b2 || b2 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b3 || b3 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b4 || b4 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b5 || b5 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b6 || b6 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b7 || b7 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b8 || b8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | 0 | + +----++----+----+----+----+----+----+----+----+----+----+ + | b9 || b9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | b0 | + +----++----+----+----+----+----+----+----+----+----+----+ - TESTS:: + TESTS: - sage: set_random_seed() - sage: n = ZZ.random_element(1,5) - sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ) - sage: all( M.is_symmetric() for M in B ) - True + We can actually construct the 27-dimensional Albert algebra, + and we get the right unit element if we recompute it:: + + sage: J = OctonionHermitianEJA(3, # long time + ....: field=QQ, # long time + ....: orthonormalize=False) # long time + sage: J.one.clear_cache() # long time + sage: J.one() # long time + b0 + b9 + b26 + sage: J.one().to_matrix() # long time + +----+----+----+ + | e0 | 0 | 0 | + +----+----+----+ + | 0 | e0 | 0 | + +----+----+----+ + | 0 | 0 | e0 | + +----+----+----+ + + The 2-by-2 algebra is isomorphic to the ten-dimensional Jordan + spin algebra, but just to be sure, we recompute its rank:: + + sage: J = OctonionHermitianEJA(2, # long time + ....: field=QQ, # long time + ....: orthonormalize=False) # long time + sage: J.rank.clear_cache() # long time + sage: J.rank() # long time + 2 + + """ + @staticmethod + def _max_random_instance_size(max_dimension): + r""" + The maximum rank of a random OctonionHermitianEJA. + """ + # There's certainly a formula for this, but with only four + # cases to worry about, I'm not that motivated to derive it. + if max_dimension >= 27: + return 3 + elif max_dimension >= 10: + return 2 + elif max_dimension >= 1: + return 1 + else: + return 0 + @classmethod + def random_instance(cls, max_dimension=None, *args, **kwargs): + """ + Return a random instance of this type of algebra. """ - Q = QuaternionAlgebra(QQ,-1,-1) - I,J,K = Q.gens() + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) + return cls(n, **kwargs) - # This is like the symmetric case, but we need to be careful: - # - # * We want conjugate-symmetry, not just symmetry. - # * The diagonal will (as a result) be real. - # - S = [] - for i in range(n): - for j in range(i+1): - Eij = matrix(Q, n, lambda k,l: k==i and l==j) - if i == j: - Sij = cls.real_embed(Eij) - S.append(Sij) - else: - # The second, third, and fourth ones have a minus - # because they're conjugated. - Sij_real = cls.real_embed(Eij + Eij.transpose()) - S.append(Sij_real) - Sij_I = cls.real_embed(I*Eij - I*Eij.transpose()) - S.append(Sij_I) - Sij_J = cls.real_embed(J*Eij - J*Eij.transpose()) - S.append(Sij_J) - Sij_K = cls.real_embed(K*Eij - K*Eij.transpose()) - S.append(Sij_K) - - # Since we embedded these, we can drop back to the "field" that we - # started with instead of the quaternion algebra "Q". - return tuple( s.change_ring(field) for s in S ) + def __init__(self, n, field=AA, **kwargs): + if n > 3: + # Otherwise we don't get an EJA. + raise ValueError("n cannot exceed 3") + # We know this is a valid EJA, but will double-check + # if the user passes check_axioms=True. + if "check_axioms" not in kwargs: kwargs["check_axioms"] = False - def __init__(self, n, field=AA, **kwargs): - basis = self._denormalized_basis(n,field) - super(QuaternionHermitianEJA, self).__init__(field, - basis, - self.jordan_product, - self.trace_inner_product, - **kwargs) - self.rank.set_cache(n) - # TODO: cache one()! + from mjo.hurwitz import OctonionMatrixAlgebra + A = OctonionMatrixAlgebra(n, scalars=field) + super().__init__(A, **kwargs) + + from mjo.eja.eja_cache import octonion_hermitian_eja_coeffs + a = octonion_hermitian_eja_coeffs(self) + if a is not None: + self.rational_algebra()._charpoly_coefficients.set_cache(a) + + +class AlbertEJA(OctonionHermitianEJA): + r""" + The Albert algebra is the algebra of three-by-three Hermitian + matrices whose entries are octonions. + + SETUP:: + + sage: from mjo.eja.eja_algebra import AlbertEJA + + EXAMPLES:: - @staticmethod - def _max_random_instance_size(): - r""" - The maximum rank of a random QuaternionHermitianEJA. - """ - return 2 # Dimension 6 + sage: AlbertEJA(field=QQ, orthonormalize=False) + Euclidean Jordan algebra of dimension 27 over Rational Field + sage: AlbertEJA() # long time + Euclidean Jordan algebra of dimension 27 over Algebraic Real Field - @classmethod - def random_instance(cls, field=AA, **kwargs): - """ - Return a random instance of this type of algebra. - """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) - return cls(n, field, **kwargs) + """ + def __init__(self, *args, **kwargs): + super().__init__(3, *args, **kwargs) -class HadamardEJA(ConcreteEuclideanJordanAlgebra): +class HadamardEJA(RationalBasisEJA, ConcreteEJA): """ - Return the Euclidean Jordan Algebra corresponding to the set - `R^n` under the Hadamard product. + Return the Euclidean Jordan algebra on `R^n` with the Hadamard + (pointwise real-number multiplication) Jordan product and the + usual inner-product. - Note: this is nothing more than the Cartesian product of ``n`` - copies of the spin algebra. Once Cartesian product algebras - are implemented, this can go. + This is nothing more than the Cartesian product of ``n`` copies of + the one-dimensional Jordan spin algebra, and is the most common + example of a non-simple Euclidean Jordan algebra. SETUP:: @@ -2242,19 +2512,19 @@ class HadamardEJA(ConcreteEuclideanJordanAlgebra): This multiplication table can be verified by hand:: sage: J = HadamardEJA(3) - sage: e0,e1,e2 = J.gens() - sage: e0*e0 - e0 - sage: e0*e1 + sage: b0,b1,b2 = J.gens() + sage: b0*b0 + b0 + sage: b0*b1 0 - sage: e0*e2 + sage: b0*b2 0 - sage: e1*e1 - e1 - sage: e1*e2 + sage: b1*b1 + b1 + sage: b1*b2 0 - sage: e2*e2 - e2 + sage: b2*b2 + b2 TESTS: @@ -2262,46 +2532,70 @@ class HadamardEJA(ConcreteEuclideanJordanAlgebra): sage: HadamardEJA(3, prefix='r').gens() (r0, r1, r2) - """ def __init__(self, n, field=AA, **kwargs): - V = VectorSpace(field, n) - basis = V.basis() - - def jordan_product(x,y): - return V([ xi*yi for (xi,yi) in zip(x,y) ]) - def inner_product(x,y): - return x.inner_product(y) - - super(HadamardEJA, self).__init__(field, - basis, - jordan_product, - inner_product, - **kwargs) - self.rank.set_cache(n) + MS = MatrixSpace(field, n, 1) if n == 0: - self.one.set_cache( self.zero() ) + jordan_product = lambda x,y: x + inner_product = lambda x,y: x else: - self.one.set_cache( sum(self.gens()) ) + def jordan_product(x,y): + return MS( xi*yi for (xi,yi) in zip(x,y) ) + + def inner_product(x,y): + return (x.T*y)[0,0] + + # New defaults for keyword arguments. Don't orthonormalize + # because our basis is already orthonormal with respect to our + # inner-product. Don't check the axioms, because we know this + # is a valid EJA... but do double-check if the user passes + # check_axioms=True. Note: we DON'T override the "check_field" + # default here, because the user can pass in a field! + if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False + if "check_axioms" not in kwargs: kwargs["check_axioms"] = False + + column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() ) + super().__init__(column_basis, + jordan_product, + inner_product, + field=field, + matrix_space=MS, + associative=True, + **kwargs) + self.rank.set_cache(n) + + self.one.set_cache( self.sum(self.gens()) ) @staticmethod - def _max_random_instance_size(): + def _max_random_instance_dimension(): r""" - The maximum dimension of a random HadamardEJA. + There's no reason to go higher than five here. That's + enough to get the point across. """ return 5 + @staticmethod + def _max_random_instance_size(max_dimension): + r""" + The maximum size (=dimension) of a random HadamardEJA. + """ + return max_dimension + @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this type of algebra. """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) - return cls(n, field, **kwargs) + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) + return cls(n, **kwargs) -class BilinearFormEJA(ConcreteEuclideanJordanAlgebra): +class BilinearFormEJA(RationalBasisEJA, ConcreteEJA): r""" The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` with the half-trace inner product and jordan product ``x*y = @@ -2361,7 +2655,6 @@ class BilinearFormEJA(ConcreteEuclideanJordanAlgebra): matrix. We opt not to orthonormalize the basis, because if we did, we would have to normalize the `s_{i}` in a similar manner:: - sage: set_random_seed() sage: n = ZZ.random_element(5) sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular') sage: B11 = matrix.identity(QQ,1) @@ -2381,72 +2674,105 @@ class BilinearFormEJA(ConcreteEuclideanJordanAlgebra): ....: for j in range(n-1) ] sage: actual == expected True + """ def __init__(self, B, field=AA, **kwargs): - if not B.is_positive_definite(): - raise ValueError("bilinear form is not positive-definite") + # The matrix "B" is supplied by the user in most cases, + # so it makes sense to check whether or not its positive- + # definite unless we are specifically asked not to... + if ("check_axioms" not in kwargs) or kwargs["check_axioms"]: + if not B.is_positive_definite(): + raise ValueError("bilinear form is not positive-definite") + + # However, all of the other data for this EJA is computed + # by us in manner that guarantees the axioms are + # satisfied. So, again, unless we are specifically asked to + # verify things, we'll skip the rest of the checks. + if "check_axioms" not in kwargs: kwargs["check_axioms"] = False n = B.nrows() - V = VectorSpace(field, n) + MS = MatrixSpace(field, n, 1) def inner_product(x,y): - return (B*x).inner_product(y) + return (y.T*B*x)[0,0] def jordan_product(x,y): - x0 = x[0] - xbar = x[1:] - y0 = y[0] - ybar = y[1:] - z0 = inner_product(x,y) + x0 = x[0,0] + xbar = x[1:,0] + y0 = y[0,0] + ybar = y[1:,0] + z0 = inner_product(y,x) zbar = y0*xbar + x0*ybar - return V([z0] + zbar.list()) + return MS([z0] + zbar.list()) - super(BilinearFormEJA, self).__init__(field, - V.basis(), - jordan_product, - inner_product, - **kwargs) + column_basis = tuple( MS(b) for b in FreeModule(field, n).basis() ) + + # TODO: I haven't actually checked this, but it seems legit. + associative = False + if n <= 2: + associative = True + + super().__init__(column_basis, + jordan_product, + inner_product, + field=field, + matrix_space=MS, + associative=associative, + **kwargs) # The rank of this algebra is two, unless we're in a # one-dimensional ambient space (because the rank is bounded # by the ambient dimension). self.rank.set_cache(min(n,2)) - if n == 0: self.one.set_cache( self.zero() ) else: self.one.set_cache( self.monomial(0) ) @staticmethod - def _max_random_instance_size(): + def _max_random_instance_dimension(): r""" - The maximum dimension of a random BilinearFormEJA. + There's no reason to go higher than five here. That's + enough to get the point across. """ return 5 + @staticmethod + def _max_random_instance_size(max_dimension): + r""" + The maximum size (=dimension) of a random BilinearFormEJA. + """ + return max_dimension + @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this algebra. """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) + if n.is_zero(): - B = matrix.identity(field, n) - return cls(B, field, **kwargs) + B = matrix.identity(ZZ, n) + return cls(B, **kwargs) - B11 = matrix.identity(field,1) - M = matrix.random(field, n-1) - I = matrix.identity(field, n-1) - alpha = field.zero() + B11 = matrix.identity(ZZ, 1) + M = matrix.random(ZZ, n-1) + I = matrix.identity(ZZ, n-1) + alpha = ZZ.zero() while alpha.is_zero(): - alpha = field.random_element().abs() + alpha = ZZ.random_element().abs() + B22 = M.transpose()*M + alpha*I from sage.matrix.special import block_matrix B = block_matrix(2,2, [ [B11, ZZ(0) ], [ZZ(0), B22 ] ]) - return cls(B, field, **kwargs) + return cls(B, **kwargs) class JordanSpinEJA(BilinearFormEJA): @@ -2465,20 +2791,20 @@ class JordanSpinEJA(BilinearFormEJA): This multiplication table can be verified by hand:: sage: J = JordanSpinEJA(4) - sage: e0,e1,e2,e3 = J.gens() - sage: e0*e0 - e0 - sage: e0*e1 - e1 - sage: e0*e2 - e2 - sage: e0*e3 - e3 - sage: e1*e2 + sage: b0,b1,b2,b3 = J.gens() + sage: b0*b0 + b0 + sage: b0*b1 + b1 + sage: b0*b2 + b2 + sage: b0*b3 + b3 + sage: b1*b2 0 - sage: e1*e3 + sage: b1*b3 0 - sage: e2*e3 + sage: b2*b3 0 We can change the generator prefix:: @@ -2490,7 +2816,6 @@ class JordanSpinEJA(BilinearFormEJA): Ensure that we have the usual inner product on `R^n`:: - sage: set_random_seed() sage: J = JordanSpinEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = x.inner_product(y) @@ -2499,31 +2824,35 @@ class JordanSpinEJA(BilinearFormEJA): True """ - def __init__(self, n, field=AA, **kwargs): - # This is a special case of the BilinearFormEJA with the identity - # matrix as its bilinear form. - B = matrix.identity(field, n) - super(JordanSpinEJA, self).__init__(B, field, **kwargs) + def __init__(self, n, *args, **kwargs): + # This is a special case of the BilinearFormEJA with the + # identity matrix as its bilinear form. + B = matrix.identity(ZZ, n) - @staticmethod - def _max_random_instance_size(): - r""" - The maximum dimension of a random JordanSpinEJA. - """ - return 5 + # Don't orthonormalize because our basis is already + # orthonormal with respect to our inner-product. + if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False + + # But also don't pass check_field=False here, because the user + # can pass in a field! + super().__init__(B, *args, **kwargs) @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this type of algebra. Needed here to override the implementation for ``BilinearFormEJA``. """ - n = ZZ.random_element(cls._max_random_instance_size() + 1) - return cls(n, field, **kwargs) + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) + return cls(n, **kwargs) -class TrivialEJA(ConcreteEuclideanJordanAlgebra): +class TrivialEJA(RationalBasisEJA, ConcreteEJA): """ The trivial Euclidean Jordan algebra consisting of only a zero element. @@ -2554,173 +2883,483 @@ class TrivialEJA(ConcreteEuclideanJordanAlgebra): """ def __init__(self, field=AA, **kwargs): jordan_product = lambda x,y: x - inner_product = lambda x,y: field(0) + inner_product = lambda x,y: field.zero() basis = () - super(TrivialEJA, self).__init__(field, - basis, - jordan_product, - inner_product, - **kwargs) + MS = MatrixSpace(field,0) + + # New defaults for keyword arguments + if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False + if "check_axioms" not in kwargs: kwargs["check_axioms"] = False + + super().__init__(basis, + jordan_product, + inner_product, + associative=True, + field=field, + matrix_space=MS, + **kwargs) + # The rank is zero using my definition, namely the dimension of the # largest subalgebra generated by any element. self.rank.set_cache(0) self.one.set_cache( self.zero() ) @classmethod - def random_instance(cls, field=AA, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): # We don't take a "size" argument so the superclass method is - # inappropriate for us. - return cls(field, **kwargs) + # inappropriate for us. The ``max_dimension`` argument is + # included so that if this method is called generically with a + # ``max_dimension=`` argument, we don't try to pass + # it on to the algebra constructor. + return cls(**kwargs) -class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): + +class CartesianProductEJA(EJA): r""" - The external (orthogonal) direct sum of two other Euclidean Jordan - algebras. Essentially the Cartesian product of its two factors. - Every Euclidean Jordan algebra decomposes into an orthogonal - direct sum of simple Euclidean Jordan algebras, so no generality - is lost by providing only this construction. + The external (orthogonal) direct sum of two or more Euclidean + Jordan algebras. Every Euclidean Jordan algebra decomposes into an + orthogonal direct sum of simple Euclidean Jordan algebras which is + then isometric to a Cartesian product, so no generality is lost by + providing only this construction. SETUP:: sage: from mjo.eja.eja_algebra import (random_eja, + ....: CartesianProductEJA, + ....: ComplexHermitianEJA, ....: HadamardEJA, - ....: RealSymmetricEJA, - ....: DirectSumEJA) + ....: JordanSpinEJA, + ....: RealSymmetricEJA) - EXAMPLES:: + EXAMPLES: + + The Jordan product is inherited from our factors and implemented by + our CombinatorialFreeModule Cartesian product superclass:: sage: J1 = HadamardEJA(2) - sage: J2 = RealSymmetricEJA(3) - sage: J = DirectSumEJA(J1,J2) - sage: J.dimension() - 8 + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: x,y = J.random_elements(2) + sage: x*y in J + True + + The ability to retrieve the original factors is implemented by our + CombinatorialFreeModule Cartesian product superclass:: + + sage: J1 = HadamardEJA(2, field=QQ) + sage: J2 = JordanSpinEJA(3, field=QQ) + sage: J = cartesian_product([J1,J2]) + sage: J.cartesian_factors() + (Euclidean Jordan algebra of dimension 2 over Rational Field, + Euclidean Jordan algebra of dimension 3 over Rational Field) + + You can provide more than two factors:: + + sage: J1 = HadamardEJA(2) + sage: J2 = JordanSpinEJA(3) + sage: J3 = RealSymmetricEJA(3) + sage: cartesian_product([J1,J2,J3]) + Euclidean Jordan algebra of dimension 2 over Algebraic Real + Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic + Real Field (+) Euclidean Jordan algebra of dimension 6 over + Algebraic Real Field + + Rank is additive on a Cartesian product:: + + sage: J1 = HadamardEJA(1) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J1.rank.clear_cache() + sage: J2.rank.clear_cache() + sage: J.rank.clear_cache() sage: J.rank() - 5 + 3 + sage: J.rank() == J1.rank() + J2.rank() + True + + The same rank computation works over the rationals, with whatever + basis you like:: + + sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False) + sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False) + sage: J = cartesian_product([J1,J2]) + sage: J1.rank.clear_cache() + sage: J2.rank.clear_cache() + sage: J.rank.clear_cache() + sage: J.rank() + 3 + sage: J.rank() == J1.rank() + J2.rank() + True + + The product algebra will be associative if and only if all of its + components are associative:: + + sage: J1 = HadamardEJA(2) + sage: J1.is_associative() + True + sage: J2 = HadamardEJA(3) + sage: J2.is_associative() + True + sage: J3 = RealSymmetricEJA(3) + sage: J3.is_associative() + False + sage: CP1 = cartesian_product([J1,J2]) + sage: CP1.is_associative() + True + sage: CP2 = cartesian_product([J1,J3]) + sage: CP2.is_associative() + False + + Cartesian products of Cartesian products work:: + + sage: J1 = JordanSpinEJA(1) + sage: J2 = JordanSpinEJA(1) + sage: J3 = JordanSpinEJA(1) + sage: J = cartesian_product([J1,cartesian_product([J2,J3])]) + sage: J.multiplication_table() + +----++----+----+----+ + | * || b0 | b1 | b2 | + +====++====+====+====+ + | b0 || b0 | 0 | 0 | + +----++----+----+----+ + | b1 || 0 | b1 | 0 | + +----++----+----+----+ + | b2 || 0 | 0 | b2 | + +----++----+----+----+ + sage: HadamardEJA(3).multiplication_table() + +----++----+----+----+ + | * || b0 | b1 | b2 | + +====++====+====+====+ + | b0 || b0 | 0 | 0 | + +----++----+----+----+ + | b1 || 0 | b1 | 0 | + +----++----+----+----+ + | b2 || 0 | 0 | b2 | + +----++----+----+----+ + + The "matrix space" of a Cartesian product always consists of + ordered pairs (or triples, or...) whose components are the + matrix spaces of its factors:: + + sage: J1 = HadamardEJA(2) + sage: J2 = ComplexHermitianEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.matrix_space() + The Cartesian product of (Full MatrixSpace of 2 by 1 dense + matrices over Algebraic Real Field, Module of 2 by 2 matrices + with entries in Algebraic Field over the scalar ring Algebraic + Real Field) + sage: J.one().to_matrix()[0] + [1] + [1] + sage: J.one().to_matrix()[1] + +---+---+ + | 1 | 0 | + +---+---+ + | 0 | 1 | + +---+---+ TESTS: - The external direct sum construction is only valid when the two factors - have the same base ring; an error is raised otherwise:: + All factors must share the same base field:: - sage: set_random_seed() - sage: J1 = random_eja(AA) - sage: J2 = random_eja(QQ,orthonormalize=False) - sage: J = DirectSumEJA(J1,J2) + sage: J1 = HadamardEJA(2, field=QQ) + sage: J2 = RealSymmetricEJA(2) + sage: CartesianProductEJA((J1,J2)) Traceback (most recent call last): ... - ValueError: algebras must share the same base field + ValueError: all factors must share the same base field + The cached unit element is the same one that would be computed:: + + sage: J1 = random_eja() # long time + sage: J2 = random_eja() # long time + sage: J = cartesian_product([J1,J2]) # long time + sage: actual = J.one() # long time + sage: J.one.clear_cache() # long time + sage: expected = J.one() # long time + sage: actual == expected # long time + True """ - def __init__(self, J1, J2, **kwargs): - if J1.base_ring() != J2.base_ring(): - raise ValueError("algebras must share the same base field") - field = J1.base_ring() - - self._factors = (J1, J2) - n1 = J1.dimension() - n2 = J2.dimension() - n = n1+n2 - V = VectorSpace(field, n) - mult_table = [ [ V.zero() for j in range(i+1) ] - for i in range(n) ] - for i in range(n1): - for j in range(i+1): - p = (J1.monomial(i)*J1.monomial(j)).to_vector() - mult_table[i][j] = V(p.list() + [field.zero()]*n2) + Element = CartesianProductEJAElement + def __init__(self, factors, **kwargs): + m = len(factors) + if m == 0: + return TrivialEJA() + + self._sets = factors + + field = factors[0].base_ring() + if not all( J.base_ring() == field for J in factors ): + raise ValueError("all factors must share the same base field") + + # Figure out the category to use. + associative = all( f.is_associative() for f in factors ) + category = EuclideanJordanAlgebras(field) + if associative: category = category.Associative() + category = category.join([category, category.CartesianProducts()]) + + # Compute my matrix space. We don't simply use the + # ``cartesian_product()`` functor here because it acts + # differently on SageMath MatrixSpaces and our custom + # MatrixAlgebras, which are CombinatorialFreeModules. We + # always want the result to be represented (and indexed) as an + # ordered tuple. This category isn't perfect, but is good + # enough for what we need to do. + MS_cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis() + MS_cat = MS_cat.Unital().CartesianProducts() + MS_factors = tuple( J.matrix_space() for J in factors ) + from sage.sets.cartesian_product import CartesianProduct + self._matrix_space = CartesianProduct(MS_factors, MS_cat) + + self._matrix_basis = [] + zero = self._matrix_space.zero() + for i in range(m): + for b in factors[i].matrix_basis(): + z = list(zero) + z[i] = b + self._matrix_basis.append(z) + + self._matrix_basis = tuple( self._matrix_space(b) + for b in self._matrix_basis ) + n = len(self._matrix_basis) + + # We already have what we need for the super-superclass constructor. + CombinatorialFreeModule.__init__(self, + field, + range(n), + prefix="b", + category=category, + bracket=False) + + # Now create the vector space for the algebra, which will have + # its own set of non-ambient coordinates (in terms of the + # supplied basis). + degree = sum( f._matrix_span.ambient_vector_space().degree() + for f in factors ) + V = VectorSpace(field, degree) + vector_basis = tuple( V(_all2list(b)) for b in self._matrix_basis ) + + # Save the span of our matrix basis (when written out as long + # vectors) because otherwise we'll have to reconstruct it + # every time we want to coerce a matrix into the algebra. + self._matrix_span = V.span_of_basis( vector_basis, check=False) + + # Since we don't (re)orthonormalize the basis, the FDEJA + # constructor is going to set self._deortho_matrix to the + # identity matrix. Here we set it to the correct value using + # the deortho matrices from our factors. + self._deortho_matrix = matrix.block_diagonal( + [J._deortho_matrix for J in factors] + ) - for i in range(n2): - for j in range(i+1): - p = (J2.monomial(i)*J2.monomial(j)).to_vector() - mult_table[n1+i][n1+j] = V([field.zero()]*n1 + p.list()) - - # TODO: build the IP table here from the two constituent IP - # matrices (it'll be block diagonal, I think). - ip_table = [ [ field.zero() for j in range(i+1) ] - for i in range(n) ] - super(DirectSumEJA, self).__init__(field, - mult_table, - ip_table, - check_axioms=False, - **kwargs) - self.rank.set_cache(J1.rank() + J2.rank()) - - - def factors(self): + self._inner_product_matrix = matrix.block_diagonal( + [J._inner_product_matrix for J in factors] + ) + self._inner_product_matrix._cache = {'hermitian': True} + self._inner_product_matrix.set_immutable() + + # Building the multiplication table is a bit more tricky + # because we have to embed the entries of the factors' + # multiplication tables into the product EJA. + zed = self.zero() + self._multiplication_table = [ [zed for j in range(i+1)] + for i in range(n) ] + + # Keep track of an offset that tallies the dimensions of all + # previous factors. If the second factor is dim=2 and if the + # first one is dim=3, then we want to skip the first 3x3 block + # when copying the multiplication table for the second factor. + offset = 0 + for f in range(m): + phi_f = self.cartesian_embedding(f) + factor_dim = factors[f].dimension() + for i in range(factor_dim): + for j in range(i+1): + f_ij = factors[f]._multiplication_table[i][j] + e = phi_f(f_ij) + self._multiplication_table[offset+i][offset+j] = e + offset += factor_dim + + self.rank.set_cache(sum(J.rank() for J in factors)) + ones = tuple(J.one().to_matrix() for J in factors) + self.one.set_cache(self(ones)) + + def _sets_keys(self): r""" - Return the pair of this algebra's factors. SETUP:: - sage: from mjo.eja.eja_algebra import (HadamardEJA, - ....: JordanSpinEJA, - ....: DirectSumEJA) + sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, + ....: RealSymmetricEJA) - EXAMPLES:: + TESTS: + + The superclass uses ``_sets_keys()`` to implement its + ``cartesian_factors()`` method:: - sage: J1 = HadamardEJA(2,QQ) - sage: J2 = JordanSpinEJA(3,QQ) - sage: J = DirectSumEJA(J1,J2) - sage: J.factors() - (Euclidean Jordan algebra of dimension 2 over Rational Field, - Euclidean Jordan algebra of dimension 3 over Rational Field) + sage: J1 = RealSymmetricEJA(2, + ....: field=QQ, + ....: orthonormalize=False, + ....: prefix="a") + sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False) + sage: J = cartesian_product([J1,J2]) + sage: x = sum(i*J.gens()[i] for i in range(len(J.gens()))) + sage: x.cartesian_factors() + (a1 + 2*a2, 3*b0 + 4*b1 + 5*b2 + 6*b3) """ - return self._factors + # Copy/pasted from CombinatorialFreeModule_CartesianProduct, + # but returning a tuple instead of a list. + return tuple(range(len(self.cartesian_factors()))) - def projections(self): + def cartesian_factors(self): + # Copy/pasted from CombinatorialFreeModule_CartesianProduct. + return self._sets + + def cartesian_factor(self, i): r""" - Return a pair of projections onto this algebra's factors. + Return the ``i``th factor of this algebra. + """ + return self._sets[i] + def _repr_(self): + # Copy/pasted from CombinatorialFreeModule_CartesianProduct. + from sage.categories.cartesian_product import cartesian_product + return cartesian_product.symbol.join("%s" % factor + for factor in self._sets) + + + @cached_method + def cartesian_projection(self, i): + r""" SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, - ....: ComplexHermitianEJA, - ....: DirectSumEJA) + sage: from mjo.eja.eja_algebra import (random_eja, + ....: JordanSpinEJA, + ....: HadamardEJA, + ....: RealSymmetricEJA, + ....: ComplexHermitianEJA) - EXAMPLES:: + EXAMPLES: + + The projection morphisms are Euclidean Jordan algebra + operators:: + + sage: J1 = HadamardEJA(2) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.cartesian_projection(0) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [1 0 0 0 0] + [0 1 0 0 0] + Domain: Euclidean Jordan algebra of dimension 2 over Algebraic + Real Field (+) Euclidean Jordan algebra of dimension 3 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic + Real Field + sage: J.cartesian_projection(1) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [0 0 1 0 0] + [0 0 0 1 0] + [0 0 0 0 1] + Domain: Euclidean Jordan algebra of dimension 2 over Algebraic + Real Field (+) Euclidean Jordan algebra of dimension 3 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic + Real Field + + The projections work the way you'd expect on the vector + representation of an element:: sage: J1 = JordanSpinEJA(2) sage: J2 = ComplexHermitianEJA(2) - sage: J = DirectSumEJA(J1,J2) - sage: (pi_left, pi_right) = J.projections() - sage: J.one().to_vector() - (1, 0, 1, 0, 0, 1) + sage: J = cartesian_product([J1,J2]) + sage: pi_left = J.cartesian_projection(0) + sage: pi_right = J.cartesian_projection(1) sage: pi_left(J.one()).to_vector() (1, 0) sage: pi_right(J.one()).to_vector() (1, 0, 0, 1) + sage: J.one().to_vector() + (1, 0, 1, 0, 0, 1) + + TESTS: + + The answer never changes:: + + sage: J1 = random_eja() + sage: J2 = random_eja() + sage: J = cartesian_product([J1,J2]) + sage: P0 = J.cartesian_projection(0) + sage: P1 = J.cartesian_projection(0) + sage: P0 == P1 + True """ - (J1,J2) = self.factors() - m = J1.dimension() - n = J2.dimension() - V_basis = self.vector_space().basis() - # Need to specify the dimensions explicitly so that we don't - # wind up with a zero-by-zero matrix when we want e.g. a - # zero-by-two matrix (important for composing things). - P1 = matrix(self.base_ring(), m, m+n, V_basis[:m]) - P2 = matrix(self.base_ring(), n, m+n, V_basis[m:]) - pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1) - pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2) - return (pi_left, pi_right) - - def inclusions(self): - r""" - Return the pair of inclusion maps from our factors into us. + offset = sum( self.cartesian_factor(k).dimension() + for k in range(i) ) + Ji = self.cartesian_factor(i) + Pi = self._module_morphism(lambda j: Ji.monomial(j - offset), + codomain=Ji) + + return EJAOperator(self,Ji,Pi.matrix()) + @cached_method + def cartesian_embedding(self, i): + r""" SETUP:: sage: from mjo.eja.eja_algebra import (random_eja, ....: JordanSpinEJA, - ....: RealSymmetricEJA, - ....: DirectSumEJA) + ....: HadamardEJA, + ....: RealSymmetricEJA) - EXAMPLES:: + EXAMPLES: + + The embedding morphisms are Euclidean Jordan algebra + operators:: + + sage: J1 = HadamardEJA(2) + sage: J2 = RealSymmetricEJA(2) + sage: J = cartesian_product([J1,J2]) + sage: J.cartesian_embedding(0) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [1 0] + [0 1] + [0 0] + [0 0] + [0 0] + Domain: Euclidean Jordan algebra of dimension 2 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 2 over + Algebraic Real Field (+) Euclidean Jordan algebra of + dimension 3 over Algebraic Real Field + sage: J.cartesian_embedding(1) + Linear operator between finite-dimensional Euclidean Jordan + algebras represented by the matrix: + [0 0 0] + [0 0 0] + [1 0 0] + [0 1 0] + [0 0 1] + Domain: Euclidean Jordan algebra of dimension 3 over + Algebraic Real Field + Codomain: Euclidean Jordan algebra of dimension 2 over + Algebraic Real Field (+) Euclidean Jordan algebra of + dimension 3 over Algebraic Real Field + + The embeddings work the way you'd expect on the vector + representation of an element:: sage: J1 = JordanSpinEJA(3) sage: J2 = RealSymmetricEJA(2) - sage: J = DirectSumEJA(J1,J2) - sage: (iota_left, iota_right) = J.inclusions() + sage: J = cartesian_product([J1,J2]) + sage: iota_left = J.cartesian_embedding(0) + sage: iota_right = J.cartesian_embedding(1) sage: iota_left(J1.zero()) == J.zero() True sage: iota_right(J2.zero()) == J.zero() @@ -2738,16 +3377,27 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): TESTS: + The answer never changes:: + + sage: J1 = random_eja() + sage: J2 = random_eja() + sage: J = cartesian_product([J1,J2]) + sage: E0 = J.cartesian_embedding(0) + sage: E1 = J.cartesian_embedding(0) + sage: E0 == E1 + True + Composing a projection with the corresponding inclusion should produce the identity map, and mismatching them should produce the zero map:: - sage: set_random_seed() sage: J1 = random_eja() sage: J2 = random_eja() - sage: J = DirectSumEJA(J1,J2) - sage: (iota_left, iota_right) = J.inclusions() - sage: (pi_left, pi_right) = J.projections() + sage: J = cartesian_product([J1,J2]) + sage: iota_left = J.cartesian_embedding(0) + sage: iota_right = J.cartesian_embedding(1) + sage: pi_left = J.cartesian_projection(0) + sage: pi_right = J.cartesian_projection(1) sage: pi_left*iota_left == J1.one().operator() True sage: pi_right*iota_right == J2.one().operator() @@ -2758,58 +3408,374 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra): True """ - (J1,J2) = self.factors() - m = J1.dimension() - n = J2.dimension() - V_basis = self.vector_space().basis() - # Need to specify the dimensions explicitly so that we don't - # wind up with a zero-by-zero matrix when we want e.g. a - # two-by-zero matrix (important for composing things). - I1 = matrix.column(self.base_ring(), m, m+n, V_basis[:m]) - I2 = matrix.column(self.base_ring(), n, m+n, V_basis[m:]) - iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1) - iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2) - return (iota_left, iota_right) + offset = sum( self.cartesian_factor(k).dimension() + for k in range(i) ) + Ji = self.cartesian_factor(i) + Ei = Ji._module_morphism(lambda j: self.monomial(j + offset), + codomain=self) + return EJAOperator(Ji,self,Ei.matrix()) - def inner_product(self, x, y): + + def subalgebra(self, basis, **kwargs): r""" - The standard Cartesian inner-product. + Create a subalgebra of this algebra from the given basis. - We project ``x`` and ``y`` onto our factors, and add up the - inner-products from the subalgebras. + Only overridden to allow us to use a special Cartesian product + subalgebra class. SETUP:: - sage: from mjo.eja.eja_algebra import (HadamardEJA, - ....: QuaternionHermitianEJA, - ....: DirectSumEJA) - - EXAMPLE:: - - sage: J1 = HadamardEJA(3,QQ) - sage: J2 = QuaternionHermitianEJA(2,QQ,orthonormalize=False) - sage: J = DirectSumEJA(J1,J2) - sage: x1 = J1.one() - sage: x2 = x1 - sage: y1 = J2.one() - sage: y2 = y1 - sage: x1.inner_product(x2) - 3 - sage: y1.inner_product(y2) - 2 - sage: J.one().inner_product(J.one()) - 5 + ....: QuaternionHermitianEJA) + + EXAMPLES: + + Subalgebras of Cartesian product EJAs have a different class + than those of non-Cartesian-product EJAs:: + + sage: J1 = HadamardEJA(2,field=QQ,orthonormalize=False) + sage: J2 = QuaternionHermitianEJA(0,field=QQ,orthonormalize=False) + sage: J = cartesian_product([J1,J2]) + sage: K1 = J1.subalgebra((J1.one(),), orthonormalize=False) + sage: K = J.subalgebra((J.one(),), orthonormalize=False) + sage: K1.__class__ is K.__class__ + False + + """ + from mjo.eja.eja_subalgebra import CartesianProductEJASubalgebra + return CartesianProductEJASubalgebra(self, basis, **kwargs) + +EJA.CartesianProduct = CartesianProductEJA + +class RationalBasisCartesianProductEJA(CartesianProductEJA, + RationalBasisEJA): + r""" + A separate class for products of algebras for which we know a + rational basis. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (EJA, + ....: HadamardEJA, + ....: JordanSpinEJA, + ....: RealSymmetricEJA) + + EXAMPLES: + + This gives us fast characteristic polynomial computations in + product algebras, too:: + + + sage: J1 = JordanSpinEJA(2) + sage: J2 = RealSymmetricEJA(3) + sage: J = cartesian_product([J1,J2]) + sage: J.characteristic_polynomial_of().degree() + 5 + sage: J.rank() + 5 + + TESTS: + + The ``cartesian_product()`` function only uses the first factor to + decide where the result will live; thus we have to be careful to + check that all factors do indeed have a ``rational_algebra()`` method + before we construct an algebra that claims to have a rational basis:: + + sage: J1 = HadamardEJA(2) + sage: jp = lambda X,Y: X*Y + sage: ip = lambda X,Y: X[0,0]*Y[0,0] + sage: b1 = matrix(QQ, [[1]]) + sage: J2 = EJA((b1,), jp, ip) + sage: cartesian_product([J2,J1]) # factor one not RationalBasisEJA + Euclidean Jordan algebra of dimension 1 over Algebraic Real + Field (+) Euclidean Jordan algebra of dimension 2 over Algebraic + Real Field + sage: cartesian_product([J1,J2]) # factor one is RationalBasisEJA + Traceback (most recent call last): + ... + ValueError: factor not a RationalBasisEJA + + """ + def __init__(self, algebras, **kwargs): + if not all( hasattr(r, "rational_algebra") for r in algebras ): + raise ValueError("factor not a RationalBasisEJA") + + CartesianProductEJA.__init__(self, algebras, **kwargs) + + @cached_method + def rational_algebra(self): + if self.base_ring() is QQ: + return self + + return cartesian_product([ + r.rational_algebra() for r in self.cartesian_factors() + ]) + + +RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA + +def random_eja(max_dimension=None, *args, **kwargs): + r""" + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS:: + + sage: n = ZZ.random_element(1,5) + sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False) + sage: J.dimension() <= n + True + + """ + # Use the ConcreteEJA default as the total upper bound (regardless + # of any whether or not any individual factors set a lower limit). + if max_dimension is None: + max_dimension = ConcreteEJA._max_random_instance_dimension() + J1 = ConcreteEJA.random_instance(max_dimension, *args, **kwargs) + + + # Roll the dice to see if we attempt a Cartesian product. + dice_roll = ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) + new_max_dimension = max_dimension - J1.dimension() + if new_max_dimension == 0 or dice_roll != 0: + # If it's already as big as we're willing to tolerate, just + # return it and don't worry about Cartesian products. + return J1 + else: + # Use random_eja() again so we can get more than two factors + # if the sub-call also Decides on a cartesian product. + J2 = random_eja(new_max_dimension, *args, **kwargs) + return cartesian_product([J1,J2]) + + +class ComplexSkewSymmetricEJA(RationalBasisEJA, ConcreteEJA): + r""" + The skew-symmetric EJA of order `2n` described in Faraut and + Koranyi's Exercise III.1.b. It has dimension `2n^2 - n`. + + It is (not obviously) isomorphic to the QuaternionHermitianEJA of + order `n`, as can be inferred by comparing rank/dimension or + explicitly from their "characteristic polynomial of" functions, + which just so happen to align nicely. + + SETUP:: + + sage: from mjo.eja.eja_algebra import (ComplexSkewSymmetricEJA, + ....: QuaternionHermitianEJA) + sage: from mjo.eja.eja_operator import EJAOperator + + EXAMPLES: + + This EJA is isomorphic to the quaternions:: + + sage: J = ComplexSkewSymmetricEJA(2, field=QQ, orthonormalize=False) + sage: K = QuaternionHermitianEJA(2, field=QQ, orthonormalize=False) + sage: jordan_isom_matrix = matrix.diagonal(QQ,[-1,1,1,1,1,-1]) + sage: phi = EJAOperator(J,K,jordan_isom_matrix) + sage: all( phi(x*y) == phi(x)*phi(y) + ....: for x in J.gens() + ....: for y in J.gens() ) + True + sage: x,y = J.random_elements(2) + sage: phi(x*y) == phi(x)*phi(y) + True + + TESTS: + + Random elements should satisfy the same conditions that the basis + elements do:: + + sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ, + ....: orthonormalize=False) + sage: x,y = K.random_elements(2) + sage: z = x*y + sage: x = x.to_matrix() + sage: y = y.to_matrix() + sage: z = z.to_matrix() + sage: all( e.is_skew_symmetric() for e in (x,y,z) ) + True + sage: J = -K.one().to_matrix() + sage: all( e*J == J*e.conjugate() for e in (x,y,z) ) + True + + The power law in Faraut & Koranyi's II.7.a is satisfied. + We're in a subalgebra of theirs, but powers are still + defined the same:: + + sage: K = ComplexSkewSymmetricEJA.random_instance(field=QQ, + ....: orthonormalize=False) + sage: x = K.random_element() + sage: k = ZZ.random_element(5) + sage: actual = x^k + sage: J = -K.one().to_matrix() + sage: expected = K(-J*(J*x.to_matrix())^k) + sage: actual == expected + True + + """ + @staticmethod + def _max_random_instance_size(max_dimension): + # Obtained by solving d = 2n^2 - n, which comes from noticing + # that, in 2x2 block form, any element of this algebra has a + # free skew-symmetric top-left block, a Hermitian top-right + # block, and two bottom blocks that are determined by the top. + # The ZZ-int-ZZ thing is just "floor." + return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4)) + + @classmethod + def random_instance(cls, max_dimension=None, *args, **kwargs): + """ + Return a random instance of this type of algebra. + """ + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) + return cls(n, **kwargs) + + @staticmethod + def _denormalized_basis(A): + """ + SETUP:: + + sage: from mjo.hurwitz import ComplexMatrixAlgebra + sage: from mjo.eja.eja_algebra import ComplexSkewSymmetricEJA + + TESTS: + + The basis elements are all skew-Hermitian:: + + sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension() + sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max) + sage: n = ZZ.random_element(n_max + 1) + sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ) + sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A) + sage: all( M.is_skew_symmetric() for M in B) + True + + The basis elements ``b`` all satisfy ``b*J == J*b.conjugate()``, + as in the definition of the algebra:: + + sage: d_max = ComplexSkewSymmetricEJA._max_random_instance_dimension() + sage: n_max = ComplexSkewSymmetricEJA._max_random_instance_size(d_max) + sage: n = ZZ.random_element(n_max + 1) + sage: A = ComplexMatrixAlgebra(2*n, scalars=QQ) + sage: I_n = matrix.identity(ZZ, n) + sage: J = matrix.block(ZZ, 2, 2, (0, I_n, -I_n, 0), subdivide=False) + sage: J = A.from_list(J.rows()) + sage: B = ComplexSkewSymmetricEJA._denormalized_basis(A) + sage: all( b*J == J*b.conjugate() for b in B ) + True """ - (pi_left, pi_right) = self.projections() - x1 = pi_left(x) - x2 = pi_right(x) - y1 = pi_left(y) - y2 = pi_right(y) + es = A.entry_algebra_gens() + gen = lambda A,m: A.monomial(m) + + basis = [] + + # The size of the blocks. We're going to treat these thing as + # 2x2 block matrices, + # + # [ x1 x2 ] + # [ -x2-conj x1-conj ] + # + # where x1 is skew-symmetric and x2 is Hermitian. + # + m = A.nrows()/2 + + # We only loop through the top half of the matrix, because the + # bottom can be constructed from the top. + for i in range(m): + # First do the top-left block, which is skew-symmetric. + # We can compute the bottom-right block in the process. + for j in range(i+1): + if i != j: + # Skew-symmetry implies zeros for (i == j). + for e in es: + # Top-left block's entry. + E_ij = gen(A, (i,j,e)) + E_ij -= gen(A, (j,i,e)) + + # Bottom-right block's entry. + F_ij = gen(A, (i+m,j+m,e)).conjugate() + F_ij -= gen(A, (j+m,i+m,e)).conjugate() + + basis.append(E_ij + F_ij) + + # Now do the top-right block, which is Hermitian, and compute + # the bottom-left block along the way. + for j in range(m,i+m+1): + if (i+m) == j: + # Hermitian matrices have real diagonal entries. + # Top-right block's entry. + E_ii = gen(A, (i,j,es[0])) + + # Bottom-left block's entry. Don't conjugate + # 'cause it's real. + E_ii -= gen(A, (i+m,j-m,es[0])) + basis.append(E_ii) + else: + for e in es: + # Top-right block's entry. BEWARE! We're not + # reflecting across the main diagonal as in + # (i,j)~(j,i). We're only reflecting across + # the diagonal for the top-right block. + E_ij = gen(A, (i,j,e)) + + # Shift it back to non-offset coords, transpose, + # conjugate, and put it back: + # + # (i,j) -> (i,j-m) -> (j-m, i) -> (j-m, i+m) + E_ij += gen(A, (j-m,i+m,e)).conjugate() + + # Bottom-left's block's below-diagonal entry. + # Just shift the top-right coords down m and + # left m. + F_ij = -gen(A, (i+m,j-m,e)).conjugate() + F_ij += -gen(A, (j,i,e)) # double-conjugate cancels - return (x1.inner_product(y1) + x2.inner_product(y2)) + basis.append(E_ij + F_ij) + return tuple( basis ) + @staticmethod + @cached_method + def _J_matrix(matrix_space): + n = matrix_space.nrows() // 2 + F = matrix_space.base_ring() + I_n = matrix.identity(F, n) + J = matrix.block(F, 2, 2, (0, I_n, -I_n, 0), subdivide=False) + return matrix_space.from_list(J.rows()) -random_eja = ConcreteEuclideanJordanAlgebra.random_instance + def J_matrix(self): + return ComplexSkewSymmetricEJA._J_matrix(self.matrix_space()) + + def __init__(self, n, field=AA, **kwargs): + # New code; always check the axioms. + #if "check_axioms" not in kwargs: kwargs["check_axioms"] = False + + from mjo.hurwitz import ComplexMatrixAlgebra + A = ComplexMatrixAlgebra(2*n, scalars=field) + J = ComplexSkewSymmetricEJA._J_matrix(A) + + def jordan_product(X,Y): + return (X*J*Y + Y*J*X)/2 + + def inner_product(X,Y): + return (X.conjugate_transpose()*Y).trace().real() + + super().__init__(self._denormalized_basis(A), + jordan_product, + inner_product, + field=field, + matrix_space=A, + **kwargs) + + # This algebra is conjectured (by me) to be isomorphic to + # the quaternion Hermitian EJA of size n, and the rank + # would follow from that. + #self.rank.set_cache(n) + self.one.set_cache( self(-J) )