eja: get the quaternions working.

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

diff --git a/mjo/eja/euclidean_jordan_algebra.py b/mjo/eja/euclidean_jordan_algebra.py
index b25839815289ba6f60630a14b790dbfe29f8969b..b8c7a912b1ded8149f7aefdf14b90b99d1056b2b 100644 (file)
--- a/mjo/eja/euclidean_jordan_algebra.py
+++ b/mjo/eja/euclidean_jordan_algebra.py
@@ -94,6 +94,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         The inner product associated with this Euclidean Jordan algebra.
 
         Will default to the trace inner product if nothing else.
         The inner product associated with this Euclidean Jordan algebra.
 
         Will default to the trace inner product if nothing else.

+
+        EXAMPLES:
+
+        The inner product must satisfy its axiom for this algebra to truly
+        be a Euclidean Jordan Algebra::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: z = J.random_element()
+            sage: (x*y).inner_product(z) == y.inner_product(x*z)
+            True
+

         """
         if (not x in self) or (not y in self):
             raise TypeError("arguments must live in this algebra")
         """
         if (not x in self) or (not y in self):
             raise TypeError("arguments must live in this algebra")


@@ -307,6 +321,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 sage: J.one().inner_product(J.one())
                 3
 
                 sage: J.one().inner_product(J.one())
                 3
 

+            Ditto for the quaternions::
+
+                sage: J = QuaternionHermitianSimpleEJA(3)
+                sage: J.one().inner_product(J.one())
+                3
+

             TESTS:
 
             Ensure that we can always compute an inner product, and that
             TESTS:
 
             Ensure that we can always compute an inner product, and that


@@ -660,6 +680,25 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 [0 0 0 0 1 0]
                 [0 0 0 0 0 1]
 
                 [0 0 0 0 1 0]
                 [0 0 0 0 0 1]
 

+            ::
+
+                sage: J = QuaternionHermitianSimpleEJA(3)
+                sage: J.one()
+                e0 + e9 + e14
+                sage: J.one().natural_representation()
+                [1 0 0 0 0 0 0 0 0 0 0 0]
+                [0 1 0 0 0 0 0 0 0 0 0 0]
+                [0 0 1 0 0 0 0 0 0 0 0 0]
+                [0 0 0 1 0 0 0 0 0 0 0 0]
+                [0 0 0 0 1 0 0 0 0 0 0 0]
+                [0 0 0 0 0 1 0 0 0 0 0 0]
+                [0 0 0 0 0 0 1 0 0 0 0 0]
+                [0 0 0 0 0 0 0 1 0 0 0 0]
+                [0 0 0 0 0 0 0 0 1 0 0 0]
+                [0 0 0 0 0 0 0 0 0 1 0 0]
+                [0 0 0 0 0 0 0 0 0 0 1 0]
+                [0 0 0 0 0 0 0 0 0 0 0 1]
+

             """
             B = self.parent().natural_basis()
             W = B[0].matrix_space()
             """
             B = self.parent().natural_basis()
             W = B[0].matrix_space()


@@ -1028,6 +1067,12 @@ def random_eja():
       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
         product.
 
       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
         product.
 

+      * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
+        in the space of ``2n``-by-``2n`` real symmetric matrices.
+
+      * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
+        in the space of ``4n``-by-``4n`` real symmetric matrices.
+

     Later this might be extended to return Cartesian products of the
     EJAs above.
 
     Later this might be extended to return Cartesian products of the
     EJAs above.
 


@@ -1041,7 +1086,8 @@ def random_eja():
     constructor = choice([eja_rn,
                           JordanSpinSimpleEJA,
                           RealSymmetricSimpleEJA,
     constructor = choice([eja_rn,
                           JordanSpinSimpleEJA,
                           RealSymmetricSimpleEJA,

-                          ComplexHermitianSimpleEJA])
+                          ComplexHermitianSimpleEJA,
+                          QuaternionHermitianSimpleEJA])

     return constructor(n, field=QQ)
 
 
     return constructor(n, field=QQ)
 
 


@@ -1102,6 +1148,48 @@ def _complex_hermitian_basis(n, field=QQ):
     return tuple(S)
 
 
     return tuple(S)
 
 

+def _quaternion_hermitian_basis(n, field=QQ):
+    """
+    Returns a basis for the space of quaternion Hermitian n-by-n matrices.
+
+    TESTS::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
+        True
+
+    """
+    Q = QuaternionAlgebra(QQ,-1,-1)
+    I,J,K = Q.gens()
+
+    # This is like the symmetric case, but we need to be careful:
+    #
+    #   * We want conjugate-symmetry, not just symmetry.
+    #   * The diagonal will (as a result) be real.
+    #
+    S = []
+    for i in xrange(n):
+        for j in xrange(i+1):
+            Eij = matrix(Q, n, lambda k,l: k==i and l==j)
+            if i == j:
+                Sij = _embed_quaternion_matrix(Eij)
+                S.append(Sij)
+            else:
+                # Beware, orthogonal but not normalized! The second,
+                # third, and fourth ones have a minus because they're
+                # conjugated.
+                Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
+                S.append(Sij_real)
+                Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
+                S.append(Sij_I)
+                Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
+                S.append(Sij_J)
+                Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
+                S.append(Sij_K)
+    return tuple(S)
+
+

 def _mat2vec(m):
         return vector(m.base_ring(), m.list())
 
 def _mat2vec(m):
         return vector(m.base_ring(), m.list())
 


@@ -1168,13 +1256,13 @@ def _embed_complex_matrix(M):
         sage: x2 = F(1 + 2*i)
         sage: x3 = F(-i)
         sage: x4 = F(6)
         sage: x2 = F(1 + 2*i)
         sage: x3 = F(-i)
         sage: x4 = F(6)

-        sage: M = matrix(F,2,[x1,x2,x3,x4])
+        sage: M = matrix(F,2,[[x1,x2],[x3,x4]])

         sage: _embed_complex_matrix(M)
         sage: _embed_complex_matrix(M)

-        [ 4  2| 1 -2]
-        [-2  4| 2  1]
+        [ 4 -2| 1  2]
+        [ 2  4|-2  1]

         [-----+-----]
         [-----+-----]

-        [ 0  1| 6  0]
-        [-1  0| 0  6]
+        [ 0 -1| 6  0]
+        [ 1  0| 0  6]

 
     """
     n = M.nrows()
 
     """
     n = M.nrows()


@@ -1185,7 +1273,7 @@ def _embed_complex_matrix(M):
     for z in M.list():
         a = z.real()
         b = z.imag()
     for z in M.list():
         a = z.real()
         b = z.imag()

-        blocks.append(matrix(field, 2, [[a,-b],[b,a]]))
+        blocks.append(matrix(field, 2, [[a,b],[-b,a]]))

 
     # We can drop the imaginaries here.
     return block_matrix(field.base_ring(), n, blocks)
 
     # We can drop the imaginaries here.
     return block_matrix(field.base_ring(), n, blocks)


@@ -1202,8 +1290,17 @@ def _unembed_complex_matrix(M):
         ....:                 [ 9,  10, 11, 12],
         ....:                 [-10, 9, -12, 11] ])
         sage: _unembed_complex_matrix(A)
         ....:                 [ 9,  10, 11, 12],
         ....:                 [-10, 9, -12, 11] ])
         sage: _unembed_complex_matrix(A)

-        [  -2*i + 1   -4*i + 3]
-        [ -10*i + 9 -12*i + 11]
+        [  2*i + 1   4*i + 3]
+        [ 10*i + 9 12*i + 11]
+
+    TESTS::
+
+        sage: set_random_seed()
+        sage: F = QuadraticField(-1, 'i')
+        sage: M = random_matrix(F, 3)
+        sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
+        True
+

     """
     n = ZZ(M.nrows())
     if M.ncols() != n:
     """
     n = ZZ(M.nrows())
     if M.ncols() != n:


@@ -1221,14 +1318,109 @@ def _unembed_complex_matrix(M):
         for j in xrange(n/2):
             submat = M[2*k:2*k+2,2*j:2*j+2]
             if submat[0,0] != submat[1,1]:
         for j in xrange(n/2):
             submat = M[2*k:2*k+2,2*j:2*j+2]
             if submat[0,0] != submat[1,1]:

-                raise ValueError('bad real submatrix')
+                raise ValueError('bad on-diagonal submatrix')

             if submat[0,1] != -submat[1,0]:
             if submat[0,1] != -submat[1,0]:

-                raise ValueError('bad imag submatrix')
-            z = submat[0,0] + submat[1,0]*i
+                raise ValueError('bad off-diagonal submatrix')
+            z = submat[0,0] + submat[0,1]*i

             elements.append(z)
 
     return matrix(F, n/2, elements)
 
             elements.append(z)
 
     return matrix(F, n/2, elements)
 

+
+def _embed_quaternion_matrix(M):
+    """
+    Embed the n-by-n quaternion matrix ``M`` into the space of real
+    matrices of size 4n-by-4n by first sending each quaternion entry
+    `z = a + bi + cj + dk` to the block-complex matrix
+    ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
+    a real matrix.
+
+    EXAMPLES::
+
+        sage: Q = QuaternionAlgebra(QQ,-1,-1)
+        sage: i,j,k = Q.gens()
+        sage: x = 1 + 2*i + 3*j + 4*k
+        sage: M = matrix(Q, 1, [[x]])
+        sage: _embed_quaternion_matrix(M)
+        [ 1  2  3  4]
+        [-2  1 -4  3]
+        [-3  4  1 -2]
+        [-4 -3  2  1]
+
+    """
+    quaternions = M.base_ring()
+    n = M.nrows()
+    if M.ncols() != n:
+        raise ValueError("the matrix 'M' must be square")
+
+    F = QuadraticField(-1, 'i')
+    i = F.gen()
+
+    blocks = []
+    for z in M.list():
+        t = z.coefficient_tuple()
+        a = t[0]
+        b = t[1]
+        c = t[2]
+        d = t[3]
+        cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
+                                    [-c + d*i, a - b*i]])
+        blocks.append(_embed_complex_matrix(cplx_matrix))
+
+    # We should have real entries by now, so use the realest field
+    # we've got for the return value.
+    return block_matrix(quaternions.base_ring(), n, blocks)
+
+
+def _unembed_quaternion_matrix(M):
+    """
+    The inverse of _embed_quaternion_matrix().
+
+    EXAMPLES::
+
+        sage: M = matrix(QQ, [[ 1,  2,  3,  4],
+        ....:                 [-2,  1, -4,  3],
+        ....:                 [-3,  4,  1, -2],
+        ....:                 [-4, -3,  2,  1]])
+        sage: _unembed_quaternion_matrix(M)
+        [1 + 2*i + 3*j + 4*k]
+
+    TESTS::
+
+        sage: set_random_seed()
+        sage: Q = QuaternionAlgebra(QQ, -1, -1)
+        sage: M = random_matrix(Q, 3)
+        sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
+        True
+
+    """
+    n = ZZ(M.nrows())
+    if M.ncols() != n:
+        raise ValueError("the matrix 'M' must be square")
+    if not n.mod(4).is_zero():
+        raise ValueError("the matrix 'M' must be a complex embedding")
+
+    Q = QuaternionAlgebra(QQ,-1,-1)
+    i,j,k = Q.gens()
+
+    # Go top-left to bottom-right (reading order), converting every
+    # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
+    # quaternion block.
+    elements = []
+    for l in xrange(n/4):
+        for m in xrange(n/4):
+            submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
+            if submat[0,0] != submat[1,1].conjugate():
+                raise ValueError('bad on-diagonal submatrix')
+            if submat[0,1] != -submat[1,0].conjugate():
+                raise ValueError('bad off-diagonal submatrix')
+            z  = submat[0,0].real() + submat[0,0].imag()*i
+            z += submat[0,1].real()*j + submat[0,1].imag()*k
+            elements.append(z)
+
+    return matrix(Q, n/4, elements)
+
+

 # The usual inner product on R^n.
 def _usual_ip(x,y):
     return x.vector().inner_product(y.vector())
 # The usual inner product on R^n.
 def _usual_ip(x,y):
     return x.vector().inner_product(y.vector())


@@ -1269,6 +1461,22 @@ def RealSymmetricSimpleEJA(n, field=QQ):
         sage: J.degree() == (n^2 + n)/2
         True
 
         sage: J.degree() == (n^2 + n)/2
         True
 

+    The Jordan multiplication is what we think it is::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = RealSymmetricSimpleEJA(n)
+        sage: x = J.random_element()
+        sage: y = J.random_element()
+        sage: actual = (x*y).natural_representation()
+        sage: X = x.natural_representation()
+        sage: Y = y.natural_representation()
+        sage: expected = (X*Y + Y*X)/2
+        sage: actual == expected
+        True
+        sage: J(expected) == x*y
+        True
+

     """
     S = _real_symmetric_basis(n, field=field)
     (Qs, T) = _multiplication_table_from_matrix_basis(S)
     """
     S = _real_symmetric_basis(n, field=field)
     (Qs, T) = _multiplication_table_from_matrix_basis(S)


@@ -1297,6 +1505,22 @@ def ComplexHermitianSimpleEJA(n, field=QQ):
         sage: J.degree() == n^2
         True
 
         sage: J.degree() == n^2
         True
 

+    The Jordan multiplication is what we think it is::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = ComplexHermitianSimpleEJA(n)
+        sage: x = J.random_element()
+        sage: y = J.random_element()
+        sage: actual = (x*y).natural_representation()
+        sage: X = x.natural_representation()
+        sage: Y = y.natural_representation()
+        sage: expected = (X*Y + Y*X)/2
+        sage: actual == expected
+        True
+        sage: J(expected) == x*y
+        True
+

     """
     S = _complex_hermitian_basis(n)
     (Qs, T) = _multiplication_table_from_matrix_basis(S)
     """
     S = _complex_hermitian_basis(n)
     (Qs, T) = _multiplication_table_from_matrix_basis(S)


@@ -1317,14 +1541,60 @@ def ComplexHermitianSimpleEJA(n, field=QQ):
                                                    inner_product=ip)
 
 
                                                    inner_product=ip)
 
 

-def QuaternionHermitianSimpleEJA(n):
+def QuaternionHermitianSimpleEJA(n, field=QQ):

     """
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
     real-part-of-trace inner product. It has dimension `2n^2 - n` over
     the reals.
     """
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
     real-part-of-trace inner product. It has dimension `2n^2 - n` over
     the reals.

+
+    TESTS:
+
+    The degree of this algebra is `n^2`::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = QuaternionHermitianSimpleEJA(n)
+        sage: J.degree() == 2*(n^2) - n
+        True
+
+    The Jordan multiplication is what we think it is::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = QuaternionHermitianSimpleEJA(n)
+        sage: x = J.random_element()
+        sage: y = J.random_element()
+        sage: actual = (x*y).natural_representation()
+        sage: X = x.natural_representation()
+        sage: Y = y.natural_representation()
+        sage: expected = (X*Y + Y*X)/2
+        sage: actual == expected
+        True
+        sage: J(expected) == x*y
+        True
+

     """
     """

-    pass
+    S = _quaternion_hermitian_basis(n)
+    (Qs, T) = _multiplication_table_from_matrix_basis(S)
+
+    # Since a+bi+cj+dk on the diagonal is represented as
+    #
+    #   a + bi +cj + dk = [  a  b  c  d]
+    #                     [ -b  a -d  c]
+    #                     [ -c  d  a -b]
+    #                     [ -d -c  b  a],
+    #
+    # we'll quadruple-count the "a" entries if we take the trace of
+    # the embedding.
+    ip = lambda X,Y: _matrix_ip(X,Y)/4
+
+    return FiniteDimensionalEuclideanJordanAlgebra(field,
+                                                   Qs,
+                                                   rank=n,
+                                                   natural_basis=T,
+                                                   inner_product=ip)
+

 
 def OctonionHermitianSimpleEJA(n):
     """
 
 def OctonionHermitianSimpleEJA(n):
     """