]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: eliminate the special element subalgebra class.
[sage.d.git] / mjo / eja / eja_algebra.py
index 081abd303c32d35c9376ad99ecf2732a14ea3e9b..83ec50e5a215417811fe09d72c61251081772160 100644 (file)
@@ -133,7 +133,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             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...
@@ -158,7 +158,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # 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) ]
 
@@ -171,15 +171,20 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 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()
@@ -315,7 +320,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         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
@@ -745,23 +750,57 @@ class FiniteDimensionalEJA(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
+
+        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()
@@ -769,6 +808,27 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             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::
@@ -780,6 +840,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             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.
@@ -1196,6 +1265,7 @@ class RationalBasisEJA(FiniteDimensionalEJA):
             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.
@@ -1256,7 +1326,14 @@ class RationalBasisEJA(FiniteDimensionalEJA):
         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
@@ -2307,10 +2384,21 @@ class BilinearFormEJA(ConcreteEJA):
         ....:              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]
@@ -2321,13 +2409,9 @@ class BilinearFormEJA(ConcreteEJA):
             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() )