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
25 from mjo
.eja
.eja_element
import FiniteDimensionalEuclideanJordanAlgebraElement
26 from mjo
.eja
.eja_utils
import _mat2vec
28 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule
):
29 # This is an ugly hack needed to prevent the category framework
30 # from implementing a coercion from our base ring (e.g. the
31 # rationals) into the algebra. First of all -- such a coercion is
32 # nonsense to begin with. But more importantly, it tries to do so
33 # in the category of rings, and since our algebras aren't
34 # associative they generally won't be rings.
35 _no_generic_basering_coercion
= True
47 sage: from mjo.eja.eja_algebra import random_eja
51 By definition, Jordan multiplication commutes::
53 sage: set_random_seed()
54 sage: J = random_eja()
55 sage: x,y = J.random_elements(2)
61 self
._natural
_basis
= natural_basis
64 category
= MagmaticAlgebras(field
).FiniteDimensional()
65 category
= category
.WithBasis().Unital()
67 fda
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
69 range(len(mult_table
)),
72 self
.print_options(bracket
='')
74 # The multiplication table we're given is necessarily in terms
75 # of vectors, because we don't have an algebra yet for
76 # anything to be an element of. However, it's faster in the
77 # long run to have the multiplication table be in terms of
78 # algebra elements. We do this after calling the superclass
79 # constructor so that from_vector() knows what to do.
80 self
._multiplication
_table
= [ map(lambda x
: self
.from_vector(x
), ls
)
81 for ls
in mult_table
]
84 def _element_constructor_(self
, elt
):
86 Construct an element of this algebra from its natural
89 This gets called only after the parent element _call_ method
90 fails to find a coercion for the argument.
94 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
95 ....: RealCartesianProductEJA,
96 ....: RealSymmetricEJA)
100 The identity in `S^n` is converted to the identity in the EJA::
102 sage: J = RealSymmetricEJA(3)
103 sage: I = matrix.identity(QQ,3)
104 sage: J(I) == J.one()
107 This skew-symmetric matrix can't be represented in the EJA::
109 sage: J = RealSymmetricEJA(3)
110 sage: A = matrix(QQ,3, lambda i,j: i-j)
112 Traceback (most recent call last):
114 ArithmeticError: vector is not in free module
118 Ensure that we can convert any element of the two non-matrix
119 simple algebras (whose natural representations are their usual
120 vector representations) back and forth faithfully::
122 sage: set_random_seed()
123 sage: J = RealCartesianProductEJA.random_instance()
124 sage: x = J.random_element()
125 sage: J(x.to_vector().column()) == x
127 sage: J = JordanSpinEJA.random_instance()
128 sage: x = J.random_element()
129 sage: J(x.to_vector().column()) == x
134 # The superclass implementation of random_element()
135 # needs to be able to coerce "0" into the algebra.
138 natural_basis
= self
.natural_basis()
139 basis_space
= natural_basis
[0].matrix_space()
140 if elt
not in basis_space
:
141 raise ValueError("not a naturally-represented algebra element")
143 # Thanks for nothing! Matrix spaces aren't vector spaces in
144 # Sage, so we have to figure out its natural-basis coordinates
145 # ourselves. We use the basis space's ring instead of the
146 # element's ring because the basis space might be an algebraic
147 # closure whereas the base ring of the 3-by-3 identity matrix
148 # could be QQ instead of QQbar.
149 V
= VectorSpace(basis_space
.base_ring(), elt
.nrows()*elt
.ncols())
150 W
= V
.span_of_basis( _mat2vec(s
) for s
in natural_basis
)
151 coords
= W
.coordinate_vector(_mat2vec(elt
))
152 return self
.from_vector(coords
)
157 Return a string representation of ``self``.
161 sage: from mjo.eja.eja_algebra import JordanSpinEJA
165 Ensure that it says what we think it says::
167 sage: JordanSpinEJA(2, field=QQ)
168 Euclidean Jordan algebra of dimension 2 over Rational Field
169 sage: JordanSpinEJA(3, field=RDF)
170 Euclidean Jordan algebra of dimension 3 over Real Double Field
173 fmt
= "Euclidean Jordan algebra of dimension {} over {}"
174 return fmt
.format(self
.dimension(), self
.base_ring())
176 def product_on_basis(self
, i
, j
):
177 return self
._multiplication
_table
[i
][j
]
179 def _a_regular_element(self
):
181 Guess a regular element. Needed to compute the basis for our
182 characteristic polynomial coefficients.
186 sage: from mjo.eja.eja_algebra import random_eja
190 Ensure that this hacky method succeeds for every algebra that we
191 know how to construct::
193 sage: set_random_seed()
194 sage: J = random_eja()
195 sage: J._a_regular_element().is_regular()
200 z
= self
.sum( (i
+1)*gs
[i
] for i
in range(len(gs
)) )
201 if not z
.is_regular():
202 raise ValueError("don't know a regular element")
207 def _charpoly_basis_space(self
):
209 Return the vector space spanned by the basis used in our
210 characteristic polynomial coefficients. This is used not only to
211 compute those coefficients, but also any time we need to
212 evaluate the coefficients (like when we compute the trace or
215 z
= self
._a
_regular
_element
()
216 # Don't use the parent vector space directly here in case this
217 # happens to be a subalgebra. In that case, we would be e.g.
218 # two-dimensional but span_of_basis() would expect three
220 V
= VectorSpace(self
.base_ring(), self
.vector_space().dimension())
221 basis
= [ (z
**k
).to_vector() for k
in range(self
.rank()) ]
222 V1
= V
.span_of_basis( basis
)
223 b
= (V1
.basis() + V1
.complement().basis())
224 return V
.span_of_basis(b
)
229 def _charpoly_coeff(self
, i
):
231 Return the coefficient polynomial "a_{i}" of this algebra's
232 general characteristic polynomial.
234 Having this be a separate cached method lets us compute and
235 store the trace/determinant (a_{r-1} and a_{0} respectively)
236 separate from the entire characteristic polynomial.
238 (A_of_x
, x
, xr
, detA
) = self
._charpoly
_matrix
_system
()
239 R
= A_of_x
.base_ring()
241 # Guaranteed by theory
244 # Danger: the in-place modification is done for performance
245 # reasons (reconstructing a matrix with huge polynomial
246 # entries is slow), but I don't know how cached_method works,
247 # so it's highly possible that we're modifying some global
248 # list variable by reference, here. In other words, you
249 # probably shouldn't call this method twice on the same
250 # algebra, at the same time, in two threads
251 Ai_orig
= A_of_x
.column(i
)
252 A_of_x
.set_column(i
,xr
)
253 numerator
= A_of_x
.det()
254 A_of_x
.set_column(i
,Ai_orig
)
256 # We're relying on the theory here to ensure that each a_i is
257 # indeed back in R, and the added negative signs are to make
258 # the whole charpoly expression sum to zero.
259 return R(-numerator
/detA
)
263 def _charpoly_matrix_system(self
):
265 Compute the matrix whose entries A_ij are polynomials in
266 X1,...,XN, the vector ``x`` of variables X1,...,XN, the vector
267 corresponding to `x^r` and the determinent of the matrix A =
268 [A_ij]. In other words, all of the fixed (cachable) data needed
269 to compute the coefficients of the characteristic polynomial.
274 # Turn my vector space into a module so that "vectors" can
275 # have multivatiate polynomial entries.
276 names
= tuple('X' + str(i
) for i
in range(1,n
+1))
277 R
= PolynomialRing(self
.base_ring(), names
)
279 # Using change_ring() on the parent's vector space doesn't work
280 # here because, in a subalgebra, that vector space has a basis
281 # and change_ring() tries to bring the basis along with it. And
282 # that doesn't work unless the new ring is a PID, which it usually
286 # Now let x = (X1,X2,...,Xn) be the vector whose entries are
290 # And figure out the "left multiplication by x" matrix in
293 monomial_matrices
= [ self
.monomial(i
).operator().matrix()
294 for i
in range(n
) ] # don't recompute these!
296 ek
= self
.monomial(k
).to_vector()
298 sum( x
[i
]*(monomial_matrices
[i
]*ek
)
299 for i
in range(n
) ) )
300 Lx
= matrix
.column(R
, lmbx_cols
)
302 # Now we can compute powers of x "symbolically"
303 x_powers
= [self
.one().to_vector(), x
]
304 for d
in range(2, r
+1):
305 x_powers
.append( Lx
*(x_powers
[-1]) )
307 idmat
= matrix
.identity(R
, n
)
309 W
= self
._charpoly
_basis
_space
()
310 W
= W
.change_ring(R
.fraction_field())
312 # Starting with the standard coordinates x = (X1,X2,...,Xn)
313 # and then converting the entries to W-coordinates allows us
314 # to pass in the standard coordinates to the charpoly and get
315 # back the right answer. Specifically, with x = (X1,X2,...,Xn),
318 # W.coordinates(x^2) eval'd at (standard z-coords)
322 # W-coords of (standard coords of x^2 eval'd at std-coords of z)
324 # We want the middle equivalent thing in our matrix, but use
325 # the first equivalent thing instead so that we can pass in
326 # standard coordinates.
327 x_powers
= [ W
.coordinate_vector(xp
) for xp
in x_powers
]
328 l2
= [idmat
.column(k
-1) for k
in range(r
+1, n
+1)]
329 A_of_x
= matrix
.column(R
, n
, (x_powers
[:r
] + l2
))
330 return (A_of_x
, x
, x_powers
[r
], A_of_x
.det())
334 def characteristic_polynomial(self
):
336 Return a characteristic polynomial that works for all elements
339 The resulting polynomial has `n+1` variables, where `n` is the
340 dimension of this algebra. The first `n` variables correspond to
341 the coordinates of an algebra element: when evaluated at the
342 coordinates of an algebra element with respect to a certain
343 basis, the result is a univariate polynomial (in the one
344 remaining variable ``t``), namely the characteristic polynomial
349 sage: from mjo.eja.eja_algebra import JordanSpinEJA
353 The characteristic polynomial in the spin algebra is given in
354 Alizadeh, Example 11.11::
356 sage: J = JordanSpinEJA(3)
357 sage: p = J.characteristic_polynomial(); p
358 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
359 sage: xvec = J.one().to_vector()
367 # The list of coefficient polynomials a_1, a_2, ..., a_n.
368 a
= [ self
._charpoly
_coeff
(i
) for i
in range(n
) ]
370 # We go to a bit of trouble here to reorder the
371 # indeterminates, so that it's easier to evaluate the
372 # characteristic polynomial at x's coordinates and get back
373 # something in terms of t, which is what we want.
375 S
= PolynomialRing(self
.base_ring(),'t')
377 S
= PolynomialRing(S
, R
.variable_names())
380 # Note: all entries past the rth should be zero. The
381 # coefficient of the highest power (x^r) is 1, but it doesn't
382 # appear in the solution vector which contains coefficients
383 # for the other powers (to make them sum to x^r).
385 a
[r
] = 1 # corresponds to x^r
387 # When the rank is equal to the dimension, trying to
388 # assign a[r] goes out-of-bounds.
389 a
.append(1) # corresponds to x^r
391 return sum( a
[k
]*(t
**k
) for k
in xrange(len(a
)) )
394 def inner_product(self
, x
, y
):
396 The inner product associated with this Euclidean Jordan algebra.
398 Defaults to the trace inner product, but can be overridden by
399 subclasses if they are sure that the necessary properties are
404 sage: from mjo.eja.eja_algebra import random_eja
408 Our inner product is "associative," which means the following for
409 a symmetric bilinear form::
411 sage: set_random_seed()
412 sage: J = random_eja()
413 sage: x,y,z = J.random_elements(3)
414 sage: (x*y).inner_product(z) == y.inner_product(x*z)
418 X
= x
.natural_representation()
419 Y
= y
.natural_representation()
420 return self
.natural_inner_product(X
,Y
)
423 def is_trivial(self
):
425 Return whether or not this algebra is trivial.
427 A trivial algebra contains only the zero element.
431 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
435 sage: J = ComplexHermitianEJA(3)
438 sage: A = J.zero().subalgebra_generated_by()
443 return self
.dimension() == 0
446 def multiplication_table(self
):
448 Return a visual representation of this algebra's multiplication
449 table (on basis elements).
453 sage: from mjo.eja.eja_algebra import JordanSpinEJA
457 sage: J = JordanSpinEJA(4)
458 sage: J.multiplication_table()
459 +----++----+----+----+----+
460 | * || e0 | e1 | e2 | e3 |
461 +====++====+====+====+====+
462 | e0 || e0 | e1 | e2 | e3 |
463 +----++----+----+----+----+
464 | e1 || e1 | e0 | 0 | 0 |
465 +----++----+----+----+----+
466 | e2 || e2 | 0 | e0 | 0 |
467 +----++----+----+----+----+
468 | e3 || e3 | 0 | 0 | e0 |
469 +----++----+----+----+----+
472 M
= list(self
._multiplication
_table
) # copy
473 for i
in xrange(len(M
)):
474 # M had better be "square"
475 M
[i
] = [self
.monomial(i
)] + M
[i
]
476 M
= [["*"] + list(self
.gens())] + M
477 return table(M
, header_row
=True, header_column
=True, frame
=True)
480 def natural_basis(self
):
482 Return a more-natural representation of this algebra's basis.
484 Every finite-dimensional Euclidean Jordan Algebra is a direct
485 sum of five simple algebras, four of which comprise Hermitian
486 matrices. This method returns the original "natural" basis
487 for our underlying vector space. (Typically, the natural basis
488 is used to construct the multiplication table in the first place.)
490 Note that this will always return a matrix. The standard basis
491 in `R^n` will be returned as `n`-by-`1` column matrices.
495 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
496 ....: RealSymmetricEJA)
500 sage: J = RealSymmetricEJA(2)
502 Finite family {0: e0, 1: e1, 2: e2}
503 sage: J.natural_basis()
505 [1 0] [ 0 1/2*sqrt2] [0 0]
506 [0 0], [1/2*sqrt2 0], [0 1]
511 sage: J = JordanSpinEJA(2)
513 Finite family {0: e0, 1: e1}
514 sage: J.natural_basis()
521 if self
._natural
_basis
is None:
522 M
= self
.natural_basis_space()
523 return tuple( M(b
.to_vector()) for b
in self
.basis() )
525 return self
._natural
_basis
528 def natural_basis_space(self
):
530 Return the matrix space in which this algebra's natural basis
533 if self
._natural
_basis
is None or len(self
._natural
_basis
) == 0:
534 return MatrixSpace(self
.base_ring(), self
.dimension(), 1)
536 return self
._natural
_basis
[0].matrix_space()
540 def natural_inner_product(X
,Y
):
542 Compute the inner product of two naturally-represented elements.
544 For example in the real symmetric matrix EJA, this will compute
545 the trace inner-product of two n-by-n symmetric matrices. The
546 default should work for the real cartesian product EJA, the
547 Jordan spin EJA, and the real symmetric matrices. The others
548 will have to be overridden.
550 return (X
.conjugate_transpose()*Y
).trace()
556 Return the unit element of this algebra.
560 sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
565 sage: J = RealCartesianProductEJA(5)
567 e0 + e1 + e2 + e3 + e4
571 The identity element acts like the identity::
573 sage: set_random_seed()
574 sage: J = random_eja()
575 sage: x = J.random_element()
576 sage: J.one()*x == x and x*J.one() == x
579 The matrix of the unit element's operator is the identity::
581 sage: set_random_seed()
582 sage: J = random_eja()
583 sage: actual = J.one().operator().matrix()
584 sage: expected = matrix.identity(J.base_ring(), J.dimension())
585 sage: actual == expected
589 # We can brute-force compute the matrices of the operators
590 # that correspond to the basis elements of this algebra.
591 # If some linear combination of those basis elements is the
592 # algebra identity, then the same linear combination of
593 # their matrices has to be the identity matrix.
595 # Of course, matrices aren't vectors in sage, so we have to
596 # appeal to the "long vectors" isometry.
597 oper_vecs
= [ _mat2vec(g
.operator().matrix()) for g
in self
.gens() ]
599 # Now we use basis linear algebra to find the coefficients,
600 # of the matrices-as-vectors-linear-combination, which should
601 # work for the original algebra basis too.
602 A
= matrix
.column(self
.base_ring(), oper_vecs
)
604 # We used the isometry on the left-hand side already, but we
605 # still need to do it for the right-hand side. Recall that we
606 # wanted something that summed to the identity matrix.
607 b
= _mat2vec( matrix
.identity(self
.base_ring(), self
.dimension()) )
609 # Now if there's an identity element in the algebra, this should work.
610 coeffs
= A
.solve_right(b
)
611 return self
.linear_combination(zip(self
.gens(), coeffs
))
614 def random_element(self
):
615 # Temporary workaround for https://trac.sagemath.org/ticket/28327
616 if self
.is_trivial():
619 s
= super(FiniteDimensionalEuclideanJordanAlgebra
, self
)
620 return s
.random_element()
622 def random_elements(self
, count
):
624 Return ``count`` random elements as a tuple.
628 sage: from mjo.eja.eja_algebra import JordanSpinEJA
632 sage: J = JordanSpinEJA(3)
633 sage: x,y,z = J.random_elements(3)
634 sage: all( [ x in J, y in J, z in J ])
636 sage: len( J.random_elements(10) ) == 10
640 return tuple( self
.random_element() for idx
in xrange(count
) )
645 Return the rank of this EJA.
649 The author knows of no algorithm to compute the rank of an EJA
650 where only the multiplication table is known. In lieu of one, we
651 require the rank to be specified when the algebra is created,
652 and simply pass along that number here.
656 sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
657 ....: RealSymmetricEJA,
658 ....: ComplexHermitianEJA,
659 ....: QuaternionHermitianEJA,
664 The rank of the Jordan spin algebra is always two::
666 sage: JordanSpinEJA(2).rank()
668 sage: JordanSpinEJA(3).rank()
670 sage: JordanSpinEJA(4).rank()
673 The rank of the `n`-by-`n` Hermitian real, complex, or
674 quaternion matrices is `n`::
676 sage: RealSymmetricEJA(4).rank()
678 sage: ComplexHermitianEJA(3).rank()
680 sage: QuaternionHermitianEJA(2).rank()
685 Ensure that every EJA that we know how to construct has a
686 positive integer rank::
688 sage: set_random_seed()
689 sage: r = random_eja().rank()
690 sage: r in ZZ and r > 0
697 def vector_space(self
):
699 Return the vector space that underlies this algebra.
703 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
707 sage: J = RealSymmetricEJA(2)
708 sage: J.vector_space()
709 Vector space of dimension 3 over...
712 return self
.zero().to_vector().parent().ambient_vector_space()
715 Element
= FiniteDimensionalEuclideanJordanAlgebraElement
718 class KnownRankEJA(object):
720 A class for algebras that we actually know we can construct. The
721 main issue is that, for most of our methods to make sense, we need
722 to know the rank of our algebra. Thus we can't simply generate a
723 "random" algebra, or even check that a given basis and product
724 satisfy the axioms; because even if everything looks OK, we wouldn't
725 know the rank we need to actuallty build the thing.
727 Not really a subclass of FDEJA because doing that causes method
728 resolution errors, e.g.
730 TypeError: Error when calling the metaclass bases
731 Cannot create a consistent method resolution
732 order (MRO) for bases FiniteDimensionalEuclideanJordanAlgebra,
737 def _max_test_case_size():
739 Return an integer "size" that is an upper bound on the size of
740 this algebra when it is used in a random test
741 case. Unfortunately, the term "size" is quite vague -- when
742 dealing with `R^n` under either the Hadamard or Jordan spin
743 product, the "size" refers to the dimension `n`. When dealing
744 with a matrix algebra (real symmetric or complex/quaternion
745 Hermitian), it refers to the size of the matrix, which is
746 far less than the dimension of the underlying vector space.
748 We default to five in this class, which is safe in `R^n`. The
749 matrix algebra subclasses (or any class where the "size" is
750 interpreted to be far less than the dimension) should override
751 with a smaller number.
756 def random_instance(cls
, field
=QQ
, **kwargs
):
758 Return a random instance of this type of algebra.
760 Beware, this will crash for "most instances" because the
761 constructor below looks wrong.
763 n
= ZZ
.random_element(cls
._max
_test
_case
_size
()) + 1
764 return cls(n
, field
, **kwargs
)
767 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra
,
770 Return the Euclidean Jordan Algebra corresponding to the set
771 `R^n` under the Hadamard product.
773 Note: this is nothing more than the Cartesian product of ``n``
774 copies of the spin algebra. Once Cartesian product algebras
775 are implemented, this can go.
779 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
783 This multiplication table can be verified by hand::
785 sage: J = RealCartesianProductEJA(3)
786 sage: e0,e1,e2 = J.gens()
802 We can change the generator prefix::
804 sage: RealCartesianProductEJA(3, prefix='r').gens()
808 def __init__(self
, n
, field
=QQ
, **kwargs
):
809 V
= VectorSpace(field
, n
)
810 mult_table
= [ [ V
.gen(i
)*(i
== j
) for j
in xrange(n
) ]
813 fdeja
= super(RealCartesianProductEJA
, self
)
814 return fdeja
.__init
__(field
, mult_table
, rank
=n
, **kwargs
)
816 def inner_product(self
, x
, y
):
818 Faster to reimplement than to use natural representations.
822 sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
826 Ensure that this is the usual inner product for the algebras
829 sage: set_random_seed()
830 sage: J = RealCartesianProductEJA.random_instance()
831 sage: x,y = J.random_elements(2)
832 sage: X = x.natural_representation()
833 sage: Y = y.natural_representation()
834 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
838 return x
.to_vector().inner_product(y
.to_vector())
843 Return a "random" finite-dimensional Euclidean Jordan Algebra.
847 For now, we choose a random natural number ``n`` (greater than zero)
848 and then give you back one of the following:
850 * The cartesian product of the rational numbers ``n`` times; this is
851 ``QQ^n`` with the Hadamard product.
853 * The Jordan spin algebra on ``QQ^n``.
855 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
858 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
859 in the space of ``2n``-by-``2n`` real symmetric matrices.
861 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
862 in the space of ``4n``-by-``4n`` real symmetric matrices.
864 Later this might be extended to return Cartesian products of the
869 sage: from mjo.eja.eja_algebra import random_eja
874 Euclidean Jordan algebra of dimension...
877 classname
= choice(KnownRankEJA
.__subclasses
__())
878 return classname
.random_instance()
885 class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra
):
887 def _max_test_case_size():
888 # Play it safe, since this will be squared and the underlying
889 # field can have dimension 4 (quaternions) too.
892 def __init__(self
, field
, basis
, rank
, normalize_basis
=True, **kwargs
):
894 Compared to the superclass constructor, we take a basis instead of
895 a multiplication table because the latter can be computed in terms
896 of the former when the product is known (like it is here).
898 # Used in this class's fast _charpoly_coeff() override.
899 self
._basis
_normalizers
= None
901 # We're going to loop through this a few times, so now's a good
902 # time to ensure that it isn't a generator expression.
905 if rank
> 1 and normalize_basis
:
906 # We'll need sqrt(2) to normalize the basis, and this
907 # winds up in the multiplication table, so the whole
908 # algebra needs to be over the field extension.
909 R
= PolynomialRing(field
, 'z')
912 if p
.is_irreducible():
913 field
= NumberField(p
, 'sqrt2', embedding
=RLF(2).sqrt())
914 basis
= tuple( s
.change_ring(field
) for s
in basis
)
915 self
._basis
_normalizers
= tuple(
916 ~
(self
.natural_inner_product(s
,s
).sqrt()) for s
in basis
)
917 basis
= tuple(s
*c
for (s
,c
) in izip(basis
,self
._basis
_normalizers
))
919 Qs
= self
.multiplication_table_from_matrix_basis(basis
)
921 fdeja
= super(MatrixEuclideanJordanAlgebra
, self
)
922 return fdeja
.__init
__(field
,
930 def _charpoly_coeff(self
, i
):
932 Override the parent method with something that tries to compute
933 over a faster (non-extension) field.
935 if self
._basis
_normalizers
is None:
936 # We didn't normalize, so assume that the basis we started
937 # with had entries in a nice field.
938 return super(MatrixEuclideanJordanAlgebra
, self
)._charpoly
_coeff
(i
)
940 basis
= ( (b
/n
) for (b
,n
) in izip(self
.natural_basis(),
941 self
._basis
_normalizers
) )
942 field
= self
.base_ring().base_ring() # yeeeaahhhhhhh
943 J
= MatrixEuclideanJordanAlgebra(field
,
946 normalize_basis
=False)
947 (_
,x
,_
,_
) = J
._charpoly
_matrix
_system
()
948 p
= J
._charpoly
_coeff
(i
)
949 # p might be missing some vars, have to substitute "optionally"
950 pairs
= izip(x
.base_ring().gens(), self
._basis
_normalizers
)
951 substitutions
= { v: v*c for (v,c) in pairs }
952 return p
.subs(substitutions
)
956 def multiplication_table_from_matrix_basis(basis
):
958 At least three of the five simple Euclidean Jordan algebras have the
959 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
960 multiplication on the right is matrix multiplication. Given a basis
961 for the underlying matrix space, this function returns a
962 multiplication table (obtained by looping through the basis
963 elements) for an algebra of those matrices.
965 # In S^2, for example, we nominally have four coordinates even
966 # though the space is of dimension three only. The vector space V
967 # is supposed to hold the entire long vector, and the subspace W
968 # of V will be spanned by the vectors that arise from symmetric
969 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
970 field
= basis
[0].base_ring()
971 dimension
= basis
[0].nrows()
973 V
= VectorSpace(field
, dimension
**2)
974 W
= V
.span_of_basis( _mat2vec(s
) for s
in basis
)
976 mult_table
= [[W
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
979 mat_entry
= (basis
[i
]*basis
[j
] + basis
[j
]*basis
[i
])/2
980 mult_table
[i
][j
] = W
.coordinate_vector(_mat2vec(mat_entry
))
988 Embed the matrix ``M`` into a space of real matrices.
990 The matrix ``M`` can have entries in any field at the moment:
991 the real numbers, complex numbers, or quaternions. And although
992 they are not a field, we can probably support octonions at some
993 point, too. This function returns a real matrix that "acts like"
994 the original with respect to matrix multiplication; i.e.
996 real_embed(M*N) = real_embed(M)*real_embed(N)
999 raise NotImplementedError
1003 def real_unembed(M
):
1005 The inverse of :meth:`real_embed`.
1007 raise NotImplementedError
1011 def natural_inner_product(cls
,X
,Y
):
1012 Xu
= cls
.real_unembed(X
)
1013 Yu
= cls
.real_unembed(Y
)
1014 tr
= (Xu
*Yu
).trace()
1016 # It's real already.
1019 # Otherwise, try the thing that works for complex numbers; and
1020 # if that doesn't work, the thing that works for quaternions.
1022 return tr
.vector()[0] # real part, imag part is index 1
1023 except AttributeError:
1024 # A quaternions doesn't have a vector() method, but does
1025 # have coefficient_tuple() method that returns the
1026 # coefficients of 1, i, j, and k -- in that order.
1027 return tr
.coefficient_tuple()[0]
1030 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1034 The identity function, for embedding real matrices into real
1040 def real_unembed(M
):
1042 The identity function, for unembedding real matrices from real
1048 class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1050 The rank-n simple EJA consisting of real symmetric n-by-n
1051 matrices, the usual symmetric Jordan product, and the trace inner
1052 product. It has dimension `(n^2 + n)/2` over the reals.
1056 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1060 sage: J = RealSymmetricEJA(2)
1061 sage: e0, e1, e2 = J.gens()
1071 The dimension of this algebra is `(n^2 + n) / 2`::
1073 sage: set_random_seed()
1074 sage: n_max = RealSymmetricEJA._max_test_case_size()
1075 sage: n = ZZ.random_element(1, n_max)
1076 sage: J = RealSymmetricEJA(n)
1077 sage: J.dimension() == (n^2 + n)/2
1080 The Jordan multiplication is what we think it is::
1082 sage: set_random_seed()
1083 sage: J = RealSymmetricEJA.random_instance()
1084 sage: x,y = J.random_elements(2)
1085 sage: actual = (x*y).natural_representation()
1086 sage: X = x.natural_representation()
1087 sage: Y = y.natural_representation()
1088 sage: expected = (X*Y + Y*X)/2
1089 sage: actual == expected
1091 sage: J(expected) == x*y
1094 We can change the generator prefix::
1096 sage: RealSymmetricEJA(3, prefix='q').gens()
1097 (q0, q1, q2, q3, q4, q5)
1099 Our natural basis is normalized with respect to the natural inner
1100 product unless we specify otherwise::
1102 sage: set_random_seed()
1103 sage: J = RealSymmetricEJA.random_instance()
1104 sage: all( b.norm() == 1 for b in J.gens() )
1107 Since our natural basis is normalized with respect to the natural
1108 inner product, and since we know that this algebra is an EJA, any
1109 left-multiplication operator's matrix will be symmetric because
1110 natural->EJA basis representation is an isometry and within the EJA
1111 the operator is self-adjoint by the Jordan axiom::
1113 sage: set_random_seed()
1114 sage: x = RealSymmetricEJA.random_instance().random_element()
1115 sage: x.operator().matrix().is_symmetric()
1120 def _denormalized_basis(cls
, n
, field
):
1122 Return a basis for the space of real symmetric n-by-n matrices.
1126 sage: from mjo.eja.eja_algebra import RealSymmetricEJA
1130 sage: set_random_seed()
1131 sage: n = ZZ.random_element(1,5)
1132 sage: B = RealSymmetricEJA._denormalized_basis(n,QQ)
1133 sage: all( M.is_symmetric() for M in B)
1137 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1141 for j
in xrange(i
+1):
1142 Eij
= matrix(field
, n
, lambda k
,l
: k
==i
and l
==j
)
1146 Sij
= Eij
+ Eij
.transpose()
1152 def _max_test_case_size():
1153 return 4 # Dimension 10
1156 def __init__(self
, n
, field
=QQ
, **kwargs
):
1157 basis
= self
._denormalized
_basis
(n
, field
)
1158 super(RealSymmetricEJA
, self
).__init
__(field
, basis
, n
, **kwargs
)
1161 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1165 Embed the n-by-n complex matrix ``M`` into the space of real
1166 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1167 bi` to the block matrix ``[[a,b],[-b,a]]``.
1171 sage: from mjo.eja.eja_algebra import \
1172 ....: ComplexMatrixEuclideanJordanAlgebra
1176 sage: F = QuadraticField(-1, 'i')
1177 sage: x1 = F(4 - 2*i)
1178 sage: x2 = F(1 + 2*i)
1181 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1182 sage: ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1191 Embedding is a homomorphism (isomorphism, in fact)::
1193 sage: set_random_seed()
1194 sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_test_case_size()
1195 sage: n = ZZ.random_element(n_max)
1196 sage: F = QuadraticField(-1, 'i')
1197 sage: X = random_matrix(F, n)
1198 sage: Y = random_matrix(F, n)
1199 sage: Xe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X)
1200 sage: Ye = ComplexMatrixEuclideanJordanAlgebra.real_embed(Y)
1201 sage: XYe = ComplexMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1208 raise ValueError("the matrix 'M' must be square")
1209 field
= M
.base_ring()
1212 a
= z
.vector()[0] # real part, I guess
1213 b
= z
.vector()[1] # imag part, I guess
1214 blocks
.append(matrix(field
, 2, [[a
,b
],[-b
,a
]]))
1216 # We can drop the imaginaries here.
1217 return matrix
.block(field
.base_ring(), n
, blocks
)
1221 def real_unembed(M
):
1223 The inverse of _embed_complex_matrix().
1227 sage: from mjo.eja.eja_algebra import \
1228 ....: ComplexMatrixEuclideanJordanAlgebra
1232 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1233 ....: [-2, 1, -4, 3],
1234 ....: [ 9, 10, 11, 12],
1235 ....: [-10, 9, -12, 11] ])
1236 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(A)
1238 [ 10*i + 9 12*i + 11]
1242 Unembedding is the inverse of embedding::
1244 sage: set_random_seed()
1245 sage: F = QuadraticField(-1, 'i')
1246 sage: M = random_matrix(F, 3)
1247 sage: Me = ComplexMatrixEuclideanJordanAlgebra.real_embed(M)
1248 sage: ComplexMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1254 raise ValueError("the matrix 'M' must be square")
1255 if not n
.mod(2).is_zero():
1256 raise ValueError("the matrix 'M' must be a complex embedding")
1259 R
= PolynomialRing(field
, 'z')
1261 F
= NumberField(z
**2 + 1,'i', embedding
=CLF(-1).sqrt())
1264 # Go top-left to bottom-right (reading order), converting every
1265 # 2-by-2 block we see to a single complex element.
1267 for k
in xrange(n
/2):
1268 for j
in xrange(n
/2):
1269 submat
= M
[2*k
:2*k
+2,2*j
:2*j
+2]
1270 if submat
[0,0] != submat
[1,1]:
1271 raise ValueError('bad on-diagonal submatrix')
1272 if submat
[0,1] != -submat
[1,0]:
1273 raise ValueError('bad off-diagonal submatrix')
1274 z
= submat
[0,0] + submat
[0,1]*i
1277 return matrix(F
, n
/2, elements
)
1281 def natural_inner_product(cls
,X
,Y
):
1283 Compute a natural inner product in this algebra directly from
1288 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1292 This gives the same answer as the slow, default method implemented
1293 in :class:`MatrixEuclideanJordanAlgebra`::
1295 sage: set_random_seed()
1296 sage: J = ComplexHermitianEJA.random_instance()
1297 sage: x,y = J.random_elements(2)
1298 sage: Xe = x.natural_representation()
1299 sage: Ye = y.natural_representation()
1300 sage: X = ComplexHermitianEJA.real_unembed(Xe)
1301 sage: Y = ComplexHermitianEJA.real_unembed(Ye)
1302 sage: expected = (X*Y).trace().vector()[0]
1303 sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
1304 sage: actual == expected
1308 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/2
1311 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra
, KnownRankEJA
):
1313 The rank-n simple EJA consisting of complex Hermitian n-by-n
1314 matrices over the real numbers, the usual symmetric Jordan product,
1315 and the real-part-of-trace inner product. It has dimension `n^2` over
1320 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1324 The dimension of this algebra is `n^2`::
1326 sage: set_random_seed()
1327 sage: n_max = ComplexHermitianEJA._max_test_case_size()
1328 sage: n = ZZ.random_element(1, n_max)
1329 sage: J = ComplexHermitianEJA(n)
1330 sage: J.dimension() == n^2
1333 The Jordan multiplication is what we think it is::
1335 sage: set_random_seed()
1336 sage: J = ComplexHermitianEJA.random_instance()
1337 sage: x,y = J.random_elements(2)
1338 sage: actual = (x*y).natural_representation()
1339 sage: X = x.natural_representation()
1340 sage: Y = y.natural_representation()
1341 sage: expected = (X*Y + Y*X)/2
1342 sage: actual == expected
1344 sage: J(expected) == x*y
1347 We can change the generator prefix::
1349 sage: ComplexHermitianEJA(2, prefix='z').gens()
1352 Our natural basis is normalized with respect to the natural inner
1353 product unless we specify otherwise::
1355 sage: set_random_seed()
1356 sage: J = ComplexHermitianEJA.random_instance()
1357 sage: all( b.norm() == 1 for b in J.gens() )
1360 Since our natural basis is normalized with respect to the natural
1361 inner product, and since we know that this algebra is an EJA, any
1362 left-multiplication operator's matrix will be symmetric because
1363 natural->EJA basis representation is an isometry and within the EJA
1364 the operator is self-adjoint by the Jordan axiom::
1366 sage: set_random_seed()
1367 sage: x = ComplexHermitianEJA.random_instance().random_element()
1368 sage: x.operator().matrix().is_symmetric()
1374 def _denormalized_basis(cls
, n
, field
):
1376 Returns a basis for the space of complex Hermitian n-by-n matrices.
1378 Why do we embed these? Basically, because all of numerical linear
1379 algebra assumes that you're working with vectors consisting of `n`
1380 entries from a field and scalars from the same field. There's no way
1381 to tell SageMath that (for example) the vectors contain complex
1382 numbers, while the scalar field is real.
1386 sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
1390 sage: set_random_seed()
1391 sage: n = ZZ.random_element(1,5)
1392 sage: field = QuadraticField(2, 'sqrt2')
1393 sage: B = ComplexHermitianEJA._denormalized_basis(n, field)
1394 sage: all( M.is_symmetric() for M in B)
1398 R
= PolynomialRing(QQ
, 'z')
1400 F
= NumberField(z
**2 + 1, 'I', embedding
=CLF(-1).sqrt())
1403 # This is like the symmetric case, but we need to be careful:
1405 # * We want conjugate-symmetry, not just symmetry.
1406 # * The diagonal will (as a result) be real.
1410 for j
in xrange(i
+1):
1411 Eij
= matrix(F
, n
, lambda k
,l
: k
==i
and l
==j
)
1413 Sij
= cls
.real_embed(Eij
)
1416 # The second one has a minus because it's conjugated.
1417 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1419 Sij_imag
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1422 # Since we embedded these, we can drop back to the "field" that we
1423 # started with instead of the complex extension "F".
1424 return ( s
.change_ring(field
) for s
in S
)
1427 def __init__(self
, n
, field
=QQ
, **kwargs
):
1428 basis
= self
._denormalized
_basis
(n
,field
)
1429 super(ComplexHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1432 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra
):
1436 Embed the n-by-n quaternion matrix ``M`` into the space of real
1437 matrices of size 4n-by-4n by first sending each quaternion entry `z
1438 = a + bi + cj + dk` to the block-complex matrix ``[[a + bi,
1439 c+di],[-c + di, a-bi]]`, and then embedding those into a real
1444 sage: from mjo.eja.eja_algebra import \
1445 ....: QuaternionMatrixEuclideanJordanAlgebra
1449 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1450 sage: i,j,k = Q.gens()
1451 sage: x = 1 + 2*i + 3*j + 4*k
1452 sage: M = matrix(Q, 1, [[x]])
1453 sage: QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1459 Embedding is a homomorphism (isomorphism, in fact)::
1461 sage: set_random_seed()
1462 sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_test_case_size()
1463 sage: n = ZZ.random_element(n_max)
1464 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1465 sage: X = random_matrix(Q, n)
1466 sage: Y = random_matrix(Q, n)
1467 sage: Xe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X)
1468 sage: Ye = QuaternionMatrixEuclideanJordanAlgebra.real_embed(Y)
1469 sage: XYe = QuaternionMatrixEuclideanJordanAlgebra.real_embed(X*Y)
1474 quaternions
= M
.base_ring()
1477 raise ValueError("the matrix 'M' must be square")
1479 F
= QuadraticField(-1, 'i')
1484 t
= z
.coefficient_tuple()
1489 cplxM
= matrix(F
, 2, [[ a
+ b
*i
, c
+ d
*i
],
1490 [-c
+ d
*i
, a
- b
*i
]])
1491 realM
= ComplexMatrixEuclideanJordanAlgebra
.real_embed(cplxM
)
1492 blocks
.append(realM
)
1494 # We should have real entries by now, so use the realest field
1495 # we've got for the return value.
1496 return matrix
.block(quaternions
.base_ring(), n
, blocks
)
1501 def real_unembed(M
):
1503 The inverse of _embed_quaternion_matrix().
1507 sage: from mjo.eja.eja_algebra import \
1508 ....: QuaternionMatrixEuclideanJordanAlgebra
1512 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1513 ....: [-2, 1, -4, 3],
1514 ....: [-3, 4, 1, -2],
1515 ....: [-4, -3, 2, 1]])
1516 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(M)
1517 [1 + 2*i + 3*j + 4*k]
1521 Unembedding is the inverse of embedding::
1523 sage: set_random_seed()
1524 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1525 sage: M = random_matrix(Q, 3)
1526 sage: Me = QuaternionMatrixEuclideanJordanAlgebra.real_embed(M)
1527 sage: QuaternionMatrixEuclideanJordanAlgebra.real_unembed(Me) == M
1533 raise ValueError("the matrix 'M' must be square")
1534 if not n
.mod(4).is_zero():
1535 raise ValueError("the matrix 'M' must be a complex embedding")
1537 # Use the base ring of the matrix to ensure that its entries can be
1538 # multiplied by elements of the quaternion algebra.
1539 field
= M
.base_ring()
1540 Q
= QuaternionAlgebra(field
,-1,-1)
1543 # Go top-left to bottom-right (reading order), converting every
1544 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1547 for l
in xrange(n
/4):
1548 for m
in xrange(n
/4):
1549 submat
= ComplexMatrixEuclideanJordanAlgebra
.real_unembed(
1550 M
[4*l
:4*l
+4,4*m
:4*m
+4] )
1551 if submat
[0,0] != submat
[1,1].conjugate():
1552 raise ValueError('bad on-diagonal submatrix')
1553 if submat
[0,1] != -submat
[1,0].conjugate():
1554 raise ValueError('bad off-diagonal submatrix')
1555 z
= submat
[0,0].vector()[0] # real part
1556 z
+= submat
[0,0].vector()[1]*i
# imag part
1557 z
+= submat
[0,1].vector()[0]*j
# real part
1558 z
+= submat
[0,1].vector()[1]*k
# imag part
1561 return matrix(Q
, n
/4, elements
)
1565 def natural_inner_product(cls
,X
,Y
):
1567 Compute a natural inner product in this algebra directly from
1572 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1576 This gives the same answer as the slow, default method implemented
1577 in :class:`MatrixEuclideanJordanAlgebra`::
1579 sage: set_random_seed()
1580 sage: J = QuaternionHermitianEJA.random_instance()
1581 sage: x,y = J.random_elements(2)
1582 sage: Xe = x.natural_representation()
1583 sage: Ye = y.natural_representation()
1584 sage: X = QuaternionHermitianEJA.real_unembed(Xe)
1585 sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
1586 sage: expected = (X*Y).trace().coefficient_tuple()[0]
1587 sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
1588 sage: actual == expected
1592 return RealMatrixEuclideanJordanAlgebra
.natural_inner_product(X
,Y
)/4
1595 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra
,
1598 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1599 matrices, the usual symmetric Jordan product, and the
1600 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1605 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1609 The dimension of this algebra is `2*n^2 - n`::
1611 sage: set_random_seed()
1612 sage: n_max = QuaternionHermitianEJA._max_test_case_size()
1613 sage: n = ZZ.random_element(1, n_max)
1614 sage: J = QuaternionHermitianEJA(n)
1615 sage: J.dimension() == 2*(n^2) - n
1618 The Jordan multiplication is what we think it is::
1620 sage: set_random_seed()
1621 sage: J = QuaternionHermitianEJA.random_instance()
1622 sage: x,y = J.random_elements(2)
1623 sage: actual = (x*y).natural_representation()
1624 sage: X = x.natural_representation()
1625 sage: Y = y.natural_representation()
1626 sage: expected = (X*Y + Y*X)/2
1627 sage: actual == expected
1629 sage: J(expected) == x*y
1632 We can change the generator prefix::
1634 sage: QuaternionHermitianEJA(2, prefix='a').gens()
1635 (a0, a1, a2, a3, a4, a5)
1637 Our natural basis is normalized with respect to the natural inner
1638 product unless we specify otherwise::
1640 sage: set_random_seed()
1641 sage: J = QuaternionHermitianEJA.random_instance()
1642 sage: all( b.norm() == 1 for b in J.gens() )
1645 Since our natural basis is normalized with respect to the natural
1646 inner product, and since we know that this algebra is an EJA, any
1647 left-multiplication operator's matrix will be symmetric because
1648 natural->EJA basis representation is an isometry and within the EJA
1649 the operator is self-adjoint by the Jordan axiom::
1651 sage: set_random_seed()
1652 sage: x = QuaternionHermitianEJA.random_instance().random_element()
1653 sage: x.operator().matrix().is_symmetric()
1658 def _denormalized_basis(cls
, n
, field
):
1660 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1662 Why do we embed these? Basically, because all of numerical
1663 linear algebra assumes that you're working with vectors consisting
1664 of `n` entries from a field and scalars from the same field. There's
1665 no way to tell SageMath that (for example) the vectors contain
1666 complex numbers, while the scalar field is real.
1670 sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
1674 sage: set_random_seed()
1675 sage: n = ZZ.random_element(1,5)
1676 sage: B = QuaternionHermitianEJA._denormalized_basis(n,QQ)
1677 sage: all( M.is_symmetric() for M in B )
1681 Q
= QuaternionAlgebra(QQ
,-1,-1)
1684 # This is like the symmetric case, but we need to be careful:
1686 # * We want conjugate-symmetry, not just symmetry.
1687 # * The diagonal will (as a result) be real.
1691 for j
in xrange(i
+1):
1692 Eij
= matrix(Q
, n
, lambda k
,l
: k
==i
and l
==j
)
1694 Sij
= cls
.real_embed(Eij
)
1697 # The second, third, and fourth ones have a minus
1698 # because they're conjugated.
1699 Sij_real
= cls
.real_embed(Eij
+ Eij
.transpose())
1701 Sij_I
= cls
.real_embed(I
*Eij
- I
*Eij
.transpose())
1703 Sij_J
= cls
.real_embed(J
*Eij
- J
*Eij
.transpose())
1705 Sij_K
= cls
.real_embed(K
*Eij
- K
*Eij
.transpose())
1708 # Since we embedded these, we can drop back to the "field" that we
1709 # started with instead of the quaternion algebra "Q".
1710 return ( s
.change_ring(field
) for s
in S
)
1713 def __init__(self
, n
, field
=QQ
, **kwargs
):
1714 basis
= self
._denormalized
_basis
(n
,field
)
1715 super(QuaternionHermitianEJA
,self
).__init
__(field
, basis
, n
, **kwargs
)
1718 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra
, KnownRankEJA
):
1720 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1721 with the usual inner product and jordan product ``x*y =
1722 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1727 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1731 This multiplication table can be verified by hand::
1733 sage: J = JordanSpinEJA(4)
1734 sage: e0,e1,e2,e3 = J.gens()
1750 We can change the generator prefix::
1752 sage: JordanSpinEJA(2, prefix='B').gens()
1756 def __init__(self
, n
, field
=QQ
, **kwargs
):
1757 V
= VectorSpace(field
, n
)
1758 mult_table
= [[V
.zero() for j
in xrange(n
)] for i
in xrange(n
)]
1768 z0
= x
.inner_product(y
)
1769 zbar
= y0
*xbar
+ x0
*ybar
1770 z
= V([z0
] + zbar
.list())
1771 mult_table
[i
][j
] = z
1773 # The rank of the spin algebra is two, unless we're in a
1774 # one-dimensional ambient space (because the rank is bounded by
1775 # the ambient dimension).
1776 fdeja
= super(JordanSpinEJA
, self
)
1777 return fdeja
.__init
__(field
, mult_table
, rank
=min(n
,2), **kwargs
)
1779 def inner_product(self
, x
, y
):
1781 Faster to reimplement than to use natural representations.
1785 sage: from mjo.eja.eja_algebra import JordanSpinEJA
1789 Ensure that this is the usual inner product for the algebras
1792 sage: set_random_seed()
1793 sage: J = JordanSpinEJA.random_instance()
1794 sage: x,y = J.random_elements(2)
1795 sage: X = x.natural_representation()
1796 sage: Y = y.natural_representation()
1797 sage: x.inner_product(y) == J.natural_inner_product(X,Y)
1801 return x
.to_vector().inner_product(y
.to_vector())