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 7c1f249b4632766e4b1808acc0275b2e55b8b506..1e5ada2188c3b6b3834165f221ab20b56d5c1f98 100644 (file)
--- a/mjo/eja/euclidean_jordan_algebra.py
+++ b/mjo/eja/euclidean_jordan_algebra.py
@@ -94,9 +94,23 @@ 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):
         """
         if (not x in self) or (not y in self):

-            raise ArgumentError("arguments must live in this algebra")
+            raise TypeError("arguments must live in this algebra")

         if self._inner_product is None:
             return x.trace_inner_product(y)
         else:
         if self._inner_product is None:
             return x.trace_inner_product(y)
         else:


@@ -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


@@ -259,7 +336,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             """
             P = self.parent()
             if not other in P:
             """
             P = self.parent()
             if not other in P:

-                raise ArgumentError("'other' must live in the same algebra")
+                raise TypeError("'other' must live in the same algebra")

 
             return P.inner_product(self, other)
 
 
             return P.inner_product(self, other)
 


@@ -294,7 +371,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             """
             if not other in self.parent():
 
             """
             if not other in self.parent():

-                raise ArgumentError("'other' must live in the same algebra")
+                raise TypeError("'other' must live in the same algebra")

 
             A = self.operator_matrix()
             B = other.operator_matrix()
 
             A = self.operator_matrix()
             B = other.operator_matrix()


@@ -384,7 +461,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             # TODO: we can do better once the call to is_invertible()
             # doesn't crash on irregular elements.
             #if not self.is_invertible():
             # TODO: we can do better once the call to is_invertible()
             # doesn't crash on irregular elements.
             #if not self.is_invertible():

-            #    raise ArgumentError('element is not invertible')
+            #    raise ValueError('element is not invertible')

 
             # We do this a little different than the usual recursive
             # call to a finite-dimensional algebra element, because we
 
             # We do this a little different than the usual recursive
             # call to a finite-dimensional algebra element, because we


@@ -740,7 +817,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             if other is None:
                 other=self
             elif not other in self.parent():
             if other is None:
                 other=self
             elif not other in self.parent():

-                raise ArgumentError("'other' must live in the same algebra")
+                raise TypeError("'other' must live in the same algebra")

 
             L = self.operator_matrix()
             M = other.operator_matrix()
 
             L = self.operator_matrix()
             M = other.operator_matrix()


@@ -903,7 +980,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             Return the trace inner product of myself and ``other``.
             """
             if not other in self.parent():
             Return the trace inner product of myself and ``other``.
             """
             if not other in self.parent():

-                raise ArgumentError("'other' must live in the same algebra")
+                raise TypeError("'other' must live in the same algebra")

 
             return (self*other).trace()
 
 
             return (self*other).trace()
 


@@ -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,24 +1182,24 @@ 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()
     if M.ncols() != n:
 
     """
     n = M.nrows()
     if M.ncols() != n:

-        raise ArgumentError("the matrix 'M' must be square")
+        raise ValueError("the matrix 'M' must be square")

     field = M.base_ring()
     blocks = []
     for z in M.list():
         a = z.real()
         b = z.imag()
     field = M.base_ring()
     blocks = []
     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,14 +1216,23 @@ 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:

-        raise ArgumentError("the matrix 'M' must be square")
+        raise ValueError("the matrix 'M' must be square")

     if not n.mod(2).is_zero():
     if not n.mod(2).is_zero():

-        raise ArgumentError("the matrix 'M' must be a complex embedding")
+        raise ValueError("the matrix 'M' must be a complex embedding")

 
     F = QuadraticField(-1, 'i')
     i = F.gen()
 
     F = QuadraticField(-1, 'i')
     i = F.gen()


@@ -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)