-from sage.functions.other import sqrt
-from sage.matrix.constructor import matrix
-from sage.modules.free_module_element import vector
-
-def _charpoly_sage_input(s):
- r"""
- Helper function that you can use on the string output from sage
- to convert a charpoly coefficient into the corresponding input
- to be cached.
-
- SETUP::
-
- sage: from mjo.eja.eja_utils import _charpoly_sage_input
-
- EXAMPLES::
-
- sage: J = JordanSpinEJA(4,QQ)
- sage: J._charpoly_coefficients()[0]
- X1^2 - X2^2 - X3^2 - X4^2
- sage: _charpoly_sage_input("X1^2 - X2^2 - X3^2 - X4^2")
- 'X[0]**2 - X[1]**2 - X[2]**2 - X[3]**2'
-
- """
- import re
-
- exponent_out = r"\^"
- exponent_in = r"**"
-
- digit_out = r"X([0-9]+)"
-
- def replace_digit(m):
- # m is a match object
- return "X[" + str(int(m.group(1)) - 1) + "]"
-
- s = re.sub(exponent_out, exponent_in, s)
- return re.sub(digit_out, replace_digit, s)
-
+from sage.structure.element import is_Matrix
def _scale(x, alpha):
r"""
# first needing to convert them to a list of octonions and
# then recursing down into the list. It also avoids the wonky
# list(x) when x is an element of a CFM. I don't know what it
- # returns but it aint the coordinates. This will fall through
- # to the iterable case the next time around.
- return _all2list(x.to_vector())
+ # returns but it aint the coordinates. We don't recurse
+ # because vectors can only contain ring elements as entries.
+ return x.to_vector().list()
+
+ if is_Matrix(x):
+ # This sucks, but for performance reasons we don't want to
+ # call _all2list recursively on the contents of a matrix
+ # when we don't have to (they only contain ring elements
+ # as entries)
+ return x.list()
try:
xl = list(x)
# Avoid the retardation of list(QQ(1)) == [1].
return [x]
- return sum(list( map(_all2list, xl) ), [])
-
-
+ return sum( map(_all2list, xl) , [])
-def _mat2vec(m):
- return vector(m.base_ring(), m.list())
-
-def _vec2mat(v):
- return matrix(v.base_ring(), sqrt(v.degree()), v.list())
def gram_schmidt(v, inner_product=None):
"""
sage: v = [v1,v2,v3]
sage: len(gram_schmidt(v)) == 2
True
- """
- if inner_product is None:
- inner_product = lambda x,y: x.inner_product(y)
- def norm(x):
- ip = inner_product(x,x)
- # Don't expand the given field; the inner-product's codomain
- # is already correct. For example QQ(2).sqrt() returns sqrt(2)
- # in SR, and that will give you weird errors about symbolics
- # when what's really going wrong is that you're trying to
- # orthonormalize in QQ.
- return ip.parent()(ip.sqrt())
-
- v = list(v) # make a copy, don't clobber the input
-
- # Drop all zero vectors before we start.
- v = [ v_i for v_i in v if not v_i.is_zero() ]
+ """
if len(v) == 0:
# cool
return v
- R = v[0].base_ring()
+ V = v[0].parent()
- # Our "zero" needs to belong to the right space for sum() to work.
- zero = v[0].parent().zero()
+ if inner_product is None:
+ inner_product = lambda x,y: x.inner_product(y)
sc = lambda x,a: a*x
- if hasattr(v[0], 'cartesian_factors'):
+ if hasattr(V, 'cartesian_factors'):
# Only use the slow implementation if necessary.
sc = _scale
def proj(x,y):
+ # project y onto the span of {x}
return sc(x, (inner_product(x,y)/inner_product(x,x)))
- # First orthogonalize...
- for i in range(1,len(v)):
- # Earlier vectors can be made into zero so we have to ignore them.
- v[i] -= sum( (proj(v[j],v[i])
- for j in range(i)
- if not v[j].is_zero() ),
- zero )
+ def normalize(x):
+ # Don't extend the given field with the necessary
+ # square roots. This will probably throw weird
+ # errors about the symbolic ring if you e.g. try
+ # to use it on a set of rational vectors that isn't
+ # already orthonormalized.
+ return sc(x, ~inner_product(x,x).sqrt())
- # And now drop all zero vectors again if they were "orthogonalized out."
- v = [ v_i for v_i in v if not v_i.is_zero() ]
+ v_out = [] # make a copy, don't clobber the input
- # Just normalize. If the algebra is missing the roots, we can't add
- # them here because then our subalgebra would have a bigger field
- # than the superalgebra.
- for i in range(len(v)):
- v[i] = sc(v[i], ~norm(v[i]))
+ for (i, v_i) in enumerate(v):
+ ortho_v_i = v_i - V.sum( proj(v_out[j],v_i) for j in range(i) )
+ if not ortho_v_i.is_zero():
+ v_out.append(normalize(ortho_v_i))
- return v
+ return v_out
projects / sage.d.git / blobdiff