mjo/eja/eja_utils.py

   1 from sage.functions.other import sqrt
   2 from sage.matrix.constructor import matrix
   3 from sage.modules.free_module_element import vector
   4 from sage.rings.number_field.number_field import NumberField
   5 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
   6 from sage.rings.real_lazy import RLF
   7
   8 def _mat2vec(m):
   9         return vector(m.base_ring(), m.list())
  10
  11 def _vec2mat(v):
  12         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
  13
  14 def gram_schmidt(v):
  15     """
  16     Perform Gram-Schmidt on the list ``v`` which are assumed to be
  17     vectors over the same base ring. Returns a list of orthonormalized
  18     vectors over the smallest extention ring containing the necessary
  19     roots.
  20
  21     SETUP::
  22
  23         sage: from mjo.eja.eja_utils import gram_schmidt
  24
  25     EXAMPLES::
  26
  27         sage: v1 = vector(QQ,(1,2,3))
  28         sage: v2 = vector(QQ,(1,-1,6))
  29         sage: v3 = vector(QQ,(2,1,-1))
  30         sage: v = [v1,v2,v3]
  31         sage: u = gram_schmidt(v)
  32         sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
  33         True
  34         sage: bool(u[0].inner_product(u[1]) == 0)
  35         True
  36         sage: bool(u[0].inner_product(u[2]) == 0)
  37         True
  38         sage: bool(u[1].inner_product(u[2]) == 0)
  39         True
  40
  41     TESTS:
  42
  43     Ensure that zero vectors don't get in the way::
  44
  45         sage: v1 = vector(QQ,(1,2,3))
  46         sage: v2 = vector(QQ,(1,-1,6))
  47         sage: v3 = vector(QQ,(0,0,0))
  48         sage: v = [v1,v2,v3]
  49         sage: len(gram_schmidt(v)) == 2
  50         True
  51
  52     """
  53     def proj(x,y):
  54         return (y.inner_product(x)/x.inner_product(x))*x
  55
  56     v = list(v) # make a copy, don't clobber the input
  57
  58     # Drop all zero vectors before we start.
  59     v = [ v_i for v_i in v if not v_i.is_zero() ]
  60
  61     if len(v) == 0:
  62         # cool
  63         return v
  64
  65     R = v[0].base_ring()
  66
  67     # First orthogonalize...
  68     for i in xrange(1,len(v)):
  69         # Earlier vectors can be made into zero so we have to ignore them.
  70         v[i] -= sum( proj(v[j],v[i]) for j in range(i) if not v[j].is_zero() )
  71
  72     # And now drop all zero vectors again if they were "orthogonalized out."
  73     v = [ v_i for v_i in v if not v_i.is_zero() ]
  74
  75     # Just normalize. If the algebra is missing the roots, we can't add
  76     # them here because then our subalgebra would have a bigger field
  77     # than the superalgebra.
  78     for i in xrange(len(v)):
  79         v[i] = v[i] / v[i].norm()
  80
  81     return v