eja: fix complex-unembedding with respect to 5537f4534.

[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 ca344edc709352dea2d1940a7646e058066d4b8d..1e5ada2188c3b6b3834165f221ab20b56d5c1f98 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")


@@ -160,6 +174,51 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         An element of a Euclidean Jordan algebra.
         """
 
         An element of a Euclidean Jordan algebra.
         """
 

+        def __init__(self, A, elt=None):
+            """
+            EXAMPLES:
+
+            The identity in `S^n` is converted to the identity in the EJA::
+
+                sage: J = RealSymmetricSimpleEJA(3)
+                sage: I = identity_matrix(QQ,3)
+                sage: J(I) == J.one()
+                True
+
+            This skew-symmetric matrix can't be represented in the EJA::
+
+                sage: J = RealSymmetricSimpleEJA(3)
+                sage: A = matrix(QQ,3, lambda i,j: i-j)
+                sage: J(A)
+                Traceback (most recent call last):
+                ...
+                ArithmeticError: vector is not in free module
+
+            """
+            # Goal: if we're given a matrix, and if it lives in our
+            # parent algebra's "natural ambient space," convert it
+            # into an algebra element.
+            #
+            # The catch is, we make a recursive call after converting
+            # the given matrix into a vector that lives in the algebra.
+            # This we need to try the parent class initializer first,
+            # to avoid recursing forever if we're given something that
+            # already fits into the algebra, but also happens to live
+            # in the parent's "natural ambient space" (this happens with
+            # vectors in R^n).
+            try:
+                FiniteDimensionalAlgebraElement.__init__(self, A, elt)
+            except ValueError:
+                natural_basis = A.natural_basis()
+                if elt in natural_basis[0].matrix_space():
+                    # Thanks for nothing! Matrix spaces aren't vector
+                    # spaces in Sage, so we have to figure out its
+                    # natural-basis coordinates ourselves.
+                    V = VectorSpace(elt.base_ring(), elt.nrows()**2)
+                    W = V.span( _mat2vec(s) for s in natural_basis )
+                    coords =  W.coordinates(_mat2vec(elt))
+                    FiniteDimensionalAlgebraElement.__init__(self, A, coords)
+

         def __pow__(self, n):
             """
             Return ``self`` raised to the power ``n``.
         def __pow__(self, n):
             """
             Return ``self`` raised to the power ``n``.


@@ -244,6 +303,24 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 sage: J(x).inner_product(J(y))
                 32
 
                 sage: J(x).inner_product(J(y))
                 32
 

+            The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
+            multiplication is the usual matrix multiplication in `S^n`,
+            so the inner product of the identity matrix with itself
+            should be the `n`::
+
+                sage: J = RealSymmetricSimpleEJA(3)
+                sage: J.one().inner_product(J.one())
+                3
+
+            Likewise, the inner product on `C^n` is `<X,Y> =
+            Re(trace(X*Y))`, where we must necessarily take the real
+            part because the product of Hermitian matrices may not be
+            Hermitian::
+
+                sage: J = ComplexHermitianSimpleEJA(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


@@ -941,13 +1018,10 @@ def eja_rn(dimension, field=QQ):
     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
            for i in xrange(dimension) ]
 
     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
            for i in xrange(dimension) ]
 

-    # The usual inner product on R^n.
-    ip = lambda x, y: x.vector().inner_product(y.vector())
-

     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=dimension,
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=dimension,

-                                                   inner_product=ip)
+                                                   inner_product=_usual_ip)

 
 
 
 
 
 


@@ -1042,6 +1116,12 @@ def _complex_hermitian_basis(n, field=QQ):
     return tuple(S)
 
 
     return tuple(S)
 
 

+def _mat2vec(m):
+        return vector(m.base_ring(), m.list())
+
+def _vec2mat(v):
+        return matrix(v.base_ring(), sqrt(v.degree()), v.list())
+

 def _multiplication_table_from_matrix_basis(basis):
     """
     At least three of the five simple Euclidean Jordan algebras have the
 def _multiplication_table_from_matrix_basis(basis):
     """
     At least three of the five simple Euclidean Jordan algebras have the


@@ -1062,19 +1142,13 @@ def _multiplication_table_from_matrix_basis(basis):
     field = basis[0].base_ring()
     dimension = basis[0].nrows()
 
     field = basis[0].base_ring()
     dimension = basis[0].nrows()
 

-    def mat2vec(m):
-        return vector(field, m.list())
-
-    def vec2mat(v):
-        return matrix(field, dimension, v.list())
-

     V = VectorSpace(field, dimension**2)
     V = VectorSpace(field, dimension**2)

-    W = V.span( mat2vec(s) for s in basis )
+    W = V.span( _mat2vec(s) for s in basis )

 
     # Taking the span above reorders our basis (thanks, jerk!) so we
     # need to put our "matrix basis" in the same order as the
     # (reordered) vector basis.
 
     # Taking the span above reorders our basis (thanks, jerk!) so we
     # need to put our "matrix basis" in the same order as the
     # (reordered) vector basis.

-    S = tuple( vec2mat(b) for b in W.basis() )
+    S = tuple( _vec2mat(b) for b in W.basis() )

 
     Qs = []
     for s in S:
 
     Qs = []
     for s in S:


@@ -1087,7 +1161,7 @@ def _multiplication_table_from_matrix_basis(basis):
         # why we're computing rows here and not columns.
         Q_rows = []
         for t in S:
         # why we're computing rows here and not columns.
         Q_rows = []
         for t in S:

-            this_row = mat2vec((s*t + t*s)/2)
+            this_row = _mat2vec((s*t + t*s)/2)

             Q_rows.append(W.coordinates(this_row))
         Q = matrix(field, W.dimension(), Q_rows)
         Qs.append(Q)
             Q_rows.append(W.coordinates(this_row))
         Q = matrix(field, W.dimension(), Q_rows)
         Qs.append(Q)


@@ -1108,13 +1182,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()


@@ -1125,7 +1199,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)


@@ -1142,8 +1216,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:


@@ -1161,14 +1244,26 @@ 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 ArgumentError('bad real submatrix')
+                raise ValueError('bad real submatrix')

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

-                raise ArgumentError('bad imag submatrix')
-            z = submat[0,0] + submat[1,0]*i
+                raise ValueError('bad imag 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)
 

+# The usual inner product on R^n.
+def _usual_ip(x,y):
+    return x.vector().inner_product(y.vector())
+
+# The inner product used for the real symmetric simple EJA.
+# We keep it as a separate function because e.g. the complex
+# algebra uses the same inner product, except divided by 2.
+def _matrix_ip(X,Y):
+    X_mat = X.natural_representation()
+    Y_mat = Y.natural_representation()
+    return (X_mat*Y_mat).trace()
+

 
 def RealSymmetricSimpleEJA(n, field=QQ):
     """
 
 def RealSymmetricSimpleEJA(n, field=QQ):
     """


@@ -1204,7 +1299,8 @@ def RealSymmetricSimpleEJA(n, field=QQ):
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=n,
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=n,

-                                                   natural_basis=T)
+                                                   natural_basis=T,
+                                                   inner_product=_matrix_ip)

 
 
 def ComplexHermitianSimpleEJA(n, field=QQ):
 
 
 def ComplexHermitianSimpleEJA(n, field=QQ):


@@ -1227,10 +1323,21 @@ def ComplexHermitianSimpleEJA(n, field=QQ):
     """
     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)

+
+    # Since a+bi on the diagonal is represented as
+    #
+    #   a + bi  = [  a  b  ]
+    #             [ -b  a  ],
+    #
+    # we'll double-count the "a" entries if we take the trace of
+    # the embedding.
+    ip = lambda X,Y: _matrix_ip(X,Y)/2
+

     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=n,
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=n,

-                                                   natural_basis=T)
+                                                   natural_basis=T,
+                                                   inner_product=ip)

 
 
 def QuaternionHermitianSimpleEJA(n):
 
 
 def QuaternionHermitianSimpleEJA(n):


@@ -1277,6 +1384,13 @@ def JordanSpinSimpleEJA(n, field=QQ):
         sage: e2*e3
         0
 
         sage: e2*e3
         0
 

+    In one dimension, this is the reals under multiplication::
+
+        sage: J1 = JordanSpinSimpleEJA(1)
+        sage: J2 = eja_rn(1)
+        sage: J1 == J2
+        True
+

     """
     Qs = []
     id_matrix = identity_matrix(field, n)
     """
     Qs = []
     id_matrix = identity_matrix(field, n)


@@ -1291,13 +1405,10 @@ def JordanSpinSimpleEJA(n, field=QQ):
         Qi[0,0] = Qi[0,0] * ~field(2)
         Qs.append(Qi)
 
         Qi[0,0] = Qi[0,0] * ~field(2)
         Qs.append(Qi)
 

-    # The usual inner product on R^n.
-    ip = lambda x, y: x.vector().inner_product(y.vector())
-

     # The rank of the spin factor algebra is two, UNLESS we're in a
     # one-dimensional ambient space (the rank is bounded by the
     # ambient dimension).
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=min(n,2),
     # The rank of the spin factor algebra is two, UNLESS we're in a
     # one-dimensional ambient space (the rank is bounded by the
     # ambient dimension).
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=min(n,2),

-                                                   inner_product=ip)
+                                                   inner_product=_usual_ip)