eja: refactor all matrix classes upwards (note: everything broken).

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 3027075374a177ba5f72bbd690a2c07dad02d136..106a0cddec06355a952e71d75699677b25dd7da9 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -172,6 +172,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         category = MagmaticAlgebras(field).FiniteDimensional()
         category = category.WithBasis().Unital().Commutative()
 
         category = MagmaticAlgebras(field).FiniteDimensional()
         category = category.WithBasis().Unital().Commutative()
 

+        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()
         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()


@@ -1732,6 +1737,86 @@ class ConcreteEJA(FiniteDimensionalEJA):
 
 
 class MatrixEJA:
 
 
 class MatrixEJA:

+    @staticmethod
+    def _denormalized_basis(A):
+        """
+        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.hurwitz import (ComplexMatrixAlgebra,
+            ....:                          QuaternionMatrixAlgebra,
+            ....:                          OctonionMatrixAlgebra)
+            sage: from mjo.eja.eja_algebra import MatrixEJA
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = MatrixSpace(QQ, n)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in  B)
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = ComplexMatrixAlgebra(n, scalars=QQ)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in  B)
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = QuaternionMatrixAlgebra(n, scalars=QQ)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in B )
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: A = OctonionMatrixAlgebra(n, scalars=QQ)
+            sage: B = MatrixEJA._denormalized_basis(A)
+            sage: all( M.is_hermitian() for M in B )
+            True
+
+        """
+        # 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])
+
+        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)
+
+        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 )
+

     @staticmethod
     def jordan_product(X,Y):
         return (X*Y + Y*X)/2
     @staticmethod
     def jordan_product(X,Y):
         return (X*Y + Y*X)/2


@@ -1741,8 +1826,56 @@ class MatrixEJA:
         r"""
         A trace inner-product for matrices that aren't embedded in the
         reals. It takes MATRICES as arguments, not EJA elements.
         r"""
         A trace inner-product for matrices that aren't embedded in the
         reals. It takes MATRICES as arguments, not EJA elements.

+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  QuaternionHermitianEJA,
+            ....:                                  OctonionHermitianEJA)
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
+        ::
+
+            sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
+        ::
+
+            sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+
+        ::
+
+            sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: I = J.one().to_matrix()
+            sage: J.trace_inner_product(I, -I)
+            -2
+

         """
         """

-        return (X*Y).trace().real()
+        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()
+

 
 
 class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
 
 
 class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):


@@ -1810,38 +1943,6 @@ class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
         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,ZZ)
-            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
     @staticmethod
     def _max_random_instance_size():
         return 4 # Dimension 10


@@ -1859,23 +1960,18 @@ class RealSymmetricEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 

-        associative = False
-        if n <= 1:
-            associative = True
-
-        super().__init__(self._denormalized_basis(n,field),
+        A = MatrixSpace(field, n)
+        super().__init__(self._denormalized_basis(A),

                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,
                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,

-                         associative=associative,

                          **kwargs)
 
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
                          **kwargs)
 
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)

-        idV = self.matrix_space().one()
-        self.one.set_cache(self(idV))
+        self.one.set_cache(self(A.one()))

 
 
 
 
 
 


@@ -1936,93 +2032,24 @@ class ComplexHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
         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: B = ComplexHermitianEJA._denormalized_basis(n,QQ)
-            sage: all( M.is_hermitian() for M in  B)
-            True
-
-        """
-        from mjo.hurwitz import ComplexMatrixAlgebra
-        A = ComplexMatrixAlgebra(n, scalars=field)
-        es = A.entry_algebra_gens()
-
-        basis = []
-        for i in range(n):
-            for j in range(i+1):
-                if i == j:
-                    E_ii = A.monomial( (i,j,es[0]) )
-                    basis.append(E_ii)
-                else:
-                    for e in es:
-                        E_ij  = A.monomial( (i,j,e)             )
-                        ec = e.conjugate()
-                        # If the conjugate has a negative sign in front
-                        # of it, (j,i,ec) won't be a monomial!
-                        if (j,i,ec) in A.indices():
-                            E_ij += A.monomial( (j,i,ec) )
-                        else:
-                            E_ij -= A.monomial( (j,i,-ec) )
-                        basis.append(E_ij)
-
-        return tuple( basis )
-
-    @staticmethod
-    def trace_inner_product(X,Y):
-        r"""
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
-
-        TESTS::
-
-            sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
-            sage: I = J.one().to_matrix()
-            sage: J.trace_inner_product(I, -I)
-            -2
-
-        """
-        return (X*Y).trace().real()
-

     def __init__(self, n, field=AA, **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
 
     def __init__(self, n, field=AA, **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
 

-        associative = False
-        if n <= 1:
-            associative = True
-
-        super().__init__(self._denormalized_basis(n,field),
+        from mjo.hurwitz import ComplexMatrixAlgebra
+        A = ComplexMatrixAlgebra(n, scalars=field)
+        super().__init__(self._denormalized_basis(A),

                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,
                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,

-                         associative=associative,

                          **kwargs)
                          **kwargs)

+

         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)

-        idV = self.matrix_space().one()
-        self.one.set_cache(self(idV))
+        self.one.set_cache(self(A.one()))

 
     @staticmethod
     def _max_random_instance_size():
 
     @staticmethod
     def _max_random_instance_size():


@@ -2094,97 +2121,24 @@ class QuaternionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
 
     """

-    @classmethod
-    def _denormalized_basis(cls, n, field):
-        """
-        Returns a basis for the space of quaternion 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 QuaternionHermitianEJA
-
-        TESTS::
-
-            sage: set_random_seed()
-            sage: n = ZZ.random_element(1,5)
-            sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
-            sage: all( M.is_hermitian() for M in B )
-            True
-
-        """
-        from mjo.hurwitz import QuaternionMatrixAlgebra
-        A = QuaternionMatrixAlgebra(n, scalars=field)
-        es = A.entry_algebra_gens()
-
-        basis = []
-        for i in range(n):
-            for j in range(i+1):
-                if i == j:
-                    E_ii = A.monomial( (i,j,es[0]) )
-                    basis.append(E_ii)
-                else:
-                    for e in es:
-                        E_ij  = A.monomial( (i,j,e)             )
-                        ec = e.conjugate()
-                        # If the conjugate has a negative sign in front
-                        # of it, (j,i,ec) won't be a monomial!
-                        if (j,i,ec) in A.indices():
-                            E_ij += A.monomial( (j,i,ec) )
-                        else:
-                            E_ij -= A.monomial( (j,i,-ec) )
-                        basis.append(E_ij)
-
-        return tuple( basis )
-
-
-    @staticmethod
-    def trace_inner_product(X,Y):
-        r"""
-        Overload the superclass method because the quaternions are weird
-        and we need to use ``coefficient_tuple()`` to get the realpart.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
-
-        TESTS::
-
-            sage: J = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
-            sage: I = J.one().to_matrix()
-            sage: J.trace_inner_product(I, -I)
-            -2
-
-        """
-        return (X*Y).trace().coefficient_tuple()[0]
-

     def __init__(self, n, field=AA, **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
 
     def __init__(self, n, field=AA, **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
 

-        associative = False
-        if n <= 1:
-            associative = True
-
-        super().__init__(self._denormalized_basis(n,field),
+        from mjo.hurwitz import QuaternionMatrixAlgebra
+        A = QuaternionMatrixAlgebra(n, scalars=field)
+        super().__init__(self._denormalized_basis(A),

                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,
                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,

-                         associative=associative,

                          **kwargs)
 
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
                          **kwargs)
 
         # TODO: this could be factored out somehow, but is left here
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)

-        idV = self.matrix_space().one()
-        self.one.set_cache(self(idV))
+        self.one.set_cache(self(A.one()))

 
 
     @staticmethod
 
 
     @staticmethod


@@ -2311,7 +2265,9 @@ class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         # if the user passes check_axioms=True.
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 

-        super().__init__(self._denormalized_basis(n,field),
+        from mjo.hurwitz import OctonionMatrixAlgebra
+        A = OctonionMatrixAlgebra(n, scalars=field)
+        super().__init__(self._denormalized_basis(A),

                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,
                          self.jordan_product,
                          self.trace_inner_product,
                          field=field,


@@ -2321,72 +2277,7 @@ class OctonionHermitianEJA(RationalBasisEJA, ConcreteEJA, MatrixEJA):
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
         # because the MatrixEJA is not presently a subclass of the
         # FDEJA class that defines rank() and one().
         self.rank.set_cache(n)

-        idV = self.matrix_space().one()
-        self.one.set_cache(self(idV))
-
-
-    @classmethod
-    def _denormalized_basis(cls, n, field):
-        """
-        Returns a basis for the space of octonion Hermitian n-by-n
-        matrices.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
-
-        EXAMPLES::
-
-            sage: B = OctonionHermitianEJA._denormalized_basis(3,QQ)
-            sage: all( M.is_hermitian() for M in B )
-            True
-            sage: len(B)
-            27
-
-        """
-        from mjo.hurwitz import OctonionMatrixAlgebra
-        A = OctonionMatrixAlgebra(n, scalars=field)
-        es = A.entry_algebra_gens()
-
-        basis = []
-        for i in range(n):
-            for j in range(i+1):
-                if i == j:
-                    E_ii = A.monomial( (i,j,es[0]) )
-                    basis.append(E_ii)
-                else:
-                    for e in es:
-                        E_ij  = A.monomial( (i,j,e)             )
-                        ec = e.conjugate()
-                        # If the conjugate has a negative sign in front
-                        # of it, (j,i,ec) won't be a monomial!
-                        if (j,i,ec) in A.indices():
-                            E_ij += A.monomial( (j,i,ec) )
-                        else:
-                            E_ij -= A.monomial( (j,i,-ec) )
-                        basis.append(E_ij)
-
-        return tuple( basis )
-
-    @staticmethod
-    def trace_inner_product(X,Y):
-        r"""
-        The octonions don't know that the reals are embedded in them,
-        so we have to take the e0 component ourselves.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import OctonionHermitianEJA
-
-        TESTS::
-
-            sage: J = OctonionHermitianEJA(2,field=QQ,orthonormalize=False)
-            sage: I = J.one().to_matrix()
-            sage: J.trace_inner_product(I, -I)
-            -2
-
-        """
-        return (X*Y).trace().coefficient(0)
+        self.one.set_cache(self(A.one()))

 
 
 class AlbertEJA(OctonionHermitianEJA):
 
 
 class AlbertEJA(OctonionHermitianEJA):

diff --git a/mjo/hurwitz.py b/mjo/hurwitz.py
index ccc8219b1a92036c6ac92f118339c160a885977b..1f7c9dc3781a844d6d77a1e63168c52340bf6770 100644 (file)
--- a/mjo/hurwitz.py
+++ b/mjo/hurwitz.py
@@ -306,6 +306,31 @@ class Octonions(CombinatorialFreeModule):
 
 
 class HurwitzMatrixAlgebraElement(MatrixAlgebraElement):
 
 
 class HurwitzMatrixAlgebraElement(MatrixAlgebraElement):

+    def conjugate_transpose(self):
+        r"""
+        Return the conjugate-transpose of this matrix.
+
+        SETUP::
+
+            sage: from mjo.hurwitz import HurwitzMatrixAlgebra
+
+        EXAMPLES::
+
+            sage: A = HurwitzMatrixAlgebra(2, QQbar, ZZ)
+            sage: M = A([ [ I,   2*I],
+            ....:         [ 3*I, 4*I] ])
+            +------+------+
+            | -1*I | -3*I |
+            +------+------+
+            | -2*I | -4*I |
+            +------+------+
+
+        """
+        entries = [ [ self[j,i].conjugate()
+                      for j in range(self.ncols())]
+                    for i in range(self.nrows()) ]
+        return self.parent()._element_constructor_(entries)
+

     def is_hermitian(self):
         r"""
 
     def is_hermitian(self):
         r"""
 


@@ -322,6 +347,8 @@ class HurwitzMatrixAlgebraElement(MatrixAlgebraElement):
             True
 
         """
             True
 
         """

+        # A tiny bit faster than checking equality with the conjugate
+        # transpose.

         return all( self[i,j] == self[j,i].conjugate()
                     for i in range(self.nrows())
                     for j in range(self.ncols()) )
         return all( self[i,j] == self[j,i].conjugate()
                     for i in range(self.nrows())
                     for j in range(self.ncols()) )