deortho_vector_basis = tuple( V(b.list()) for b in basis )
from mjo.eja.eja_utils import gram_schmidt
- basis = gram_schmidt(basis, inner_product)
+ basis = tuple(gram_schmidt(basis, inner_product))
# Save the (possibly orthonormalized) matrix basis for
# later...
# Now we actually compute the multiplication and inner-product
# tables/matrices using the possibly-orthonormalized basis.
- self._inner_product_matrix = matrix.zero(field, n)
+ self._inner_product_matrix = matrix.identity(field, n)
self._multiplication_table = [ [0 for j in range(i+1)]
for i in range(n) ]
q_i = basis[i]
q_j = basis[j]
- elt = jordan_product(q_i, q_j)
- ip = inner_product(q_i, q_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 = W.coordinate_vector(V(elt.list()))
self._multiplication_table[i][j] = self.from_vector(elt)
- self._inner_product_matrix[i,j] = ip
- self._inner_product_matrix[j,i] = ip
+
+ 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()
This method should of course always return ``True``, unless
this algebra was constructed with ``check_axioms=False`` and
- passed an invalid multiplication table.
+ passed an invalid Jordan or inner-product.
"""
# Used to check whether or not something is zero in an inexact
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
+
+ The unit element in the Hadamard EJA is inherited in the
+ subalgebras generated by its elements::
sage: J = HadamardEJA(5)
sage: J.one()
e0 + e1 + e2 + e3 + e4
+ sage: x = sum(J.gens())
+ sage: A = x.subalgebra_generated_by(orthonormalize=False)
+ sage: A.one()
+ f0
+ sage: A.one().superalgebra_element()
+ e0 + e1 + e2 + e3 + e4
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()
+ 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::
+ ::
+
+ sage: set_random_seed()
+ 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,
+ regardless of the base field and whether or not we
+ orthonormalize::
sage: set_random_seed()
sage: J = random_eja()
sage: expected = matrix.identity(J.base_ring(), J.dimension())
sage: actual == expected
True
+ sage: x = J.random_element()
+ sage: A = x.subalgebra_generated_by()
+ sage: actual = A.one().operator().matrix()
+ sage: expected = matrix.identity(A.base_ring(), A.dimension())
+ sage: actual == expected
+ True
+
+ ::
+
+ sage: set_random_seed()
+ 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: J.one() == cached
True
+ ::
+
+ sage: set_random_seed()
+ 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.
jordan_product,
inner_product,
field=AA,
- orthonormalize=True,
check_field=True,
- check_axioms=True,
**kwargs):
if check_field:
if not all( all(b_i in QQ for b_i in b.list()) for b in basis ):
raise TypeError("basis not rational")
+ self._rational_algebra = None
if field is not QQ:
# There's no point in constructing the extra algebra if this
# one is already rational.
field=QQ,
orthonormalize=False,
check_field=False,
- check_axioms=False,
- **kwargs)
+ check_axioms=False)
super().__init__(basis,
jordan_product,
inner_product,
field=field,
check_field=check_field,
- check_axioms=check_axioms,
**kwargs)
@cached_method
a = ( a_i.change_ring(self.base_ring())
for a_i in self._rational_algebra._charpoly_coefficients() )
- # Now convert the coordinate variables back to the
+ if self._deortho_matrix is None:
+ # This can happen if our base ring was, say, AA and we
+ # chose not to (or didn't need to) orthonormalize. It's
+ # still faster to do the computations over QQ even if
+ # the numbers in the boxes stay the same.
+ return tuple(a)
+
+ # Otherwise, convert the coordinate variables back to the
# deorthonormalized ones.
R = self.coordinate_polynomial_ring()
from sage.modules.free_module_element import vector
class ComplexMatrixEJA(MatrixEJA):
+ # A manual dictionary-cache for the complex_extension() method,
+ # since apparently @classmethods can't also be @cached_methods.
+ _complex_extension = {}
+
+ @classmethod
+ def complex_extension(cls,field):
+ r"""
+ The complex field that we embed/unembed, as an extension
+ of the given ``field``.
+ """
+ if field in cls._complex_extension:
+ return cls._complex_extension[field]
+
+ # Sage doesn't know how to adjoin the complex "i" (the root of
+ # x^2 + 1) to a field in a general way. Here, we just enumerate
+ # all of the cases that I have cared to support so far.
+ 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
+ elif not field.is_exact():
+ # RDF or RR
+ F = field.complex_field()
+ else:
+ # Works for QQ and... maybe some other fields.
+ R = PolynomialRing(field, 'z')
+ z = R.gen()
+ F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+
+ cls._complex_extension[field] = F
+ return F
+
@staticmethod
def dimension_over_reals():
return 2
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]]))
+ a = z.real()
+ b = z.imag()
+ blocks.append(matrix(field, 2, [ [ a, b],
+ [-b, a] ]))
return matrix.block(field, n, blocks)
super(ComplexMatrixEJA,cls).real_unembed(M)
n = ZZ(M.nrows())
d = cls.dimension_over_reals()
-
- # 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()
-
- # Sage doesn't know how to adjoin the complex "i" (the root of
- # x^2 + 1) to a field in a general way. Here, we just enumerate
- # all of the cases that I have cared to support so far.
- 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
- elif not field.is_exact():
- # RDF or RR
- F = field.complex_field()
- else:
- # Works for QQ and... maybe some other fields.
- F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt())
+ F = cls.complex_extension(M.base_ring())
i = F.gen()
# Go top-left to bottom-right (reading order), converting every
sage: set_random_seed()
sage: n = ZZ.random_element(1,5)
- sage: field = QuadraticField(2, 'sqrt2')
sage: B = ComplexHermitianEJA._denormalized_basis(n)
sage: all( M.is_symmetric() for M in B)
True
# * The diagonal will (as a result) be real.
#
S = []
+ Eij = matrix.zero(F,n)
for i in range(n):
for j in range(i+1):
- Eij = matrix(F, n, lambda k,l: k==i and l==j)
+ # "build" E_ij
+ Eij[i,j] = 1
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())
+ Eij[j,i] = 1 # Eij = Eij + Eij.transpose()
+ Sij_real = cls.real_embed(Eij)
S.append(Sij_real)
- Sij_imag = cls.real_embed(I*Eij - I*Eij.transpose())
+ # Eij = I*Eij - I*Eij.transpose()
+ Eij[i,j] = I
+ Eij[j,i] = -I
+ Sij_imag = cls.real_embed(Eij)
S.append(Sij_imag)
+ Eij[j,i] = 0
+ # "erase" E_ij
+ Eij[i,j] = 0
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the complex extension "F".
return cls(n, **kwargs)
class QuaternionMatrixEJA(MatrixEJA):
+
+ # A manual dictionary-cache for the quaternion_extension() method,
+ # since apparently @classmethods can't also be @cached_methods.
+ _quaternion_extension = {}
+
+ @classmethod
+ def quaternion_extension(cls,field):
+ r"""
+ The quaternion field that we embed/unembed, as an extension
+ of the given ``field``.
+ """
+ if field in cls._quaternion_extension:
+ return cls._quaternion_extension[field]
+
+ Q = QuaternionAlgebra(field,-1,-1)
+
+ cls._quaternion_extension[field] = Q
+ return Q
+
@staticmethod
def dimension_over_reals():
return 4
# 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)
+ Q = cls.quaternion_extension(M.base_ring())
i,j,k = Q.gens()
# Go top-left to bottom-right (reading order), converting every
# * The diagonal will (as a result) be real.
#
S = []
+ Eij = matrix.zero(Q,n)
for i in range(n):
for j in range(i+1):
- Eij = matrix(Q, n, lambda k,l: k==i and l==j)
+ # "build" E_ij
+ Eij[i,j] = 1
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())
+ # Eij = Eij + Eij.transpose()
+ Eij[j,i] = 1
+ Sij_real = cls.real_embed(Eij)
S.append(Sij_real)
- Sij_I = cls.real_embed(I*Eij - I*Eij.transpose())
+ # Eij = I*(Eij - Eij.transpose())
+ Eij[i,j] = I
+ Eij[j,i] = -I
+ Sij_I = cls.real_embed(Eij)
S.append(Sij_I)
- Sij_J = cls.real_embed(J*Eij - J*Eij.transpose())
+ # Eij = J*(Eij - Eij.transpose())
+ Eij[i,j] = J
+ Eij[j,i] = -J
+ Sij_J = cls.real_embed(Eij)
S.append(Sij_J)
- Sij_K = cls.real_embed(K*Eij - K*Eij.transpose())
+ # Eij = K*(Eij - Eij.transpose())
+ Eij[i,j] = K
+ Eij[j,i] = -K
+ Sij_K = cls.real_embed(Eij)
S.append(Sij_K)
+ Eij[j,i] = 0
+ # "erase" E_ij
+ Eij[i,j] = 0
# Since we embedded these, we can drop back to the "field" that we
# started with instead of the quaternion algebra "Q".
....: for j in range(n-1) ]
sage: actual == expected
True
+
"""
def __init__(self, B, **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
def inner_product(x,y):
return (y.T*B*x)[0,0]
xbar = x[1:,0]
y0 = y[0,0]
ybar = y[1:,0]
- z0 = (y.T*x)[0,0]
+ z0 = inner_product(y,x)
zbar = y0*xbar + x0*ybar
- return P([0] + zbar.list())
-
- # 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
+ return P([z0] + zbar.list())
n = B.nrows()
column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )