X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=2b769ac447bce5b149782bbf96ada23631f9b2c9;hb=823fb2a587e26436f46854fe44be0e8df46a6715;hp=4b8d466c3456f9ec1daf6ca2d701f3078a3524d1;hpb=ef4a0674cdc17613ce65103fc679506e7c7008de;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 4b8d466..2b769ac 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -17,7 +17,7 @@ from sage.misc.lazy_import import lazy_import from sage.misc.prandom import choice from sage.misc.table import table from sage.modules.free_module import FreeModule, VectorSpace -from sage.rings.all import (ZZ, QQ, RR, RLF, CLF, +from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF, PolynomialRing, QuadraticField) from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement @@ -54,7 +54,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): def __init__(self, field, mult_table, - rank, prefix='e', category=None, natural_basis=None, @@ -91,7 +90,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): # a real embedding. raise ValueError('field is not real') - self._rank = rank self._natural_basis = natural_basis if category is None: @@ -194,6 +192,24 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): coords = W.coordinate_vector(_mat2vec(elt)) return self.from_vector(coords) + @staticmethod + def _max_test_case_size(): + """ + 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 quite vague -- 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. + + We default to five in this class, which is safe in `R^n`. The + matrix algebra subclasses (or any class where the "size" is + interpreted to be far less than the dimension) should override + with a smaller number. + """ + return 5 def _repr_(self): """ @@ -207,8 +223,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Ensure that it says what we think it says:: - sage: JordanSpinEJA(2, field=QQ) - Euclidean Jordan algebra of dimension 2 over Rational Field + 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 @@ -267,8 +283,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return V.span_of_basis(b) - - @cached_method def _charpoly_coeff(self, i): """ Return the coefficient polynomial "a_{i}" of this algebra's @@ -278,102 +292,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): store the trace/determinant (a_{r-1} and a_{0} respectively) separate from the entire characteristic polynomial. """ - (A_of_x, x, xr, detA) = self._charpoly_matrix_system() - R = A_of_x.base_ring() - - if i == self.rank(): - return R.one() - if i > self.rank(): - # Guaranteed by theory - return R.zero() - - # Danger: the in-place modification is done for performance - # reasons (reconstructing a matrix with huge polynomial - # entries is slow), but I don't know how cached_method works, - # so it's highly possible that we're modifying some global - # list variable by reference, here. In other words, you - # probably shouldn't call this method twice on the same - # algebra, at the same time, in two threads - Ai_orig = A_of_x.column(i) - A_of_x.set_column(i,xr) - numerator = A_of_x.det() - A_of_x.set_column(i,Ai_orig) - - # We're relying on the theory here to ensure that each a_i is - # indeed back in R, and the added negative signs are to make - # the whole charpoly expression sum to zero. - return R(-numerator/detA) - - - @cached_method - def _charpoly_matrix_system(self): - """ - Compute the matrix whose entries A_ij are polynomials in - X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector - corresponding to `x^r` and the determinent of the matrix A = - [A_ij]. In other words, all of the fixed (cachable) data needed - to compute the coefficients of the characteristic polynomial. - """ - r = self.rank() - n = self.dimension() - - # Turn my vector space into a module so that "vectors" can - # have multivatiate polynomial entries. - names = tuple('X' + str(i) for i in range(1,n+1)) - R = PolynomialRing(self.base_ring(), names) - - # Using change_ring() on the parent's vector space doesn't work - # here because, in a subalgebra, that vector space has a basis - # and change_ring() tries to bring the basis along with it. And - # that doesn't work unless the new ring is a PID, which it usually - # won't be. - V = FreeModule(R,n) - - # Now let x = (X1,X2,...,Xn) be the vector whose entries are - # indeterminates... - x = V(names) - - # And figure out the "left multiplication by x" matrix in - # that setting. - lmbx_cols = [] - monomial_matrices = [ self.monomial(i).operator().matrix() - for i in range(n) ] # don't recompute these! - for k in range(n): - ek = self.monomial(k).to_vector() - lmbx_cols.append( - sum( x[i]*(monomial_matrices[i]*ek) - for i in range(n) ) ) - Lx = matrix.column(R, lmbx_cols) - - # Now we can compute powers of x "symbolically" - x_powers = [self.one().to_vector(), x] - for d in range(2, r+1): - x_powers.append( Lx*(x_powers[-1]) ) - - idmat = matrix.identity(R, n) - - W = self._charpoly_basis_space() - W = W.change_ring(R.fraction_field()) - - # Starting with the standard coordinates x = (X1,X2,...,Xn) - # and then converting the entries to W-coordinates allows us - # to pass in the standard coordinates to the charpoly and get - # back the right answer. Specifically, with x = (X1,X2,...,Xn), - # we have - # - # W.coordinates(x^2) eval'd at (standard z-coords) - # = - # W-coords of (z^2) - # = - # W-coords of (standard coords of x^2 eval'd at std-coords of z) - # - # We want the middle equivalent thing in our matrix, but use - # the first equivalent thing instead so that we can pass in - # standard coordinates. - x_powers = [ W.coordinate_vector(xp) for xp in x_powers ] - l2 = [idmat.column(k-1) for k in range(r+1, n+1)] - A_of_x = matrix.column(R, n, (x_powers[:r] + l2)) - return (A_of_x, x, x_powers[r], A_of_x.det()) + return self._charpoly_coefficients()[i] @cached_method @@ -420,20 +339,21 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): r = self.rank() n = self.dimension() - # The list of coefficient polynomials a_0, a_1, a_2, ..., a_n. - a = [ self._charpoly_coeff(i) for i in range(r+1) ] + # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1). + a = self._charpoly_coefficients() # We go to a bit of trouble here to reorder the # indeterminates, so that it's easier to evaluate the # characteristic polynomial at x's coordinates and get back # something in terms of t, which is what we want. - R = a[0].parent() S = PolynomialRing(self.base_ring(),'t') t = S.gen(0) - S = PolynomialRing(S, R.variable_names()) - t = S(t) + if r > 0: + R = a[0].parent() + S = PolynomialRing(S, R.variable_names()) + t = S(t) - return sum( a[k]*(t**k) for k in range(len(a)) ) + return (t**r + sum( a[k]*(t**k) for k in range(r) )) def inner_product(self, x, y): @@ -551,8 +471,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Finite family {0: e0, 1: e1, 2: e2} sage: J.natural_basis() ( - [1 0] [ 0 1/2*sqrt2] [0 0] - [0 0], [1/2*sqrt2 0], [0 1] + [1 0] [ 0 0.7071067811865475?] [0 0] + [0 0], [0.7071067811865475? 0], [0 1] ) :: @@ -757,7 +677,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half J1 = trivial # eigenvalue one - for (eigval, eigspace) in c.operator().matrix().left_eigenspaces(): + for (eigval, eigspace) in c.operator().matrix().right_eigenspaces(): if eigval == ~(self.base_ring()(2)): J5 = eigspace else: @@ -773,65 +693,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): return (J0, J5, J1) - def orthogonal_idempotents(self): - r""" - Generate a set of `r` orthogonal idempotents for this algebra, - where `r` is its rank. - - This implementation is based on the so-called "central - orthogonal idempotents" implemented for (semisimple) centers - of SageMath ``FiniteDimensionalAlgebrasWithBasis``. Since all - Euclidean Jordan algebas are commutative (and thus equal to - their own centers) and semisimple, the method should work more - or less as implemented, if it ever worked in the first place. - (I don't know the justification for the original implementation. - yet). - - How it works: we loop through the algebras generators, looking - for their eigenspaces. If there's more than one eigenspace, - and if they result in more than one subalgebra, then we split - those subalgebras recursively until we get to subalgebras of - dimension one (whose idempotent is the unit element). Why does - some generator have to produce at least two subalgebras? I - dunno. But it seems to work. - - Beware that Koecher defines the "center" of a Jordan algebra to - be something else, because the usual definition is stupid in a - (necessarily commutative) Jordan algebra. - """ - if self.dimension() == 1: - return [self.one()] - - for g in self.gens(): - eigenpairs = g.operator().matrix().right_eigenspaces() - if len(eigenpairs) >= 2: - subalgebras = [] - for eigval, eigspace in eigenpairs: - # Make sub-EJAs from the matrix eigenspaces... - sb = tuple( self.from_vector(b) for b in eigspace.basis() ) - try: - # This will fail if e.g. the eigenspace basis - # contains two elements and their product - # isn't a linear combination of the two of - # them (i.e. the generated EJA isn't actually - # two dimensional). - s = FiniteDimensionalEuclideanJordanSubalgebra(self, sb) - subalgebras.append(s) - except: - pass - if len(subalgebras) >= 2: - # apply this method recursively. - return tuple( c.superalgebra_element() - for subalgebra in subalgebras - for c in subalgebra.orthogonal_idempotents() ) - - # If we got here, the algebra didn't decompose, at least not when we looked at - # the eigenspaces corresponding only to basis elements of the algebra. The - # implementation I stole says that this should work because of Schur's Lemma, - # so I personally blame Schur's Lemma if it does not. - raise Exception("Schur's Lemma didn't work!") - - def random_elements(self, count): """ Return ``count`` random elements as a tuple. @@ -850,23 +711,75 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): True """ - return tuple( self.random_element() for idx in range(count) ) + return tuple( self.random_element() for idx in range(count) ) + @classmethod + def random_instance(cls, field=AA, **kwargs): + """ + Return a random instance of this type of algebra. - def rank(self): + Beware, this will crash for "most instances" because the + constructor below looks wrong. """ - Return the rank of this EJA. + if cls is TrivialEJA: + # The TrivialEJA class doesn't take an "n" argument because + # there's only one. + return cls(field) - ALGORITHM: + n = ZZ.random_element(cls._max_test_case_size()) + 1 + return cls(n, field, **kwargs) - The author knows of no algorithm to compute the rank of an EJA - where only the multiplication table is known. In lieu of one, we - require the rank to be specified when the algebra is created, - and simply pass along that number here. + @cached_method + def _charpoly_coefficients(self): + r""" + The `r` polynomial coefficients of the "characteristic polynomial + of" function. + """ + n = self.dimension() + var_names = [ "X" + str(z) for z in range(1,n+1) ] + R = PolynomialRing(self.base_ring(), var_names) + 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) + # Compute an extra power in case the rank is equal to + # the dimension (otherwise, we would stop at x^(r-1)). + x_powers = [ (L_x**k)*self.one().to_vector() + for k in range(n+1) ] + A = matrix.column(F, x_powers[:n]) + AE = A.extended_echelon_form() + E = AE[:,n:] + A_rref = AE[:,:n] + r = A_rref.rank() + b = x_powers[r] + + # 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] + + @cached_method + def rank(self): + r""" + Return the rank of this EJA. + + This is a cached method because we know the rank a priori for + all of the algebras we can construct. Thus we can avoid the + expensive ``_charpoly_coefficients()`` call unless we truly + need to compute the whole characteristic polynomial. SETUP:: - sage: from mjo.eja.eja_algebra import (JordanSpinEJA, + sage: from mjo.eja.eja_algebra import (HadamardEJA, + ....: JordanSpinEJA, ....: RealSymmetricEJA, ....: ComplexHermitianEJA, ....: QuaternionHermitianEJA, @@ -907,8 +820,43 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): sage: r > 0 or (r == 0 and J.is_trivial()) True + Ensure that computing the rank actually works, since the ranks + of all simple algebras are known and will be cached by default:: + + sage: J = HadamardEJA(4) + sage: J.rank.clear_cache() + sage: J.rank() + 4 + + :: + + sage: J = JordanSpinEJA(4) + sage: J.rank.clear_cache() + sage: J.rank() + 2 + + :: + + sage: J = RealSymmetricEJA(3) + sage: J.rank.clear_cache() + sage: J.rank() + 3 + + :: + + sage: J = ComplexHermitianEJA(2) + sage: J.rank.clear_cache() + sage: J.rank() + 2 + + :: + + sage: J = QuaternionHermitianEJA(2) + sage: J.rank.clear_cache() + sage: J.rank() + 2 """ - return self._rank + return len(self._charpoly_coefficients()) def vector_space(self): @@ -932,61 +880,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule): Element = FiniteDimensionalEuclideanJordanAlgebraElement -class KnownRankEJA(object): - """ - A class for algebras that we actually know we can construct. The - main issue is that, for most of our methods to make sense, we need - to know the rank of our algebra. Thus we can't simply generate a - "random" algebra, or even check that a given basis and product - satisfy the axioms; because even if everything looks OK, we wouldn't - know the rank we need to actuallty build the thing. - - Not really a subclass of FDEJA because doing that causes method - resolution errors, e.g. - - TypeError: Error when calling the metaclass bases - Cannot create a consistent method resolution - order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra, - KnownRankEJA - - """ - @staticmethod - def _max_test_case_size(): - """ - 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 quite vague -- 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. - - We default to five in this class, which is safe in `R^n`. The - matrix algebra subclasses (or any class where the "size" is - interpreted to be far less than the dimension) should override - with a smaller number. - """ - return 5 - - @classmethod - def random_instance(cls, field=QQ, **kwargs): - """ - Return a random instance of this type of algebra. - - Beware, this will crash for "most instances" because the - constructor below looks wrong. - """ - if cls is TrivialEJA: - # The TrivialEJA class doesn't take an "n" argument because - # there's only one. - return cls(field) - - n = ZZ.random_element(cls._max_test_case_size()) + 1 - return cls(n, field, **kwargs) - - -class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): +class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra): """ Return the Euclidean Jordan Algebra corresponding to the set `R^n` under the Hadamard product. @@ -1026,13 +920,14 @@ class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): (r0, r1, r2) """ - def __init__(self, n, field=QQ, **kwargs): + def __init__(self, n, field=AA, **kwargs): V = VectorSpace(field, n) mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ] for i in range(n) ] fdeja = super(HadamardEJA, self) - return fdeja.__init__(field, mult_table, rank=n, **kwargs) + fdeja.__init__(field, mult_table, **kwargs) + self.rank.set_cache(n) def inner_product(self, x, y): """ @@ -1059,7 +954,7 @@ class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): return x.to_vector().inner_product(y.to_vector()) -def random_eja(field=QQ, nontrivial=False): +def random_eja(field=AA, nontrivial=False): """ Return a "random" finite-dimensional Euclidean Jordan Algebra. @@ -1073,9 +968,13 @@ def random_eja(field=QQ, nontrivial=False): Euclidean Jordan algebra of dimension... """ - eja_classes = KnownRankEJA.__subclasses__() - if nontrivial: - eja_classes.remove(TrivialEJA) + eja_classes = [HadamardEJA, + JordanSpinEJA, + RealSymmetricEJA, + ComplexHermitianEJA, + QuaternionHermitianEJA] + if not nontrivial: + eja_classes.append(TrivialEJA) classname = choice(eja_classes) return classname.random_instance(field=field) @@ -1091,20 +990,20 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): # field can have dimension 4 (quaternions) too. return 2 - def __init__(self, field, basis, rank, normalize_basis=True, **kwargs): + def __init__(self, field, basis, normalize_basis=True, **kwargs): """ Compared to the superclass constructor, we take a basis instead of a multiplication table because the latter can be computed in terms of the former when the product is known (like it is here). """ - # Used in this class's fast _charpoly_coeff() override. + # Used in this class's fast _charpoly_coefficients() override. self._basis_normalizers = None # We're going to loop through this a few times, so now's a good # time to ensure that it isn't a generator expression. basis = tuple(basis) - if rank > 1 and normalize_basis: + if len(basis) > 1 and normalize_basis: # We'll need sqrt(2) to normalize the basis, and this # winds up in the multiplication table, so the whole # algebra needs to be over the field extension. @@ -1121,45 +1020,31 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra): Qs = self.multiplication_table_from_matrix_basis(basis) fdeja = super(MatrixEuclideanJordanAlgebra, self) - return fdeja.__init__(field, - Qs, - rank=rank, - natural_basis=basis, - **kwargs) + fdeja.__init__(field, Qs, natural_basis=basis, **kwargs) + return @cached_method - def _charpoly_coeff(self, i): - """ + def _charpoly_coefficients(self): + r""" Override the parent method with something that tries to compute over a faster (non-extension) field. """ if self._basis_normalizers is None: # We didn't normalize, so assume that the basis we started # with had entries in a nice field. - return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i) + return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients() else: basis = ( (b/n) for (b,n) in zip(self.natural_basis(), self._basis_normalizers) ) # Do this over the rationals and convert back at the end. + # Only works because we know the entries of the basis are + # integers. J = MatrixEuclideanJordanAlgebra(QQ, basis, - self.rank(), normalize_basis=False) - (_,x,_,_) = J._charpoly_matrix_system() - p = J._charpoly_coeff(i) - # p might be missing some vars, have to substitute "optionally" - pairs = zip(x.base_ring().gens(), self._basis_normalizers) - substitutions = { v: v*c for (v,c) in pairs } - result = p.subs(substitutions) - - # The result of "subs" can be either a coefficient-ring - # element or a polynomial. Gotta handle both cases. - if result in QQ: - return self.base_ring()(result) - else: - return result.change_ring(self.base_ring()) + return J._charpoly_coefficients() @staticmethod @@ -1256,7 +1141,7 @@ class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return M -class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): +class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner @@ -1279,8 +1164,8 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): In theory, our "field" can be any subfield of the reals:: - sage: RealSymmetricEJA(2, AA) - Euclidean Jordan algebra of dimension 3 over Algebraic Real Field + sage: RealSymmetricEJA(2, RDF) + Euclidean Jordan algebra of dimension 3 over Real Double Field sage: RealSymmetricEJA(2, RR) Euclidean Jordan algebra of dimension 3 over Real Field with 53 bits of precision @@ -1372,9 +1257,10 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA): return 4 # Dimension 10 - def __init__(self, n, field=QQ, **kwargs): + def __init__(self, n, field=AA, **kwargs): basis = self._denormalized_basis(n, field) - super(RealSymmetricEJA, self).__init__(field, basis, n, **kwargs) + super(RealSymmetricEJA, self).__init__(field, basis, **kwargs) + self.rank.set_cache(n) class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @@ -1392,7 +1278,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): EXAMPLES:: - sage: F = QuadraticField(-1, 'i') + sage: F = QuadraticField(-1, 'I') sage: x1 = F(4 - 2*i) sage: x2 = F(1 + 2*i) sage: x3 = F(-i) @@ -1412,7 +1298,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): sage: set_random_seed() sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size() sage: n = ZZ.random_element(n_max) - sage: F = QuadraticField(-1, 'i') + sage: F = QuadraticField(-1, 'I') sage: X = random_matrix(F, n) sage: Y = random_matrix(F, n) sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X) @@ -1455,15 +1341,15 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): ....: [ 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] + [ 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: F = QuadraticField(-1, 'I') sage: M = random_matrix(F, 3) sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M) sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M @@ -1481,7 +1367,12 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): field = M.base_ring() R = PolynomialRing(field, 'z') z = R.gen() - F = field.extension(z**2 + 1, 'i', embedding=CLF(-1).sqrt()) + 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 @@ -1522,7 +1413,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): sage: Ye = y.natural_representation() sage: X = ComplexHermitianEJA.real_unembed(Xe) sage: Y = ComplexHermitianEJA.real_unembed(Ye) - sage: expected = (X*Y).trace().vector()[0] + sage: expected = (X*Y).trace().real() sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye) sage: actual == expected True @@ -1531,7 +1422,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2 -class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): +class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, @@ -1546,8 +1437,8 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): In theory, our "field" can be any subfield of the reals:: - sage: ComplexHermitianEJA(2, AA) - Euclidean Jordan algebra of dimension 4 over Algebraic Real Field + 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 @@ -1657,9 +1548,10 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA): return ( s.change_ring(field) for s in S ) - def __init__(self, n, field=QQ, **kwargs): + def __init__(self, n, field=AA, **kwargs): basis = self._denormalized_basis(n,field) - super(ComplexHermitianEJA,self).__init__(field, basis, n, **kwargs) + super(ComplexHermitianEJA,self).__init__(field, basis, **kwargs) + self.rank.set_cache(n) class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): @@ -1709,7 +1601,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): if M.ncols() != n: raise ValueError("the matrix 'M' must be square") - F = QuadraticField(-1, 'i') + F = QuadraticField(-1, 'I') i = F.gen() blocks = [] @@ -1785,10 +1677,10 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): raise ValueError('bad on-diagonal submatrix') if submat[0,1] != -submat[1,0].conjugate(): raise ValueError('bad off-diagonal submatrix') - z = submat[0,0].vector()[0] # real part - z += submat[0,0].vector()[1]*i # imag part - z += submat[0,1].vector()[0]*j # real part - z += submat[0,1].vector()[1]*k # imag part + 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) @@ -1825,8 +1717,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra): return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4 -class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, - KnownRankEJA): +class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra): """ The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the @@ -1841,8 +1732,8 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, In theory, our "field" can be any subfield of the reals:: - sage: QuaternionHermitianEJA(2, AA) - Euclidean Jordan algebra of dimension 6 over Algebraic Real Field + sage: QuaternionHermitianEJA(2, RDF) + Euclidean Jordan algebra of dimension 6 over Real Double Field sage: QuaternionHermitianEJA(2, RR) Euclidean Jordan algebra of dimension 6 over Real Field with 53 bits of precision @@ -1953,12 +1844,13 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra, return ( s.change_ring(field) for s in S ) - def __init__(self, n, field=QQ, **kwargs): + def __init__(self, n, field=AA, **kwargs): basis = self._denormalized_basis(n,field) - super(QuaternionHermitianEJA,self).__init__(field, basis, n, **kwargs) + super(QuaternionHermitianEJA,self).__init__(field, basis, **kwargs) + self.rank.set_cache(n) -class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): +class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra): 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 = @@ -2013,7 +1905,7 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): sage: actual == expected True """ - def __init__(self, n, field=QQ, B=None, **kwargs): + def __init__(self, n, field=AA, B=None, **kwargs): if B is None: self._B = matrix.identity(field, max(0,n-1)) else: @@ -2038,7 +1930,8 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): # one-dimensional ambient space (because the rank is bounded # by the ambient dimension). fdeja = super(BilinearFormEJA, self) - return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs) + fdeja.__init__(field, mult_table, **kwargs) + self.rank.set_cache(min(n,2)) def inner_product(self, x, y): r""" @@ -2126,13 +2019,13 @@ class JordanSpinEJA(BilinearFormEJA): True """ - def __init__(self, n, field=QQ, **kwargs): + def __init__(self, n, field=AA, **kwargs): # This is a special case of the BilinearFormEJA with the identity # matrix as its bilinear form. return super(JordanSpinEJA, self).__init__(n, field, **kwargs) -class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): +class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra): """ The trivial Euclidean Jordan algebra consisting of only a zero element. @@ -2156,14 +2049,15 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA): sage: J.one().norm() 0 sage: J.one().subalgebra_generated_by() - Euclidean Jordan algebra of dimension 0 over Rational Field + Euclidean Jordan algebra of dimension 0 over Algebraic Real Field sage: J.rank() 0 """ - def __init__(self, field=QQ, **kwargs): + def __init__(self, field=AA, **kwargs): mult_table = [] fdeja = super(TrivialEJA, self) # The rank is zero using my definition, namely the dimension of the # largest subalgebra generated by any element. - return fdeja.__init__(field, mult_table, rank=0, **kwargs) + fdeja.__init__(field, mult_table, **kwargs) + self.rank.set_cache(0)