mjo/ldlt.py: begin fast block-LDLT implementation.

[sage.d.git] / mjo / eja / eja_algebra.py

diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py
index 689a3db016437d1e6eda5c6372e52a3513896671..8bee7297fa170c0faf47610dc0ac435a881531fb 100644 (file)
--- a/mjo/eja/eja_algebra.py
+++ b/mjo/eja/eja_algebra.py
@@ -1994,23 +1994,20 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
 
         TESTS:
 
-        Ensure that this is one-half of the trace inner-product::
+        Ensure that this is one-half of the trace inner-product when
+        the algebra isn't just the reals (when ``n`` isn't one). This
+        is in Faraut and Koranyi, and also my "On the symmetry..."
+        paper::
 
             sage: set_random_seed()
-            sage: n = ZZ.random_element(5)
-            sage: M = matrix.random(QQ, n-1, algorithm='unimodular')
+            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: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
-            sage: V = J.vector_space()
-            sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
-            ....:         for ei in eis ]
-            sage: actual = [ sis[i]*sis[j]
-            ....:            for i in range(n-1)
-            ....:            for j in range(n-1) ]
-            sage: expected = [ J.one() if i == j else J.zero()
-            ....:              for i in range(n-1)
-            ....:              for j in range(n-1) ]
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: x.inner_product(y) == (x*y).trace()/2
+            True
 
         """
         xvec = x.to_vector()
@@ -2019,7 +2016,8 @@ class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
         ybar = yvec[1:]
         return x[0]*y[0] + (self._B*xbar).inner_product(ybar)
 
-class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+
+class JordanSpinEJA(BilinearFormEJA):
     """
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the usual inner product and jordan product ``x*y =
@@ -2056,42 +2054,9 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
         sage: JordanSpinEJA(2, prefix='B').gens()
         (B0, B1)
 
-    """
-    def __init__(self, n, field=QQ, **kwargs):
-        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:]
-                # z = x*y
-                z0 = x.inner_product(y)
-                zbar = y0*xbar + x0*ybar
-                z = V([z0] + zbar.list())
-                mult_table[i][j] = z
-
-        # The rank of the spin algebra is two, unless we're in a
-        # one-dimensional ambient space (because the rank is bounded by
-        # the ambient dimension).
-        fdeja = super(JordanSpinEJA, self)
-        return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
-
-    def inner_product(self, x, y):
-        """
-        Faster to reimplement than to use natural representations.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import JordanSpinEJA
-
-        TESTS:
+    TESTS:
 
-        Ensure that this is the usual inner product for the algebras
-        over `R^n`::
+        Ensure that we have the usual inner product on `R^n`::
 
             sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()
@@ -2101,8 +2066,11 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
             sage: x.inner_product(y) == J.natural_inner_product(X,Y)
             True
 
-        """
-        return x.to_vector().inner_product(y.to_vector())
+    """
+    def __init__(self, n, field=QQ, **kwargs):
+        # This is a special case of the BilinearFormEJA with the identity
+        # matrix as its bilinear form.
+        return super(JordanSpinEJA, self).__init__(n, field, **kwargs)
 
 
 class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):