]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: use new constructor for BilinearFormEJA.
[sage.d.git] / mjo / eja / eja_algebra.py
index 9d53bc5fdfd8fca177e43e49a187387fbe5a6085..9eba6699819956048747f648fd18bc80847fbbb6 100644 (file)
@@ -134,10 +134,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # 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.
-        self._multiplication_table = [
-            list(map(lambda x: self.from_vector(x), ls))
-            for ls in mult_table
-        ]
+        self._multiplication_table = [ [ self.vector_space().zero()
+                                         for i in range(n) ]
+                                       for j in range(n) ]
+        # take advantage of symmetry
+        for i in range(n):
+            for j in range(i+1):
+                elt = self.from_vector(mult_table[i][j])
+                self._multiplication_table[i][j] = elt
+                self._multiplication_table[j][i] = elt
 
         if check_axioms:
             if not self._is_commutative():
@@ -989,38 +994,13 @@ 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: 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: set_random_seed()    # long time
+            sage: J = random_eja()     # long time
+            sage: caches = J.rank()    # long time
+            sage: J.rank.clear_cache() # long time
+            sage: J.rank() == cached   # long time
+            True
 
-            sage: J = QuaternionHermitianEJA(2)
-            sage: J.rank.clear_cache()
-            sage: J.rank()
-            2
         """
         return len(self._charpoly_coefficients())
 
@@ -1045,6 +1025,103 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
+class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra):
+    def __init__(self,
+                 field,
+                 basis,
+                 jordan_product,
+                 inner_product,
+                 orthonormalize=True,
+                 prefix='e',
+                 category=None,
+                 check_field=True,
+                 check_axioms=True):
+
+        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
+                vector_basis = tuple( map(_mat2vec,basis) )
+                degree = basis[0].nrows()**2
+            else:
+                degree = basis[0].degree()
+
+        V = VectorSpace(field, degree)
+
+        # Compute this from "Q" (obtained from Gram-Schmidt) below as
+        # R = Q.solve_right(A), where the rows of "Q" are the
+        # orthonormalized vector_basis and and the rows of "A" are the
+        # original vector_basis.
+        self._deorthonormalization_matrix = None
+
+        if orthonormalize:
+            from mjo.eja.eja_utils import gram_schmidt
+            A = matrix(field, vector_basis)
+            vector_basis = gram_schmidt(vector_basis, inner_product)
+            W = V.span_of_basis( vector_basis )
+            Q = matrix(field, vector_basis)
+            # A = QR <==> A.T == R.T*Q.T
+            # So, Q.solve_right() is equivalent to the Q.T.solve_left()
+            # that we want.
+            self._deorthonormalization_matrix = Q.solve_right(A)
+
+            if basis_is_matrices:
+                from mjo.eja.eja_utils import _vec2mat
+                basis = tuple( map(_vec2mat,vector_basis) )
+
+        W = V.span_of_basis( vector_basis )
+
+        mult_table = [ [0 for i in range(n)] for j in range(n) ]
+        ip_table = [ [0 for i in range(n)] for j 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
+                mult_table[j][i] = elt
+                ip_table[i][j] = ip
+                ip_table[j][i] = ip
+
+        self._inner_product_matrix = matrix(field,ip_table)
+
+        if basis_is_matrices:
+            for m in basis:
+                m.set_immutable()
+        else:
+            basis = tuple( x.column() for x in basis )
+
+        super().__init__(field,
+                         mult_table,
+                         prefix,
+                         category,
+                         basis, # matrix basis
+                         check_field,
+                         check_axioms)
 
 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     r"""
@@ -1090,7 +1167,7 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
             return superclass._charpoly_coefficients()
 
         mult_table = tuple(
-            map(lambda x: x.to_vector(), ls)
+            tuple(map(lambda x: x.to_vector(), ls))
             for ls in self._multiplication_table
         )
 
@@ -2060,7 +2137,7 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         return cls(n, field, **kwargs)
 
 
-class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
+class HadamardEJA(RationalBasisEuclideanJordanAlgebraNg,
                   ConcreteEuclideanJordanAlgebra):
     """
     Return the Euclidean Jordan Algebra corresponding to the set
@@ -2103,22 +2180,17 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
     """
     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) ]
+        basis = V.basis()
 
-        # Inner products are real numbers and not algebra
-        # elements, so once we turn the algebra element
-        # into a vector in inner_product(), we never go
-        # back. As a result -- contrary to what we do with
-        # self._multiplication_table -- we store the inner
-        # product table as a plain old matrix and not as
-        # an algebra operator.
-        ip_table = matrix.identity(field,n)
-        self._inner_product_matrix = ip_table
+        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,
-                                          mult_table,
-                                          check_axioms=False,
+                                          basis,
+                                          jordan_product,
+                                          inner_product,
                                           **kwargs)
         self.rank.set_cache(n)
 
@@ -2143,7 +2215,7 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
         return cls(n, field, **kwargs)
 
 
-class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
+class BilinearFormEJA(RationalBasisEuclideanJordanAlgebraNg,
                       ConcreteEuclideanJordanAlgebra):
     r"""
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
@@ -2224,39 +2296,28 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
         True
     """
     def __init__(self, B, field=AA, **kwargs):
-        n = B.nrows()
-
         if not B.is_positive_definite():
             raise ValueError("bilinear form is not positive-definite")
 
+        n = B.nrows()
         V = VectorSpace(field, n)
-        mult_table = [[V.zero() for j in range(n)] for i in range(n)]
-        for i in range(n):
-            for j in range(n):
-                x = V.gen(i)
-                y = V.gen(j)
-                x0 = x[0]
-                xbar = x[1:]
-                y0 = y[0]
-                ybar = y[1:]
-                z0 = (B*x).inner_product(y)
-                zbar = y0*xbar + x0*ybar
-                z = V([z0] + zbar.list())
-                mult_table[i][j] = z
-
-        # Inner products are real numbers and not algebra
-        # elements, so once we turn the algebra element
-        # into a vector in inner_product(), we never go
-        # back. As a result -- contrary to what we do with
-        # self._multiplication_table -- we store the inner
-        # product table as a plain old matrix and not as
-        # an algebra operator.
-        ip_table = B
-        self._inner_product_matrix = ip_table
+
+        def inner_product(x,y):
+            return (B*x).inner_product(y)
+
+        def jordan_product(x,y):
+            x0 = x[0]
+            xbar = x[1:]
+            y0 = y[0]
+            ybar = y[1:]
+            z0 = inner_product(x,y)
+            zbar = y0*xbar + x0*ybar
+            return V([z0] + zbar.list())
 
         super(BilinearFormEJA, self).__init__(field,
-                                              mult_table,
-                                              check_axioms=False,
+                                              V.basis(),
+                                              jordan_product,
+                                              inner_product,
                                               **kwargs)
 
         # The rank of this algebra is two, unless we're in a