1 from sage
.matrix
.constructor
import matrix
2 from sage
.rings
.all
import QQ
4 from mjo
.eja
.eja_subalgebra
import FiniteDimensionalEuclideanJordanSubalgebra
7 class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclideanJordanSubalgebra
):
8 def __init__(self
, elt
, orthonormalize_basis
):
9 self
._superalgebra
= elt
.parent()
10 category
= self
._superalgebra
.category().Associative()
11 V
= self
._superalgebra
.vector_space()
12 field
= self
._superalgebra
.base_ring()
14 # This list is guaranteed to contain all independent powers,
15 # because it's the maximal set of powers that could possibly
16 # be independent (by a dimension argument).
17 powers
= [ elt
**k
for k
in range(V
.dimension()) ]
18 power_vectors
= [ p
.to_vector() for p
in powers
]
19 P
= matrix(field
, power_vectors
)
21 if orthonormalize_basis
== False:
22 # Echelonize the matrix ourselves, because otherwise the
23 # call to P.pivot_rows() below can choose a non-optimal
24 # row-reduction algorithm. In particular, scaling can
25 # help over AA because it avoids the RecursionError that
26 # gets thrown when we have to look too hard for a root.
28 # Beware: QQ supports an entirely different set of "algorithm"
29 # keywords than do AA and RR.
32 algo
= "scaled_partial_pivoting"
33 P
.echelonize(algorithm
=algo
)
35 # In this case, we just need to figure out which elements
36 # of the "powers" list are redundant... First compute the
37 # vector subspace spanned by the powers of the given
40 # Figure out which powers form a linearly-independent set.
41 ind_rows
= P
.pivot_rows()
43 # Pick those out of the list of all powers.
44 superalgebra_basis
= tuple(map(powers
.__getitem
__, ind_rows
))
46 # If we're going to orthonormalize the basis anyway, we
47 # might as well just do Gram-Schmidt on the whole list of
48 # powers. The redundant ones will get zero'd out. If this
49 # looks like a roundabout way to orthonormalize, it is.
50 # But converting everything from algebra elements to vectors
51 # to matrices and then back again turns out to be about
52 # as fast as reimplementing our own Gram-Schmidt that
54 G
,_
= P
.gram_schmidt(orthonormal
=True)
55 basis_vectors
= [ g
for g
in G
.rows() if not g
.is_zero() ]
56 superalgebra_basis
= [ self
._superalgebra
.from_vector(b
)
57 for b
in basis_vectors
]
59 fdeja
= super(FiniteDimensionalEuclideanJordanElementSubalgebra
, self
)
60 fdeja
.__init
__(self
._superalgebra
,
65 # The rank is the highest possible degree of a minimal
66 # polynomial, and is bounded above by the dimension. We know
67 # in this case that there's an element whose minimal
68 # polynomial has the same degree as the space's dimension
69 # (remember how we constructed the space?), so that must be
71 self
.rank
.set_cache(self
.dimension())
76 Return the multiplicative identity element of this algebra.
78 The superclass method computes the identity element, which is
79 beyond overkill in this case: the superalgebra identity
80 restricted to this algebra is its identity. Note that we can't
81 count on the first basis element being the identity -- it migth
82 have been scaled if we orthonormalized the basis.
86 sage: from mjo.eja.eja_algebra import (HadamardEJA,
91 sage: J = HadamardEJA(5)
93 e0 + e1 + e2 + e3 + e4
94 sage: x = sum(J.gens())
95 sage: A = x.subalgebra_generated_by()
98 sage: A.one().superalgebra_element()
99 e0 + e1 + e2 + e3 + e4
103 The identity element acts like the identity over the rationals::
105 sage: set_random_seed()
106 sage: x = random_eja(field=QQ).random_element()
107 sage: A = x.subalgebra_generated_by()
108 sage: x = A.random_element()
109 sage: A.one()*x == x and x*A.one() == x
112 The identity element acts like the identity over the algebraic
113 reals with an orthonormal basis::
115 sage: set_random_seed()
116 sage: x = random_eja().random_element()
117 sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
118 sage: x = A.random_element()
119 sage: A.one()*x == x and x*A.one() == x
122 The matrix of the unit element's operator is the identity over
125 sage: set_random_seed()
126 sage: x = random_eja(field=QQ).random_element()
127 sage: A = x.subalgebra_generated_by()
128 sage: actual = A.one().operator().matrix()
129 sage: expected = matrix.identity(A.base_ring(), A.dimension())
130 sage: actual == expected
133 The matrix of the unit element's operator is the identity over
134 the algebraic reals with an orthonormal basis::
136 sage: set_random_seed()
137 sage: x = random_eja().random_element()
138 sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
139 sage: actual = A.one().operator().matrix()
140 sage: expected = matrix.identity(A.base_ring(), A.dimension())
141 sage: actual == expected
145 if self
.dimension() == 0:
148 sa_one
= self
.superalgebra().one().to_vector()
149 # The extra hackery is because foo.to_vector() might not
150 # live in foo.parent().vector_space()!
151 coords
= sum( a
*b
for (a
,b
)
153 self
.superalgebra().vector_space().basis()) )
154 return self
.from_vector(self
.vector_space().coordinate_vector(coords
))