2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
4 are used in optimization, and have some additional nice methods beyond
5 what can be supported in a general Jordan Algebra.
8 from itertools
import izip
, repeat
10 from sage
.algebras
.quatalg
.quaternion_algebra
import QuaternionAlgebra
11 from sage
.categories
.magmatic_algebras
import MagmaticAlgebras
12 from sage
.combinat
.free_module
import CombinatorialFreeModule
13 from sage
.matrix
.constructor
import matrix
14 from sage
.matrix
.matrix_space
import MatrixSpace
15 from sage
.misc
.cachefunc
import cached_method
16 from sage
.misc
.prandom
import choice
17 from sage
.misc
.table
import table
18 from sage
.modules
.free_module
import FreeModule
, VectorSpace
19 from sage
.rings
.integer_ring
import ZZ
20 from sage
.rings
.number_field
.number_field
import NumberField
, QuadraticField
21 from sage
.rings
.polynomial
.polynomial_ring_constructor
import PolynomialRing
22 from sage
.rings
.rational_field
import QQ
23 from sage
.rings
.real_lazy
import CLF
, RLF
24 from sage
.structure
.element
import is_Matrix
26 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
27 from mjo
.eja
.eja_utils
import _mat2vec
29 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
30 # This is an ugly hack needed to prevent the category framework
31 # from implementing a coercion from our base ring (e.g. the
32 # rationals) into the algebra. First of all -- such a coercion is
33 # nonsense to begin with. But more importantly, it tries to do so
34 # in the category of rings, and since our algebras aren't
35 # associative they generally won't be rings.
36 _no_generic_basering_coercion
= True
48 sage: from mjo.eja.eja_algebra import random_eja
52 By definition, Jordan multiplication commutes::
54 sage: set_random_seed()
55 sage: J = random_eja()
56 sage: x,y = J.random_elements(2)
62 self
._natural
_basis
= natural_basis
65 category
= MagmaticAlgebras(field
).FiniteDimensional()
66 category
= category
.WithBasis().Unital()
68 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
70 range(len(mult_table
)),
73 self
.print_options(bracket
='')
75 # The multiplication table we're given is necessarily in terms
76 # of vectors, because we don't have an algebra yet for
77 # anything to be an element of. However, it's faster in the
78 # long run to have the multiplication table be in terms of
79 # algebra elements. We do this after calling the superclass
80 # constructor so that from_vector() knows what to do.
81 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
82 for ls
in mult_table
]
85 def _element_constructor_(self
, elt
):
87 Construct an element of this algebra from its natural
90 This gets called only after the parent element _call_ method
91 fails to find a coercion for the argument.
95 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
96 ....: RealCartesianProductEJA,
97 ....: RealSymmetricEJA)
101 The identity in `S^n` is converted to the identity in the EJA::
103 sage: J = RealSymmetricEJA(3)
104 sage: I = matrix.identity(QQ,3)
105 sage: J(I) == J.one()
108 This skew-symmetric matrix can't be represented in the EJA::
110 sage: J = RealSymmetricEJA(3)
111 sage: A = matrix(QQ,3, lambda i,j: i-j)
113 Traceback (most recent call last):
115 ArithmeticError: vector is not in free module
119 Ensure that we can convert any element of the two non-matrix
120 simple algebras (whose natural representations are their usual
121 vector representations) back and forth faithfully::
123 sage: set_random_seed()
124 sage: J = RealCartesianProductEJA.random_instance()
125 sage: x = J.random_element()
126 sage: J(x.to_vector().column()) == x
128 sage: J = JordanSpinEJA.random_instance()
129 sage: x = J.random_element()
130 sage: J(x.to_vector().column()) == x
135 # The superclass implementation of random_element()
136 # needs to be able to coerce "0" into the algebra.
139 natural_basis
= self
.natural_basis()
140 basis_space
= natural_basis
[0].matrix_space()
141 if elt
not in basis_space
:
142 raise ValueError("not a naturally-represented algebra element")
144 # Thanks for nothing! Matrix spaces aren't vector spaces in
145 # Sage, so we have to figure out its natural-basis coordinates
146 # ourselves. We use the basis space's ring instead of the
147 # element's ring because the basis space might be an algebraic
148 # closure whereas the base ring of the 3-by-3 identity matrix
149 # could be QQ instead of QQbar.
150 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
151 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
152 coords
= W
.coordinate_vector(_mat2vec(elt
))
153 return self
.from_vector(coords
)
157 def _max_test_case_size():
159 Return an integer "size" that is an upper bound on the size of
160 this algebra when it is used in a random test
161 case. Unfortunately, the term "size" is quite vague -- when
162 dealing with `R^n` under either the Hadamard or Jordan spin
163 product, the "size" refers to the dimension `n`. When dealing
164 with a matrix algebra (real symmetric or complex/quaternion
165 Hermitian), it refers to the size of the matrix, which is
166 far less than the dimension of the underlying vector space.
168 We default to five in this class, which is safe in `R^n`. The
169 matrix algebra subclasses (or any class where the "size" is
170 interpreted to be far less than the dimension) should override
171 with a smaller number.
178 Return a string representation of ``self``.
182 sage: from mjo.eja.eja_algebra import JordanSpinEJA
186 Ensure that it says what we think it says::
188 sage: JordanSpinEJA(2, field=QQ)
189 Euclidean Jordan algebra of dimension 2 over Rational Field
190 sage: JordanSpinEJA(3, field=RDF)
191 Euclidean Jordan algebra of dimension 3 over Real Double Field
194 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
195 return fmt
.format(self
.dimension(), self
.base_ring())
197 def product_on_basis(self
, i
, j
):
198 return self
._multiplication
_table
[i
][j
]
200 def _a_regular_element(self
):
202 Guess a regular element. Needed to compute the basis for our
203 characteristic polynomial coefficients.
207 sage: from mjo.eja.eja_algebra import random_eja
211 Ensure that this hacky method succeeds for every algebra that we
212 know how to construct::
214 sage: set_random_seed()
215 sage: J = random_eja()
216 sage: J._a_regular_element().is_regular()
221 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
222 if not z
.is_regular():
223 raise ValueError("don't know a regular element")
228 def _charpoly_basis_space(self
):
230 Return the vector space spanned by the basis used in our
231 characteristic polynomial coefficients. This is used not only to
232 compute those coefficients, but also any time we need to
233 evaluate the coefficients (like when we compute the trace or
236 z
= self
._a
_regular
_element
()
237 # Don't use the parent vector space directly here in case this
238 # happens to be a subalgebra. In that case, we would be e.g.
239 # two-dimensional but span_of_basis() would expect three
241 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
242 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
243 V1
= V
.span_of_basis( basis
)
244 b
= (V1
.basis() + V1
.complement().basis())
245 return V
.span_of_basis(b
)
250 def _charpoly_coeff(self
, i
):
252 Return the coefficient polynomial "a_{i}" of this algebra's
253 general characteristic polynomial.
255 Having this be a separate cached method lets us compute and
256 store the trace/determinant (a_{r-1} and a_{0} respectively)
257 separate from the entire characteristic polynomial.
259 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
260 R
= A_of_x
.base_ring()
262 # Guaranteed by theory
265 # Danger: the in-place modification is done for performance
266 # reasons (reconstructing a matrix with huge polynomial
267 # entries is slow), but I don't know how cached_method works,
268 # so it's highly possible that we're modifying some global
269 # list variable by reference, here. In other words, you
270 # probably shouldn't call this method twice on the same
271 # algebra, at the same time, in two threads
272 Ai_orig
= A_of_x
.column(i
)
273 A_of_x
.set_column(i
,xr
)
274 numerator
= A_of_x
.det()
275 A_of_x
.set_column(i
,Ai_orig
)
277 # We're relying on the theory here to ensure that each a_i is
278 # indeed back in R, and the added negative signs are to make
279 # the whole charpoly expression sum to zero.
280 return R(-numerator
/detA
)
284 def _charpoly_matrix_system(self
):
286 Compute the matrix whose entries A_ij are polynomials in
287 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
288 corresponding to `x^r` and the determinent of the matrix A =
289 [A_ij]. In other words, all of the fixed (cachable) data needed
290 to compute the coefficients of the characteristic polynomial.
295 # Turn my vector space into a module so that "vectors" can
296 # have multivatiate polynomial entries.
297 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
298 R
= PolynomialRing(self
.base_ring(), names
)
300 # Using change_ring() on the parent's vector space doesn't work
301 # here because, in a subalgebra, that vector space has a basis
302 # and change_ring() tries to bring the basis along with it. And
303 # that doesn't work unless the new ring is a PID, which it usually
307 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
311 # And figure out the "left multiplication by x" matrix in
314 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
315 for i
in range(n
) ] # don't recompute these!
317 ek
= self
.monomial(k
).to_vector()
319 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
320 for i
in range(n
) ) )
321 Lx
= matrix
.column(R
, lmbx_cols
)
323 # Now we can compute powers of x "symbolically"
324 x_powers
= [self
.one().to_vector(), x
]
325 for d
in range(2, r
+1):
326 x_powers
.append( Lx
*(x_powers
[-1]) )
328 idmat
= matrix
.identity(R
, n
)
330 W
= self
._charpoly
_basis
_space
()
331 W
= W
.change_ring(R
.fraction_field())
333 # Starting with the standard coordinates x = (X1,X2,...,Xn)
334 # and then converting the entries to W-coordinates allows us
335 # to pass in the standard coordinates to the charpoly and get
336 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
339 # W.coordinates(x^2) eval'd at (standard z-coords)
343 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
345 # We want the middle equivalent thing in our matrix, but use
346 # the first equivalent thing instead so that we can pass in
347 # standard coordinates.
348 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
349 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
350 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
351 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
355 def characteristic_polynomial(self
):
357 Return a characteristic polynomial that works for all elements
360 The resulting polynomial has `n+1` variables, where `n` is the
361 dimension of this algebra. The first `n` variables correspond to
362 the coordinates of an algebra element: when evaluated at the
363 coordinates of an algebra element with respect to a certain
364 basis, the result is a univariate polynomial (in the one
365 remaining variable ``t``), namely the characteristic polynomial
370 sage: from mjo.eja.eja_algebra import JordanSpinEJA
374 The characteristic polynomial in the spin algebra is given in
375 Alizadeh, Example 11.11::
377 sage: J = JordanSpinEJA(3)
378 sage: p = J.characteristic_polynomial(); p
379 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
380 sage: xvec = J.one().to_vector()
388 # The list of coefficient polynomials a_1, a_2, ..., a_n.
389 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
391 # We go to a bit of trouble here to reorder the
392 # indeterminates, so that it's easier to evaluate the
393 # characteristic polynomial at x's coordinates and get back
394 # something in terms of t, which is what we want.
396 S
= PolynomialRing(self
.base_ring(),'t')
398 S
= PolynomialRing(S
, R
.variable_names())
401 # Note: all entries past the rth should be zero. The
402 # coefficient of the highest power (x^r) is 1, but it doesn't
403 # appear in the solution vector which contains coefficients
404 # for the other powers (to make them sum to x^r).
406 a
[r
] = 1 # corresponds to x^r
408 # When the rank is equal to the dimension, trying to
409 # assign a[r] goes out-of-bounds.
410 a
.append(1) # corresponds to x^r
412 return sum( a
[k
]*(t
**k
) for k
in xrange(len(a
)) )
415 def inner_product(self
, x
, y
):
417 The inner product associated with this Euclidean Jordan algebra.
419 Defaults to the trace inner product, but can be overridden by
420 subclasses if they are sure that the necessary properties are
425 sage: from mjo.eja.eja_algebra import random_eja
429 Our inner product satisfies the Jordan axiom, which is also
430 referred to as "associativity" for a symmetric bilinear form::
432 sage: set_random_seed()
433 sage: J = random_eja()
434 sage: x,y,z = J.random_elements(3)
435 sage: (x*y).inner_product(z) == y.inner_product(x*z)
439 X
= x
.natural_representation()
440 Y
= y
.natural_representation()
441 return self
.natural_inner_product(X
,Y
)
444 def is_trivial(self
):
446 Return whether or not this algebra is trivial.
448 A trivial algebra contains only the zero element.
452 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
456 sage: J = ComplexHermitianEJA(3)
459 sage: A = J.zero().subalgebra_generated_by()
464 return self
.dimension() == 0
467 def multiplication_table(self
):
469 Return a visual representation of this algebra's multiplication
470 table (on basis elements).
474 sage: from mjo.eja.eja_algebra import JordanSpinEJA
478 sage: J = JordanSpinEJA(4)
479 sage: J.multiplication_table()
480 +----++----+----+----+----+
481 | * || e0 | e1 | e2 | e3 |
482 +====++====+====+====+====+
483 | e0 || e0 | e1 | e2 | e3 |
484 +----++----+----+----+----+
485 | e1 || e1 | e0 | 0 | 0 |
486 +----++----+----+----+----+
487 | e2 || e2 | 0 | e0 | 0 |
488 +----++----+----+----+----+
489 | e3 || e3 | 0 | 0 | e0 |
490 +----++----+----+----+----+
493 M
= list(self
._multiplication
_table
) # copy
494 for i
in xrange(len(M
)):
495 # M had better be "square"
496 M
[i
] = [self
.monomial(i
)] + M
[i
]
497 M
= [["*"] + list(self
.gens())] + M
498 return table(M
, header_row
=True, header_column
=True, frame
=True)
501 def natural_basis(self
):
503 Return a more-natural representation of this algebra's basis.
505 Every finite-dimensional Euclidean Jordan Algebra is a direct
506 sum of five simple algebras, four of which comprise Hermitian
507 matrices. This method returns the original "natural" basis
508 for our underlying vector space. (Typically, the natural basis
509 is used to construct the multiplication table in the first place.)
511 Note that this will always return a matrix. The standard basis
512 in `R^n` will be returned as `n`-by-`1` column matrices.
516 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
517 ....: RealSymmetricEJA)
521 sage: J = RealSymmetricEJA(2)
523 Finite family {0: e0, 1: e1, 2: e2}
524 sage: J.natural_basis()
526 [1 0] [ 0 1/2*sqrt2] [0 0]
527 [0 0], [1/2*sqrt2 0], [0 1]
532 sage: J = JordanSpinEJA(2)
534 Finite family {0: e0, 1: e1}
535 sage: J.natural_basis()
542 if self
._natural
_basis
is None:
543 M
= self
.natural_basis_space()
544 return tuple( M(b
.to_vector()) for b
in self
.basis() )
546 return self
._natural
_basis
549 def natural_basis_space(self
):
551 Return the matrix space in which this algebra's natural basis
554 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
555 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
557 return self
._natural
_basis
[0].matrix_space()
561 def natural_inner_product(X
,Y
):
563 Compute the inner product of two naturally-represented elements.
565 For example in the real symmetric matrix EJA, this will compute
566 the trace inner-product of two n-by-n symmetric matrices. The
567 default should work for the real cartesian product EJA, the
568 Jordan spin EJA, and the real symmetric matrices. The others
569 will have to be overridden.
571 return (X
.conjugate_transpose()*Y
).trace()
577 Return the unit element of this algebra.
581 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
586 sage: J = RealCartesianProductEJA(5)
588 e0 + e1 + e2 + e3 + e4
592 The identity element acts like the identity::
594 sage: set_random_seed()
595 sage: J = random_eja()
596 sage: x = J.random_element()
597 sage: J.one()*x == x and x*J.one() == x
600 The matrix of the unit element's operator is the identity::
602 sage: set_random_seed()
603 sage: J = random_eja()
604 sage: actual = J.one().operator().matrix()
605 sage: expected = matrix.identity(J.base_ring(), J.dimension())
606 sage: actual == expected
610 # We can brute-force compute the matrices of the operators
611 # that correspond to the basis elements of this algebra.
612 # If some linear combination of those basis elements is the
613 # algebra identity, then the same linear combination of
614 # their matrices has to be the identity matrix.
616 # Of course, matrices aren't vectors in sage, so we have to
617 # appeal to the "long vectors" isometry.
618 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
620 # Now we use basis linear algebra to find the coefficients,
621 # of the matrices-as-vectors-linear-combination, which should
622 # work for the original algebra basis too.
623 A
= matrix
.column(self
.base_ring(), oper_vecs
)
625 # We used the isometry on the left-hand side already, but we
626 # still need to do it for the right-hand side. Recall that we
627 # wanted something that summed to the identity matrix.
628 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
630 # Now if there's an identity element in the algebra, this should work.
631 coeffs
= A
.solve_right(b
)
632 return self
.linear_combination(zip(self
.gens(), coeffs
))
635 def random_element(self
):
636 # Temporary workaround for https://trac.sagemath.org/ticket/28327
637 if self
.is_trivial():
640 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
641 return s
.random_element()
643 def random_elements(self
, count
):
645 Return ``count`` random elements as a tuple.
649 sage: from mjo.eja.eja_algebra import JordanSpinEJA
653 sage: J = JordanSpinEJA(3)
654 sage: x,y,z = J.random_elements(3)
655 sage: all( [ x in J, y in J, z in J ])
657 sage: len( J.random_elements(10) ) == 10
661 return tuple( self
.random_element() for idx
in xrange(count
) )
664 def random_instance(cls
, field
=QQ
, **kwargs
):
666 Return a random instance of this type of algebra.
668 In subclasses for algebras that we know how to construct, this
669 is a shortcut for constructing test cases and examples.
671 if cls
is FiniteDimensionalEuclideanJordanAlgebra
:
672 # Red flag! But in theory we could do this I guess. The
673 # only finite-dimensional exceptional EJA is the
674 # octononions. So, we could just create an EJA from an
675 # associative matrix algebra (generated by a subset of
676 # elements) with the symmetric product. Or, we could punt
677 # to random_eja() here, override it in our subclasses, and
678 # not worry about it.
679 raise NotImplementedError
681 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
682 return cls(n
, field
, **kwargs
)
687 Return the rank of this EJA.
691 The author knows of no algorithm to compute the rank of an EJA
692 where only the multiplication table is known. In lieu of one, we
693 require the rank to be specified when the algebra is created,
694 and simply pass along that number here.
698 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
699 ....: RealSymmetricEJA,
700 ....: ComplexHermitianEJA,
701 ....: QuaternionHermitianEJA,
706 The rank of the Jordan spin algebra is always two::
708 sage: JordanSpinEJA(2).rank()
710 sage: JordanSpinEJA(3).rank()
712 sage: JordanSpinEJA(4).rank()
715 The rank of the `n`-by-`n` Hermitian real, complex, or
716 quaternion matrices is `n`::
718 sage: RealSymmetricEJA(4).rank()
720 sage: ComplexHermitianEJA(3).rank()
722 sage: QuaternionHermitianEJA(2).rank()
727 Ensure that every EJA that we know how to construct has a
728 positive integer rank::
730 sage: set_random_seed()
731 sage: r = random_eja().rank()
732 sage: r in ZZ and r > 0
739 def vector_space(self
):
741 Return the vector space that underlies this algebra.
745 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
749 sage: J = RealSymmetricEJA(2)
750 sage: J.vector_space()
751 Vector space of dimension 3 over...
754 return self
.zero().to_vector().parent().ambient_vector_space()
757 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
760 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
):
762 Return the Euclidean Jordan Algebra corresponding to the set
763 `R^n` under the Hadamard product.
765 Note: this is nothing more than the Cartesian product of ``n``
766 copies of the spin algebra. Once Cartesian product algebras
767 are implemented, this can go.
771 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
775 This multiplication table can be verified by hand::
777 sage: J = RealCartesianProductEJA(3)
778 sage: e0,e1,e2 = J.gens()
794 We can change the generator prefix::
796 sage: RealCartesianProductEJA(3, prefix='r').gens()
800 def __init__(self
, n
, field
=QQ
, **kwargs
):
801 V
= VectorSpace(field
, n
)
802 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in xrange(n
) ]
805 fdeja
= super(RealCartesianProductEJA
, self
)
806 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
808 def inner_product(self
, x
, y
):
810 Faster to reimplement than to use natural representations.
814 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
818 Ensure that this is the usual inner product for the algebras
821 sage: set_random_seed()
822 sage: J = RealCartesianProductEJA.random_instance()
823 sage: x,y = J.random_elements(2)
824 sage: X = x.natural_representation()
825 sage: Y = y.natural_representation()
826 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
830 return x
.to_vector().inner_product(y
.to_vector())
835 Return a "random" finite-dimensional Euclidean Jordan Algebra.
839 For now, we choose a random natural number ``n`` (greater than zero)
840 and then give you back one of the following:
842 * The cartesian product of the rational numbers ``n`` times; this is
843 ``QQ^n`` with the Hadamard product.
845 * The Jordan spin algebra on ``QQ^n``.
847 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
850 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
851 in the space of ``2n``-by-``2n`` real symmetric matrices.
853 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
854 in the space of ``4n``-by-``4n`` real symmetric matrices.
856 Later this might be extended to return Cartesian products of the
861 sage: from mjo.eja.eja_algebra import random_eja
866 Euclidean Jordan algebra of dimension...
869 classname
= choice([RealCartesianProductEJA
,
873 QuaternionHermitianEJA
])
874 return classname
.random_instance()
881 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
883 def _max_test_case_size():
884 # Play it safe, since this will be squared and the underlying
885 # field can have dimension 4 (quaternions) too.
888 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
890 Compared to the superclass constructor, we take a basis instead of
891 a multiplication table because the latter can be computed in terms
892 of the former when the product is known (like it is here).
894 # Used in this class's fast _charpoly_coeff() override.
895 self
._basis
_normalizers
= None
897 # We're going to loop through this a few times, so now's a good
898 # time to ensure that it isn't a generator expression.
901 if rank
> 1 and normalize_basis
:
902 # We'll need sqrt(2) to normalize the basis, and this
903 # winds up in the multiplication table, so the whole
904 # algebra needs to be over the field extension.
905 R
= PolynomialRing(field
, 'z')
908 if p
.is_irreducible():
909 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
910 basis
= tuple( s
.change_ring(field
) for s
in basis
)
911 self
._basis
_normalizers
= tuple(
912 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
913 basis
= tuple(s
*c
for (s
,c
) in izip(basis
,self
._basis
_normalizers
))
915 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
917 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
918 return fdeja
.__init
__(field
,
926 def _charpoly_coeff(self
, i
):
928 Override the parent method with something that tries to compute
929 over a faster (non-extension) field.
931 if self
._basis
_normalizers
is None:
932 # We didn't normalize, so assume that the basis we started
933 # with had entries in a nice field.
934 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
936 basis
= ( (b
/n
) for (b
,n
) in izip(self
.natural_basis(),
937 self
._basis
_normalizers
) )
938 field
= self
.base_ring().base_ring() # yeeeaahhhhhhh
939 J
= MatrixEuclideanJordanAlgebra(field
,
942 normalize_basis
=False)
943 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
944 p
= J
._charpoly
_coeff
(i
)
945 # p might be missing some vars, have to substitute "optionally"
946 pairs
= izip(x
.base_ring().gens(), self
._basis
_normalizers
)
947 substitutions
= { v: v*c for (v,c) in pairs }
948 return p
.subs(substitutions
)
952 def multiplication_table_from_matrix_basis(basis
):
954 At least three of the five simple Euclidean Jordan algebras have the
955 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
956 multiplication on the right is matrix multiplication. Given a basis
957 for the underlying matrix space, this function returns a
958 multiplication table (obtained by looping through the basis
959 elements) for an algebra of those matrices.
961 # In S^2, for example, we nominally have four coordinates even
962 # though the space is of dimension three only. The vector space V
963 # is supposed to hold the entire long vector, and the subspace W
964 # of V will be spanned by the vectors that arise from symmetric
965 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
966 field
= basis
[0].base_ring()
967 dimension
= basis
[0].nrows()
969 V
= VectorSpace(field
, dimension
**2)
970 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
972 mult_table
= [[W
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
975 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
976 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
984 Embed the matrix ``M`` into a space of real matrices.
986 The matrix ``M`` can have entries in any field at the moment:
987 the real numbers, complex numbers, or quaternions. And although
988 they are not a field, we can probably support octonions at some
989 point, too. This function returns a real matrix that "acts like"
990 the original with respect to matrix multiplication; i.e.
992 real_embed(M*N) = real_embed(M)*real_embed(N)
995 raise NotImplementedError
1001 The inverse of :meth:`real_embed`.
1003 raise NotImplementedError
1007 def natural_inner_product(cls
,X
,Y
):
1008 Xu
= cls
.real_unembed(X
)
1009 Yu
= cls
.real_unembed(Y
)
1010 tr
= (Xu
*Yu
).trace()
1012 # It's real already.
1015 # Otherwise, try the thing that works for complex numbers; and
1016 # if that doesn't work, the thing that works for quaternions.
1018 return tr
.vector()[0] # real part, imag part is index 1
1019 except AttributeError:
1020 # A quaternions doesn't have a vector() method, but does
1021 # have coefficient_tuple() method that returns the
1022 # coefficients of 1, i, j, and k -- in that order.
1023 return tr
.coefficient_tuple()[0]
1026 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1030 The identity function, for embedding real matrices into real
1036 def real_unembed(M
):
1038 The identity function, for unembedding real matrices from real
1044 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
):
1046 The rank-n simple EJA consisting of real symmetric n-by-n
1047 matrices, the usual symmetric Jordan product, and the trace inner
1048 product. It has dimension `(n^2 + n)/2` over the reals.
1052 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1056 sage: J = RealSymmetricEJA(2)
1057 sage: e0, e1, e2 = J.gens()
1067 The dimension of this algebra is `(n^2 + n) / 2`::
1069 sage: set_random_seed()
1070 sage: n_max = RealSymmetricEJA._max_test_case_size()
1071 sage: n = ZZ.random_element(1, n_max)
1072 sage: J = RealSymmetricEJA(n)
1073 sage: J.dimension() == (n^2 + n)/2
1076 The Jordan multiplication is what we think it is::
1078 sage: set_random_seed()
1079 sage: J = RealSymmetricEJA.random_instance()
1080 sage: x,y = J.random_elements(2)
1081 sage: actual = (x*y).natural_representation()
1082 sage: X = x.natural_representation()
1083 sage: Y = y.natural_representation()
1084 sage: expected = (X*Y + Y*X)/2
1085 sage: actual == expected
1087 sage: J(expected) == x*y
1090 We can change the generator prefix::
1092 sage: RealSymmetricEJA(3, prefix='q').gens()
1093 (q0, q1, q2, q3, q4, q5)
1095 Our natural basis is normalized with respect to the natural inner
1096 product unless we specify otherwise::
1098 sage: set_random_seed()
1099 sage: J = RealSymmetricEJA.random_instance()
1100 sage: all( b.norm() == 1 for b in J.gens() )
1103 Since our natural basis is normalized with respect to the natural
1104 inner product, and since we know that this algebra is an EJA, any
1105 left-multiplication operator's matrix will be symmetric because
1106 natural->EJA basis representation is an isometry and within the EJA
1107 the operator is self-adjoint by the Jordan axiom::
1109 sage: set_random_seed()
1110 sage: x = RealSymmetricEJA.random_instance().random_element()
1111 sage: x.operator().matrix().is_symmetric()
1116 def _denormalized_basis(cls
, n
, field
):
1118 Return a basis for the space of real symmetric n-by-n matrices.
1122 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1126 sage: set_random_seed()
1127 sage: n = ZZ.random_element(1,5)
1128 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1129 sage: all( M.is_symmetric() for M in B)
1133 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1137 for j
in xrange(i
+1):
1138 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1142 Sij
= Eij
+ Eij
.transpose()
1148 def _max_test_case_size():
1149 return 4 # Dimension 10
1152 def __init__(self
, n
, field
=QQ
, **kwargs
):
1153 basis
= self
._denormalized
_basis
(n
, field
)
1154 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1157 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1161 Embed the n-by-n complex matrix ``M`` into the space of real
1162 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1163 bi` to the block matrix ``[[a,b],[-b,a]]``.
1167 sage: from mjo.eja.eja_algebra import \
1168 ....: ComplexMatrixEuclideanJordanAlgebra
1172 sage: F = QuadraticField(-1, 'i')
1173 sage: x1 = F(4 - 2*i)
1174 sage: x2 = F(1 + 2*i)
1177 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1178 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1187 Embedding is a homomorphism (isomorphism, in fact)::
1189 sage: set_random_seed()
1190 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1191 sage: n = ZZ.random_element(n_max)
1192 sage: F = QuadraticField(-1, 'i')
1193 sage: X = random_matrix(F, n)
1194 sage: Y = random_matrix(F, n)
1195 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1196 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1197 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1204 raise ValueError("the matrix 'M' must be square")
1205 field
= M
.base_ring()
1208 a
= z
.vector()[0] # real part, I guess
1209 b
= z
.vector()[1] # imag part, I guess
1210 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1212 # We can drop the imaginaries here.
1213 return matrix
.block(field
.base_ring(), n
, blocks
)
1217 def real_unembed(M
):
1219 The inverse of _embed_complex_matrix().
1223 sage: from mjo.eja.eja_algebra import \
1224 ....: ComplexMatrixEuclideanJordanAlgebra
1228 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1229 ....: [-2, 1, -4, 3],
1230 ....: [ 9, 10, 11, 12],
1231 ....: [-10, 9, -12, 11] ])
1232 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1234 [ 10*i + 9 12*i + 11]
1238 Unembedding is the inverse of embedding::
1240 sage: set_random_seed()
1241 sage: F = QuadraticField(-1, 'i')
1242 sage: M = random_matrix(F, 3)
1243 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1244 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1250 raise ValueError("the matrix 'M' must be square")
1251 if not n
.mod(2).is_zero():
1252 raise ValueError("the matrix 'M' must be a complex embedding")
1254 field
= M
.base_ring() # This should already have sqrt2
1255 R
= PolynomialRing(field
, 'z')
1257 F
= NumberField(z
**2 + 1,'i', embedding
=CLF(-1).sqrt())
1260 # Go top-left to bottom-right (reading order), converting every
1261 # 2-by-2 block we see to a single complex element.
1263 for k
in xrange(n
/2):
1264 for j
in xrange(n
/2):
1265 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1266 if submat
[0,0] != submat
[1,1]:
1267 raise ValueError('bad on-diagonal submatrix')
1268 if submat
[0,1] != -submat
[1,0]:
1269 raise ValueError('bad off-diagonal submatrix')
1270 z
= submat
[0,0] + submat
[0,1]*i
1273 return matrix(F
, n
/2, elements
)
1277 def natural_inner_product(cls
,X
,Y
):
1279 Compute a natural inner product in this algebra directly from
1284 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1288 This gives the same answer as the slow, default method implemented
1289 in :class:`MatrixEuclideanJordanAlgebra`::
1291 sage: set_random_seed()
1292 sage: J = ComplexHermitianEJA.random_instance()
1293 sage: x,y = J.random_elements(2)
1294 sage: Xe = x.natural_representation()
1295 sage: Ye = y.natural_representation()
1296 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1297 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1298 sage: expected = (X*Y).trace().vector()[0]
1299 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1300 sage: actual == expected
1304 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1307 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
):
1309 The rank-n simple EJA consisting of complex Hermitian n-by-n
1310 matrices over the real numbers, the usual symmetric Jordan product,
1311 and the real-part-of-trace inner product. It has dimension `n^2` over
1316 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1320 The dimension of this algebra is `n^2`::
1322 sage: set_random_seed()
1323 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1324 sage: n = ZZ.random_element(1, n_max)
1325 sage: J = ComplexHermitianEJA(n)
1326 sage: J.dimension() == n^2
1329 The Jordan multiplication is what we think it is::
1331 sage: set_random_seed()
1332 sage: J = ComplexHermitianEJA.random_instance()
1333 sage: x,y = J.random_elements(2)
1334 sage: actual = (x*y).natural_representation()
1335 sage: X = x.natural_representation()
1336 sage: Y = y.natural_representation()
1337 sage: expected = (X*Y + Y*X)/2
1338 sage: actual == expected
1340 sage: J(expected) == x*y
1343 We can change the generator prefix::
1345 sage: ComplexHermitianEJA(2, prefix='z').gens()
1348 Our natural basis is normalized with respect to the natural inner
1349 product unless we specify otherwise::
1351 sage: set_random_seed()
1352 sage: J = ComplexHermitianEJA.random_instance()
1353 sage: all( b.norm() == 1 for b in J.gens() )
1356 Since our natural basis is normalized with respect to the natural
1357 inner product, and since we know that this algebra is an EJA, any
1358 left-multiplication operator's matrix will be symmetric because
1359 natural->EJA basis representation is an isometry and within the EJA
1360 the operator is self-adjoint by the Jordan axiom::
1362 sage: set_random_seed()
1363 sage: x = ComplexHermitianEJA.random_instance().random_element()
1364 sage: x.operator().matrix().is_symmetric()
1370 def _denormalized_basis(cls
, n
, field
):
1372 Returns a basis for the space of complex Hermitian n-by-n matrices.
1374 Why do we embed these? Basically, because all of numerical linear
1375 algebra assumes that you're working with vectors consisting of `n`
1376 entries from a field and scalars from the same field. There's no way
1377 to tell SageMath that (for example) the vectors contain complex
1378 numbers, while the scalar field is real.
1382 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1386 sage: set_random_seed()
1387 sage: n = ZZ.random_element(1,5)
1388 sage: field = QuadraticField(2, 'sqrt2')
1389 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1390 sage: all( M.is_symmetric() for M in B)
1394 R
= PolynomialRing(field
, 'z')
1396 F
= NumberField(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1399 # This is like the symmetric case, but we need to be careful:
1401 # * We want conjugate-symmetry, not just symmetry.
1402 # * The diagonal will (as a result) be real.
1406 for j
in xrange(i
+1):
1407 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1409 Sij
= cls
.real_embed(Eij
)
1412 # The second one has a minus because it's conjugated.
1413 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1415 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1418 # Since we embedded these, we can drop back to the "field" that we
1419 # started with instead of the complex extension "F".
1420 return ( s
.change_ring(field
) for s
in S
)
1423 def __init__(self
, n
, field
=QQ
, **kwargs
):
1424 basis
= self
._denormalized
_basis
(n
,field
)
1425 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1428 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1432 Embed the n-by-n quaternion matrix ``M`` into the space of real
1433 matrices of size 4n-by-4n by first sending each quaternion entry `z
1434 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1435 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1440 sage: from mjo.eja.eja_algebra import \
1441 ....: QuaternionMatrixEuclideanJordanAlgebra
1445 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1446 sage: i,j,k = Q.gens()
1447 sage: x = 1 + 2*i + 3*j + 4*k
1448 sage: M = matrix(Q, 1, [[x]])
1449 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1455 Embedding is a homomorphism (isomorphism, in fact)::
1457 sage: set_random_seed()
1458 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1459 sage: n = ZZ.random_element(n_max)
1460 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1461 sage: X = random_matrix(Q, n)
1462 sage: Y = random_matrix(Q, n)
1463 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1464 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1465 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1470 quaternions
= M
.base_ring()
1473 raise ValueError("the matrix 'M' must be square")
1475 F
= QuadraticField(-1, 'i')
1480 t
= z
.coefficient_tuple()
1485 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1486 [-c
+ d
*i
, a
- b
*i
]])
1487 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1488 blocks
.append(realM
)
1490 # We should have real entries by now, so use the realest field
1491 # we've got for the return value.
1492 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1497 def real_unembed(M
):
1499 The inverse of _embed_quaternion_matrix().
1503 sage: from mjo.eja.eja_algebra import \
1504 ....: QuaternionMatrixEuclideanJordanAlgebra
1508 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1509 ....: [-2, 1, -4, 3],
1510 ....: [-3, 4, 1, -2],
1511 ....: [-4, -3, 2, 1]])
1512 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1513 [1 + 2*i + 3*j + 4*k]
1517 Unembedding is the inverse of embedding::
1519 sage: set_random_seed()
1520 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1521 sage: M = random_matrix(Q, 3)
1522 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1523 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1529 raise ValueError("the matrix 'M' must be square")
1530 if not n
.mod(4).is_zero():
1531 raise ValueError("the matrix 'M' must be a complex embedding")
1533 # Use the base ring of the matrix to ensure that its entries can be
1534 # multiplied by elements of the quaternion algebra.
1535 field
= M
.base_ring()
1536 Q
= QuaternionAlgebra(field
,-1,-1)
1539 # Go top-left to bottom-right (reading order), converting every
1540 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1543 for l
in xrange(n
/4):
1544 for m
in xrange(n
/4):
1545 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1546 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1547 if submat
[0,0] != submat
[1,1].conjugate():
1548 raise ValueError('bad on-diagonal submatrix')
1549 if submat
[0,1] != -submat
[1,0].conjugate():
1550 raise ValueError('bad off-diagonal submatrix')
1551 z
= submat
[0,0].vector()[0] # real part
1552 z
+= submat
[0,0].vector()[1]*i
# imag part
1553 z
+= submat
[0,1].vector()[0]*j
# real part
1554 z
+= submat
[0,1].vector()[1]*k
# imag part
1557 return matrix(Q
, n
/4, elements
)
1561 def natural_inner_product(cls
,X
,Y
):
1563 Compute a natural inner product in this algebra directly from
1568 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1572 This gives the same answer as the slow, default method implemented
1573 in :class:`MatrixEuclideanJordanAlgebra`::
1575 sage: set_random_seed()
1576 sage: J = QuaternionHermitianEJA.random_instance()
1577 sage: x,y = J.random_elements(2)
1578 sage: Xe = x.natural_representation()
1579 sage: Ye = y.natural_representation()
1580 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1581 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1582 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1583 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1584 sage: actual == expected
1588 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1591 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
):
1593 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1594 matrices, the usual symmetric Jordan product, and the
1595 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1600 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1604 The dimension of this algebra is `2*n^2 - n`::
1606 sage: set_random_seed()
1607 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1608 sage: n = ZZ.random_element(1, n_max)
1609 sage: J = QuaternionHermitianEJA(n)
1610 sage: J.dimension() == 2*(n^2) - n
1613 The Jordan multiplication is what we think it is::
1615 sage: set_random_seed()
1616 sage: J = QuaternionHermitianEJA.random_instance()
1617 sage: x,y = J.random_elements(2)
1618 sage: actual = (x*y).natural_representation()
1619 sage: X = x.natural_representation()
1620 sage: Y = y.natural_representation()
1621 sage: expected = (X*Y + Y*X)/2
1622 sage: actual == expected
1624 sage: J(expected) == x*y
1627 We can change the generator prefix::
1629 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1630 (a0, a1, a2, a3, a4, a5)
1632 Our natural basis is normalized with respect to the natural inner
1633 product unless we specify otherwise::
1635 sage: set_random_seed()
1636 sage: J = QuaternionHermitianEJA.random_instance()
1637 sage: all( b.norm() == 1 for b in J.gens() )
1640 Since our natural basis is normalized with respect to the natural
1641 inner product, and since we know that this algebra is an EJA, any
1642 left-multiplication operator's matrix will be symmetric because
1643 natural->EJA basis representation is an isometry and within the EJA
1644 the operator is self-adjoint by the Jordan axiom::
1646 sage: set_random_seed()
1647 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1648 sage: x.operator().matrix().is_symmetric()
1653 def _denormalized_basis(cls
, n
, field
):
1655 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1657 Why do we embed these? Basically, because all of numerical
1658 linear algebra assumes that you're working with vectors consisting
1659 of `n` entries from a field and scalars from the same field. There's
1660 no way to tell SageMath that (for example) the vectors contain
1661 complex numbers, while the scalar field is real.
1665 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1669 sage: set_random_seed()
1670 sage: n = ZZ.random_element(1,5)
1671 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1672 sage: all( M.is_symmetric() for M in B )
1676 Q
= QuaternionAlgebra(QQ
,-1,-1)
1679 # This is like the symmetric case, but we need to be careful:
1681 # * We want conjugate-symmetry, not just symmetry.
1682 # * The diagonal will (as a result) be real.
1686 for j
in xrange(i
+1):
1687 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1689 Sij
= cls
.real_embed(Eij
)
1692 # The second, third, and fourth ones have a minus
1693 # because they're conjugated.
1694 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1696 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1698 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1700 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1705 def __init__(self
, n
, field
=QQ
, **kwargs
):
1706 basis
= self
._denormalized
_basis
(n
,field
)
1707 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1710 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
):
1712 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1713 with the usual inner product and jordan product ``x*y =
1714 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1719 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1723 This multiplication table can be verified by hand::
1725 sage: J = JordanSpinEJA(4)
1726 sage: e0,e1,e2,e3 = J.gens()
1742 We can change the generator prefix::
1744 sage: JordanSpinEJA(2, prefix='B').gens()
1748 def __init__(self
, n
, field
=QQ
, **kwargs
):
1749 V
= VectorSpace(field
, n
)
1750 mult_table
= [[V
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
1760 z0
= x
.inner_product(y
)
1761 zbar
= y0
*xbar
+ x0
*ybar
1762 z
= V([z0
] + zbar
.list())
1763 mult_table
[i
][j
] = z
1765 # The rank of the spin algebra is two, unless we're in a
1766 # one-dimensional ambient space (because the rank is bounded by
1767 # the ambient dimension).
1768 fdeja
= super(JordanSpinEJA
, self
)
1769 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1771 def inner_product(self
, x
, y
):
1773 Faster to reimplement than to use natural representations.
1777 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1781 Ensure that this is the usual inner product for the algebras
1784 sage: set_random_seed()
1785 sage: J = JordanSpinEJA.random_instance()
1786 sage: x,y = J.random_elements(2)
1787 sage: X = x.natural_representation()
1788 sage: Y = y.natural_representation()
1789 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1793 return x
.to_vector().inner_product(y
.to_vector())