]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: add random_instance() for BilinearFormEJA.
[sage.d.git] / mjo / eja / eja_algebra.py
index aec6b50e7f5c31bfb73cb8de06ec4863cc454a66..8ca501d6ba944a539901d384a78a5c8a0f45a7ca 100644 (file)
@@ -233,7 +233,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         interpreted to be far less than the dimension) should override
         with a smaller number.
         """
-        return 5
+        raise NotImplementedError
 
     def _repr_(self):
         """
@@ -599,19 +599,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # appeal to the "long vectors" isometry.
         oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ]
 
-        # Now we use basis linear algebra to find the coefficients,
+        # Now we use basic linear algebra to find the coefficients,
         # of the matrices-as-vectors-linear-combination, which should
         # work for the original algebra basis too.
-        A = matrix.column(self.base_ring(), oper_vecs)
+        A = matrix(self.base_ring(), oper_vecs)
 
         # We used the isometry on the left-hand side already, but we
         # still need to do it for the right-hand side. Recall that we
         # wanted something that summed to the identity matrix.
         b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) )
 
-        # Now if there's an identity element in the algebra, this should work.
-        coeffs = A.solve_right(b)
-        return self.linear_combination(zip(self.gens(), coeffs))
+        # Now if there's an identity element in the algebra, this
+        # should work. We solve on the left to avoid having to
+        # transpose the matrix "A".
+        return self.from_vector(A.solve_left(b))
 
 
     def peirce_decomposition(self, c):
@@ -839,11 +840,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         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_random_instance_size() + 1)
         return cls(n, field, **kwargs)
 
@@ -2056,6 +2052,10 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra):
                                           **kwargs)
         self.rank.set_cache(n)
 
+    @staticmethod
+    def _max_random_instance_size():
+        return 5
+
     def inner_product(self, x, y):
         """
         Faster to reimplement than to use natural representations.
@@ -2085,10 +2085,15 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
     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 =
-    (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
-    symmetric positive-definite "bilinear form" matrix. It has
-    dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
-    when ``B`` is the identity matrix of order ``n-1``.
+    (<Bx,y>,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is
+    a symmetric positive-definite "bilinear form" matrix. Its
+    dimension is the size of `B`, and it has rank two in dimensions
+    larger than two. It reduces to the ``JordanSpinEJA`` when `B` is
+    the identity matrix of order ``n``.
+
+    We insist that the one-by-one upper-left identity block of `B` be
+    passed in as well so that we can be passed a matrix of size zero
+    to construct a trivial algebra.
 
     SETUP::
 
@@ -2100,7 +2105,8 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
     When no bilinear form is specified, the identity matrix is used,
     and the resulting algebra is the Jordan spin algebra::
 
-        sage: J0 = BilinearFormEJA(3)
+        sage: B = matrix.identity(AA,3)
+        sage: J0 = BilinearFormEJA(B)
         sage: J1 = JordanSpinEJA(3)
         sage: J0.multiplication_table() == J0.multiplication_table()
         True
@@ -2109,7 +2115,8 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
 
     We can create a zero-dimensional algebra::
 
-        sage: J = BilinearFormEJA(0)
+        sage: B = matrix.identity(AA,0)
+        sage: J = BilinearFormEJA(B)
         sage: J.basis()
         Finite family {}
 
@@ -2121,8 +2128,11 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
         sage: set_random_seed()
         sage: n = ZZ.random_element(5)
         sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
-        sage: B = M.transpose()*M
-        sage: J = BilinearFormEJA(n, B=B)
+        sage: B11 = matrix.identity(QQ,1)
+        sage: B22 = M.transpose()*M
+        sage: B = block_matrix(2,2,[ [B11,0  ],
+        ....:                        [0, B22 ] ])
+        sage: J = BilinearFormEJA(B)
         sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
         sage: V = J.vector_space()
         sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
@@ -2136,11 +2146,12 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
         sage: actual == expected
         True
     """
-    def __init__(self, n, field=AA, B=None, **kwargs):
-        if B is None:
-            self._B = matrix.identity(field, max(0,n-1))
-        else:
-            self._B = B
+    def __init__(self, B, field=AA, **kwargs):
+        self._B = B
+        n = B.nrows()
+
+        if not B.is_positive_definite():
+            raise TypeError("matrix B is not positive-definite")
 
         V = VectorSpace(field, n)
         mult_table = [[V.zero() for j in range(n)] for i in range(n)]
@@ -2152,7 +2163,7 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
                 xbar = x[1:]
                 y0 = y[0]
                 ybar = y[1:]
-                z0 = x0*y0 + (self._B*xbar).inner_product(ybar)
+                z0 = (B*x).inner_product(y)
                 zbar = y0*xbar + x0*ybar
                 z = V([z0] + zbar.list())
                 mult_table[i][j] = z
@@ -2166,6 +2177,35 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
                                               **kwargs)
         self.rank.set_cache(min(n,2))
 
+    @staticmethod
+    def _max_random_instance_size():
+        return 5
+
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        """
+        Return a random instance of this algebra.
+        """
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        if n == 0:
+            # Special case needed since we use (n-1) below.
+            B = matrix.identity(field, 0)
+            return cls(B, field, **kwargs)
+
+        B11 = matrix.identity(field,1)
+        M = matrix.random(field, n-1)
+        I = matrix.identity(field, n-1)
+        alpha = field.zero()
+        while alpha.is_zero():
+            alpha = field.random_element().abs()
+        B22 = M.transpose()*M + alpha*I
+
+        from sage.matrix.special import block_matrix
+        B = block_matrix(2,2, [ [B11,   ZZ(0) ],
+                                [ZZ(0), B22 ] ])
+
+        return cls(B, field, **kwargs)
+
     def inner_product(self, x, y):
         r"""
         Half of the trace inner product.
@@ -2185,21 +2225,15 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
         paper::
 
             sage: set_random_seed()
-            sage: n = ZZ.random_element(2,5)
-            sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
-            sage: B = M.transpose()*M
-            sage: J = BilinearFormEJA(n, B=B)
+            sage: J = BilinearFormEJA.random_instance()
+            sage: n = J.dimension()
             sage: x = J.random_element()
             sage: y = J.random_element()
-            sage: x.inner_product(y) == (x*y).trace()/2
+            sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
             True
 
         """
-        xvec = x.to_vector()
-        xbar = xvec[1:]
-        yvec = y.to_vector()
-        ybar = yvec[1:]
-        return x[0]*y[0] + (self._B*xbar).inner_product(ybar)
+        return (self._B*x.to_vector()).inner_product(y.to_vector())
 
 
 class JordanSpinEJA(BilinearFormEJA):
@@ -2255,7 +2289,8 @@ class JordanSpinEJA(BilinearFormEJA):
     def __init__(self, n, field=AA, **kwargs):
         # This is a special case of the BilinearFormEJA with the identity
         # matrix as its bilinear form.
-        super(JordanSpinEJA, self).__init__(n, field, **kwargs)
+        B = matrix.identity(field, n)
+        super(JordanSpinEJA, self).__init__(B, field, **kwargs)
 
 
 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
@@ -2297,6 +2332,11 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
         # largest subalgebra generated by any element.
         self.rank.set_cache(0)
 
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        # We don't take a "size" argument so the superclass method is
+        # inappropriate for us.
+        return cls(field, **kwargs)
 
 class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
     r"""