Add features to quinary_module_l

quinary_module_l.sage

+#load('heckes/quinary2.sage')
+#load('ext.sage')
+
+
 ###############################################################################
 ###############################################################################
 ########################### Dimension formulas ################################
@@ -208,6 +212,29 @@ def qf_find_primitive_zeros_mod_p_regular(Q, p):
         v = i*k*v1 + (-Q(v5)-i*j)*v2 + k^2*v3 + j*k*v4 + k*v5
         yield primitivize(v, p)
 
+
+def qf_find_primitive_zeros_mod_p_regular_i(Q, p, B, i):
+
+  """
+  Given a quinary regular quadratic form over ZZ, and an odd prime p,
+  returns the primitive zeros mod p of Q.
+  """
+
+  v1, v2, v3, v4, v5 = B
+  if i == 0:
+    yield primitivize(v1, p)
+  yield primitivize(i*v1+v2, p)
+  yield primitivize(i*v1+v3, p)
+  for j in range(p):
+    yield primitivize((-i*j)*v1+i*v2+j*v3+v4, p)
+    if j != 0:
+      v = -Q(v5)*v1+i*j*v2+j^2*v4+j*v5
+      yield primitivize(v, p)
+    for k in range(1, p):
+      v = i*k*v1 + (-Q(v5)-i*j)*v2 + k^2*v3 + j*k*v4 + k*v5
+      yield primitivize(v, p)
+
+
 def qf_find_primitive_zeros_mod_p_nonsingular(Q, p):
 
   """
@@ -268,6 +295,76 @@ def qf_find_p_neighbor_from_vec(Q, p, v, return_matrix = False):
   good_basis[1]  = good_basis[1] * p                          ## Multiply w by p
   ## Return the associated quadratic form
   M = matrix(good_basis)
+  #new_Q = deepcopy(Q)                        ## Note: This avoids a circular import of QuadraticForm!
+  #if hasattr(new_Q, '__theta_vec'):
+  #    delattr(new_Q, '__theta_vec')
+  #if hasattr(new_Q, '__automorphisms_pari'):
+  #    delattr(new_Q, '__automorphisms_pari')
+  #if hasattr(new_Q, '__number_of_automorphisms'):
+  #    delattr(new_Q, '__number_of_automorphisms')
+  new_Q = QuadraticForm(R, M*Q.matrix()*M.transpose())
+  #new_Q.__init__(R, M * Q.matrix() * M.transpose())
+  if return_matrix:
+      return new_Q, M.transpose()
+  else:
+      return new_Q
+
+
+#está fallando algo acá, el good basis no tiene determinante 1
+def qf_find_p_neighbor_from_vec2(Q, p, v, return_matrix = False):
+
+  """
+  Given a quinary quadratic form over ZZ, a prime p, and a vector
+  v such that Q(v) = 0 (mod p^2), and v != 0 mod p, returns the
+  p-neighbor associated to v, with radical at some place v0.
+  Copied from the qf package of sage, but with bugs fixed.
+  """
+
+  v0 = Q.matrix().change_ring(GF(p)).kernel().0.change_ring(ZZ)
+  M0 = matrix(ZZ, extend_to_primitive(v0))
+  Q = Q(M0.transpose())
+  v0 = vector((1, 0, 0, 0, 0))
+  v = (M0.inverse().change_ring(ZZ)).transpose()*v
+  
+
+  R = Q.base_ring()
+  n = Q.dim()
+  B2 = Q.matrix()
+
+  ## Find a (dual) vector w with B(v,w) != 0 (mod p)
+  v_dual = B2 * vector(v)     ## We want the dot product with this to not be divisible by 2*p.
+  y_ind = 0
+  while ((y_ind < n) and (v_dual[y_ind] % p) == 0):   ## Check the dot product for the std basis vectors!
+      y_ind += 1
+  if y_ind == n:
+      raise RuntimeError("Oops!  One of the standard basis vectors should have worked.")
+  w = vector([R(i == y_ind)  for i in range(n)])
+  vw_prod = (v * Q.matrix()).dot_product(w)
+  s = Q(v)
+  if (s % p**2 != 0):
+      al = (-s / (p * vw_prod)) % p
+      v1 = v + p * al * w
+      v1w_prod = (v1 * Q.matrix()).dot_product(w)
+  else:
+      v1 = v
+      v1w_prod = vw_prod
+  print(Q.bilinear_map(v0, v1)%p, Q.bilinear_map(v0, w)%p, Q.bilinear_map(v1, w)%p)
+  print(matrix([v0, v1, w]))
+  while gcd(matrix([v0, v1, w]).minors(3)) != 1:
+    v1[randint(0, 4)]+=(-1)^randint(0,1)*p
+
+  good_basis = extend_to_primitive([v0, v1, w])
+  print(matrix(good_basis).det())
+  for i in range(3,n):
+      ith_prod = (good_basis[i] * Q.matrix()).dot_product(v1)
+      c = (ith_prod / v1w_prod) % p
+      good_basis[i] = good_basis[i] - c * w  ## Ensures that this extension has <v_i, v> = 0 (mod p)
+  ## Perform the p-neighbor switch
+  good_basis_old = deepcopy(good_basis)
+  good_basis[1]  = vector([x/p  for x in good_basis[1]])    ## Divide v1 by p
+  good_basis[2]  = good_basis[2] * p                          ## Multiply w by p
+  ## Return the associated quadratic form
+  M = matrix(good_basis)
   new_Q = deepcopy(Q)                        ## Note: This avoids a circular import of QuadraticForm!
   if hasattr(new_Q, '__theta_vec'):
       delattr(new_Q, '__theta_vec')
@@ -277,11 +374,12 @@ def qf_find_p_neighbor_from_vec(Q, p, v, return_matrix = False):
       delattr(new_Q, '__number_of_automorphisms')
   new_Q.__init__(R, M * Q.matrix() * M.transpose())
   if return_matrix:
-      return new_Q, M.transpose()
+      return new_Q, (M*M0).transpose()
   else:
       return new_Q
   return QuadraticForm(R, M * Q.matrix() * M.transpose())
 
+
 def qf_find_p_neighbors(Q, p, return_matrix = False):
 
   """
@@ -298,6 +396,23 @@ def qf_find_p_neighbors(Q, p, return_matrix = False):
       Q1 = qf_find_p_neighbor_from_vec(Q, p, v, return_matrix = return_matrix)
       yield qf_reduction_pari(Q1)
 
+def qf_find_p_neighbors_i(Q, p, B, i, return_matrix = False):
+
+  """
+  Given a quinary quadratic form over ZZ, and a prime p,
+  returns the p-neighbors of Q.
+  """
+
+  for v in qf_find_primitive_zeros_mod_p_regular_i(Q, p, B, i):
+    if return_matrix:
+      Q1, M1 = qf_find_p_neighbor_from_vec(Q, p, v, return_matrix = return_matrix)
+      Q2, M2 = qf_reduction_pari(Q1, return_matrix = 1)
+      yield Q2, M1*M2
+    else:
+      Q1 = qf_find_p_neighbor_from_vec(Q, p, v, return_matrix = return_matrix)
+      yield qf_reduction_pari(Q1)
+
+
 def qf_find_genus(Q, p):
 
   """
@@ -561,7 +676,50 @@ def qf_symmetry(Q, v):
   of the symmetry of v.
   """
   
-  return identity_matrix(5) - v.column()*matrix(v)*Q.matrix()/Q.change_ring(QQ)(v)
+  n = Q.dim()
+  return identity_matrix(n) - v.column()*matrix(v)*Q.matrix()/Q.change_ring(QQ)(v)
+
+def qf_automorphism_symmetries(Q, A):
+
+  n = Q.dim()
+  if A == identity_matrix(QQ, n):
+    return []
+  bs = (A - 1).columns()
+  for b1 in bs:
+    if b1 != 0:
+      break
+  A1 = qf_symmetry(Q, b1)*A
+  if A1 == identity_matrix(QQ, n):
+    return [b1]
+  bs = (A1 - 1).columns()
+  for b2 in bs:
+    if b2 != 0:
+      break
+  A1 = qf_symmetry(Q, b2)*A1
+  if A1 == identity_matrix(QQ, n):
+    return [b1, b2]
+  bs = (A1 - 1).columns()
+  for b3 in bs:
+    if b3 != 0:
+      break
+  A1 = qf_symmetry(Q, b3)*A1
+  if A1 == identity_matrix(QQ, n):
+    return [b1, b2, b3]
+  bs = (A1 - 1).columns()
+  for b4 in bs:
+    if b4 != 0:
+      break
+  A1 = qf_symmetry(Q, b4)*A1
+  if A1 == identity_matrix(QQ, n):
+    return [b1, b2, b3, b4]
+  bs = (A1 - 1).columns()
+  for b5 in bs:
+    if b5 != 0:
+      break
+  return [b1, b2, b3, b4, b5]
+
+  
+
 
 def qf_automorphism_symmetries_proper(Q, A):
 
@@ -604,6 +762,42 @@ def qf_spin_norm(Q, A):
 
   return prod([Q.change_ring(QQ)(v) for v in qf_automorphism_symmetries_proper(Q, A)]).squarefree_part()
 
+def qf_theta_pol(Q, P, N):
+
+  L = Q.vectors_by_length(N)
+  T = [0]
+  auts = Q.number_of_automorphisms()
+  for i in range(1, len(L)):
+    Lt = L[i]
+    T.append(sum([QQ(P(x1 = v[0])(x2 = v[1])(x3 = v[2])(x4 = v[3])(x5 = v[4])) for v in Lt])*i^2)
+  return T
+
+def qf_theta_pols(Qs, Ps, N):
+
+  h = sum([len(Ps[i]) for i in range(len(Ps))])
+  vt = [qf_theta_pol(Qs[i], Ps[i][j], N) for i in range(len(Qs)) for j in range(len(Ps[i]))]
+  M = matrix(QQ, h, [vt[i][j] for i in range(h) for j in range(N)])
+  return M
+
+def qf_spin_norm_vector(Q, A, v, p):
+
+  A = A.change_ring(GF(p))
+  s = (v\(A*v)).det()
+  if s == 1:
+    return 1
+  elif s == -1: 
+    return -1
+  else:
+    print('error')
+
+def qf_radical(Q, p):
+
+  return Q.matrix().change_ring(GF(p)).kernel().matrix().transpose()
+  #return vector([ZZ(gp.matkermod(Q, N)[i + 1, 1]) for i in range(5)])
+
+
+ 
+
 
 #######################################
 #######################################
@@ -624,17 +818,46 @@ def nu(d, n):
 class quinary_module():
   
   """
-  Class of quinary qfs module used to compute
-  quinary orthogonal modular forms.
+  Class of quinary qfs module used to compute quinary orthogonal modular forms.
+  The class can be initialized with a quadratic form, or with a nonsquare level N.
+  If it is initialized with a level N, it uses a quinary quadratic form of level N, with maximal Hasse conductor.
   """
 
   def __init__(self, *args):
 
     self._thetaN = 10
-    if isinstance(args[0], sage.quadratic_forms.quadratic_form.QuadraticForm):
+    if parent(args[0]) == ZZ:
+      N = args[0]
+      if is_square(N):
+        print("N cannot be square")
+      else:
+        gens = [QuadraticForm(g.representative()) for g in QuadraticForm.genera((5,0), 2*N, even=True) if g.level() == 2*N]
+        aux = -1
+        for Q_aux in gens:
+          xi_con = prod([p for p in prime_divisors(N) if Q_aux.xi_rec(p) == -1])
+          if xi_con > aux:
+            aux = xi_con
+            Q = Q_aux
+        self._Q = Q
+        self._disc = self._Q.disc()
+        self._tpinit = False
+        self._vectors = {}
+        for p in self._disc.prime_divisors():
+          if p == 2:
+            self._vectors[2] = qf_radical(self._Q, 4)
+          else:
+            self._vectors[p] = qf_radical(self._Q, p)
+      
+    elif isinstance(args[0], sage.quadratic_forms.quadratic_form.QuadraticForm):
       self._Q = args[0]
       self._disc = self._Q.disc()
       self._tpinit = False
+      self._vectors = {}
+      for p in self._Q.disc().prime_divisors():
+        if p == 2:
+          self._vectors[2] = qf_radical(self._Q, 4)
+        else:
+          self._vectors[p] = qf_radical(self._Q, p)
     else:
       Qs = args[0]
       self._Q = Qs[0]
@@ -656,11 +879,19 @@ class quinary_module():
     self._tpsd = {}
     self._tp2s = {}
     self._tp2sd = {}
+    self._tpns = {}
+    self._tpn2s = {}
+    self._tpdns = {}
+    self._tpdn2s = {}
+    self._tpj2s = {}
+    self._tpdj2s = {}
     self._eigen_tp = {}
     self._eigen_tp2 = {}
     self._charpols = {}
     self._ambiguous = {}
     self._not_ambiguous = {}
+    self._invariants = {}
+    self._invariantsj2 = {}
     self.PB = 6
     self.Q0 = self._Q
     #self.spins = [x for x in self._disc.divisors() if is_squarefree(x)]
@@ -685,6 +916,14 @@ class quinary_module():
     M = matrix(h, [vt[i][j] for i in range(h) for j in range(h)])
     return M.kernel()
 
+  def theta_kernel_pol(self, n):
+
+    Tp = self.Tp_d_n(2, 1, n)
+    Qs = qmod._iso_dict.values()
+    Ps = [[sum([basis_n(n)[i]*qmod._invariants[n][1][k][i, j] for i in range(len(basis_n(n)))]) for j in range(qmod._invariants[n][1][k].ncols())] for k in range(len(Qs))]
+    h = sum([len(x) for x in Ps])
+    return qf_theta_pols(Qs, Ps, h+1).kernel()
+
   def theta_kernel_dimension(self):
 
     """
@@ -693,7 +932,19 @@ class quinary_module():
 
     return self.theta_kernel().dimension()
 
-  def _Tp_init(self, p):
+  def _spin_norm(self, A):
+
+    s = 1
+    for p in self._vectors:
+      if p == 2:
+        dp = qf_spin_norm_vector(self._Q, A, self._vectors[p], 4)
+      else:
+        dp = qf_spin_norm_vector(self._Q, A, self._vectors[p], p)
+      if dp == -1:
+        s*= p
+    return s
+
+  def _Tp_init(self, p, spin_method = 1):
 
     """
     Initialize the module, finds the genus using p-neighbors
@@ -746,7 +997,10 @@ class quinary_module():
             X*= X.det()
         else:
             X = identity_matrix(5)
+        if spin_method == 1:
           s = qf_spin_norm(Q, Q_transformations[R]*A*X*Q_transformations[iso_dict[j]]^-1)
+        elif spin_method == 2:
+          s = self._spin_norm(Q_transformations[R]*A*X*Q_transformations[iso_dict[j]]^-1)
         if (i, j) in tp:
           if s in tp[(i,j)]:
             tp[(i,j)][s]+= 1
@@ -763,7 +1017,7 @@ class quinary_module():
     self._transformations = Q_transformations
     return tp
 
-  def Tp(self, p):
+  def Tp(self, p, spin_method = 1):
 
     """
     Returns the data to compute the Hecke operator. 
@@ -771,7 +1025,7 @@ class quinary_module():
     """
 
     if not self._tpinit:
-      return self._Tp_init(p)
+      return self._Tp_init(p, spin_method = spin_method)
 
     if p in self._tps:
       return self._tps[p]
@@ -799,7 +1053,10 @@ class quinary_module():
               X*= X.det()
               j = iso_dict_inv[Qiso]
               break
+        if spin_method == 1:
           s = qf_spin_norm(self.Q0, self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
+        elif spin_method == 2:
+          s = self._spin_norm(self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
         if s in tp[(i,j)]:
           tp[(i,j)][s]+= 1
         else:
@@ -807,7 +1064,98 @@ class quinary_module():
     self._tps[p] = tp
     return tp
 
-  def Tp_d(self, p, d):
+  def Tp_j2(self, p, k = 1, spin_method = 1):
+
+    if not self._tpinit:
+      self._Tp_init(2)
+
+    if (p, k) in self._tpj2s:
+      return self._tpj2s[(p, k)]
+    iso_dict = self._iso_dict
+    iso_dict_theta = self._iso_dict_theta
+    iso_dict_inv = self._iso_dict_inv
+    g = len(iso_dict)
+    tpj2 = {(i, j):{} for i in range(g) for j in range(g)}
+    for i in range(g):
+      Q = iso_dict[i]
+      if k == 1:
+        neighQ = qf_find_p_neighbors(Q, p, 1)
+      else:
+        neighQ = qf_find_p2_neighbors(Q, p, 1)
+      for R, A in neighQ:
+        Rtheta = R.theta_series(self._thetaN)
+        Rpossible_isos = iso_dict_theta[Rtheta]
+        if len(Rpossible_isos) == 1:
+          Qiso = Rpossible_isos[0]
+          X = R.is_globally_equivalent_to(Qiso, 1)
+          X*= X.det()
+          j = iso_dict_inv[Qiso]
+        else:
+          for Qiso in Rpossible_isos:
+            X = R.is_globally_equivalent_to(Qiso, 1)
+            if not X == False:
+              X*= X.det()
+              j = iso_dict_inv[Qiso]
+              break
+        if spin_method == 1:
+          s = qf_spin_norm(self.Q0, self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
+        elif spin_method == 2:
+          s = self._spin_norm(self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
+        if s in tpj2[(i,j)]:
+          tpj2[(i,j)][s]+= actionext((A*X)^-1)
+        else:
+          tpj2[(i,j)][s] = actionext((A*X)^-1)
+    self._tpj2s[(p, k)] = tpj2
+    return tpj2
+
+
+  def Tp_n(self, p, n, k = 1, spin_method = 1):
+
+    if not self._tpinit:
+      self._Tp_init(2)
+
+    if (p, n, k) in self._tpns:
+      return self._tpns[(p, n, k)]
+
+    iso_dict = self._iso_dict
+    iso_dict_theta = self._iso_dict_theta
+    iso_dict_inv = self._iso_dict_inv
+    g = len(iso_dict)
+    tpn = {(i, j):{} for i in range(g) for j in range(g)}
+    for i in range(g):
+      Q = iso_dict[i]
+      if k == 1:
+        neighQ = qf_find_p_neighbors(Q, p, 1)
+      else:
+        neighQ = qf_find_p2_neighbors(Q, p, 1)
+      for R, A in neighQ:
+        Rtheta = R.theta_series(self._thetaN)
+        Rpossible_isos = iso_dict_theta[Rtheta]
+        if len(Rpossible_isos) == 1:
+          Qiso = Rpossible_isos[0]
+          X = R.is_globally_equivalent_to(Qiso, 1)
+          X*= X.det()
+          j = iso_dict_inv[Qiso]
+        else:
+          for Qiso in Rpossible_isos:
+            X = R.is_globally_equivalent_to(Qiso, 1)
+            if not X == False:
+              X*= X.det()
+              j = iso_dict_inv[Qiso]
+              break
+        if spin_method == 1:
+          s = qf_spin_norm(self.Q0, self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
+        elif spin_method == 2:
+          s = self._spin_norm(self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
+        if s in tpn[(i,j)]:
+          tpn[(i,j)][s]+= action_n((A*X)^-1, n)
+        else:
+          tpn[(i,j)][s] = action_n((A*X)^-1, n)
+    self._tps[(p, n, k)] = tpn
+    return tpn
+
+
+  def Tp_d(self, p, d, spin_method = 1):
 
       """
       Returns the Hecke operator Tp1 at p with character nu_d, with nu_1 the
@@ -817,7 +1165,7 @@ class quinary_module():
       if p in self._tpsd:
           if d in self._tpsd[p]:
               return self._tpsd[p][d]
-      tp = self.Tp(p)
+      tp = self.Tp(p, spin_method = spin_method)
       ncols = len(self._iso_dict)
       tpd = matrix(ncols)
       for (i, j) in tp:
@@ -832,8 +1180,13 @@ class quinary_module():
           not_amb = list(range(ncols))
           for i in self._iso_dict:
               R = self._iso_dict[i]
+              B = self._transformations[R]
               for A in R.automorphisms():
                   if A.det() == 1:
+                      #s = qf_spin_norm(R, A)
+                      if spin_method == 2:
+                        s = self._spin_norm(B*A*B^-1)
+                      elif spin_method == 1:
                         s = qf_spin_norm(R, A)
                       if nu(d, s) == -1:
                           amb.append(i)
@@ -850,8 +1203,73 @@ class quinary_module():
           self._tpsd[p] = {d:tpd}
       return tpd
 
+  #TODO: falta ver que pasa cuando d|p
 
-  def Tp2(self, p):
+  def Tp_d_n(self, p, d, n, k = 1):
+    
+    """
+    Asumo que ya se hizo la inicialización.
+    """
+
+    if (p, d, n, k) in self._tpdns:
+      return self._tpdns[(p, d, n, k)]
+    tpn = self.Tp_n(p, n, k)
+    if n not in self._invariants:
+      self._invariants[n] = {}
+    if d not in self._invariants[n]:
+      self._invariants[n][d] = {}
+      for i in self._iso_dict:
+        q = self._iso_dict[i]
+        inv = qf_invariant3(q, d, n)
+        self._invariants[n][d][i] = inv.matrix().transpose()
+    ncols = sum([m.ncols() for m in self._invariants[n][d].values()])
+    tpdn = matrix(QQ, ncols)
+    i0 = 0
+    for i in range(len(self._iso_dict)):
+      j0 = 0
+      for j in range(len(self._iso_dict)):
+        Mij = self._invariants[n][d][i]\(sum([tpn[(i, j)][kk]*nu(d, kk) for kk in tpn[(i, j)]])*self._invariants[n][d][j]*p^n)
+        tpdn[i0:i0+Mij.nrows(), j0:j0+Mij.ncols()] = Mij
+        j0+= Mij.ncols()
+      i0+= Mij.nrows()
+    if self.disc()%p == 0 and (self.disc()//p)%p != 0:
+      tpdn*= nu(d, p)
+      tpdn+= nu(d, p)*p^n
+    self._tpdns[(p, d, n, k)] = tpdn
+    return tpdn
+    
+  def Tp_d_j2(self, p, d, k = 1):
+
+    """
+    Asumo que ya se hizo la inicialización.
+    """
+
+    if (p, d, k) in self._tpdj2s:
+      return self._tpdj2s[(p, d, k)]
+    tpj2 = self.Tp_j2(p, k)
+    if d not in self._invariantsj2:
+      self._invariantsj2[d] = {}
+      for i in self._iso_dict:
+        q = self._iso_dict[i]
+        inv = qf_invariantj2(q, d)
+        self._invariantsj2[d][i] = inv.matrix().transpose()
+    ncols = sum([m.ncols() for m in self._invariantsj2[d].values()])
+    tpdj2 = matrix(QQ, ncols)
+    i0 = 0
+    for i in range(len(self._iso_dict)):
+      j0 = 0
+      for j in range(len(self._iso_dict)):
+        Mij = self._invariantsj2[d][i]\(sum([tpj2[(i, j)][kk]*nu(d, kk) for kk in tpj2[(i, j)]])*self._invariantsj2[d][j]*p^k)
+        tpdj2[i0:i0+Mij.nrows(), j0:j0+Mij.ncols()] = Mij
+        j0+= Mij.ncols()
+      i0+= Mij.nrows()
+    if self.disc()%p == 0 and (self.disc()//p)%p != 0:
+      tpdj2*= nu(d, p)
+    self._tpdj2s[(p, d, k)] = tpdj2
+    return tpdj2
+
+
+  def Tp2(self, p, spin_method = 1):
 
     """
     Returns the data to compute the Hecke operator.
@@ -888,7 +1306,10 @@ class quinary_module():
               X*= X.det()
               j = iso_dict_inv[Qiso]
               break
+        if spin_method == 1:
           s = qf_spin_norm(self.Q0, self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
+        elif spin_method == 2:
+          s = self._spin_norm(self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
         if s in tp2[(i,j)]:
           tp2[(i,j)][s]+= 1
         else:
@@ -896,7 +1317,7 @@ class quinary_module():
     self._tp2s[p] = tp2
     return tp2
 
-  def Tp2_d(self, p, d):
+  def Tp2_d(self, p, d, spin_method = 1):
 
       """
       Returns the Hecke operator Tp2 at p with character nu_d, with nu_1 the
@@ -907,7 +1328,7 @@ class quinary_module():
       if p in self._tp2sd:
           if d in self._tp2sd[p]:
               return self._tp2sd[p][d]
-      tp2 = self.Tp2(p)
+      tp2 = self.Tp2(p, spin_method = spin_method)
       ncols = len(self._iso_dict)
       tp2d = matrix(ncols)
       for (i, j) in tp2:
@@ -938,7 +1359,7 @@ class quinary_module():
           self._tp2sd[p] = {d:tp2d}
       return tp2d
 
-  def eigenvalue_tp_d(self, p, v, d):
+  def eigenvalue_tp_d(self, p, v, d, spin_method = 1):
 
     """
     Given a vector v that is a Hecke-eigenvector for
@@ -972,7 +1393,10 @@ class quinary_module():
       if j in not_amb:
         X = R.is_globally_equivalent_to(Qiso, 1)
         X*= X.det()
+        if spin_method == 1:
           s = qf_spin_norm(self.Q0, self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
+        elif spin_method == 2:
+          s = s = self._spin_norm(self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
         w[change(j)]+= nu(d, s)
     v_dual = vector([v[change(j)]*ZZ(self._iso_dict[j].number_of_automorphisms()) for j in range(len(self._iso_dict)) if j in not_amb])
     eigen = v_dual*w/v_dual[i]*nu(d, p)
@@ -1036,7 +1460,7 @@ class quinary_module():
     return eigen
 
 
-  def eigenvalue_tp2_d(self, p, v, d):
+  def eigenvalue_tp2_d(self, p, v, d, spin_method = 1):
 
     """
     Given a vector v that is a Hecke-eigenvector for
@@ -1070,7 +1494,10 @@ class quinary_module():
       if j in not_amb:
         X = R.is_globally_equivalent_to(Qiso, 1)
         X*= X.det()
+        if spin_method == 1:
           s = qf_spin_norm(self.Q0, self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
+        elif spin_method == 2:
+          s = self._spin_norm(self._transformations[Q]*A*X*self._transformations[Qiso]^-1)
         sk = nu(d, s)
         w[change(j)]+= sk
     v_dual = vector([v[change(j)]*ZZ(self._iso_dict[j].number_of_automorphisms()) for j in range(len(self._iso_dict)) if j in not_amb])
@@ -1116,17 +1543,17 @@ class quinary_module():
     self._ambiguous = ambs
     self._not_ambiguous = not_ambs
 
-  def decompose(self, d = 1, bound = 40):
+  def decompose(self, d = 1, k = 1, bound = 40):
 
     """
     Finds the decomposition in Hecke-eigenspaces.
     """
     p = 2
     T = self.Tp_d(p, d)
-    todo = [ ZZ**len(self._not_ambiguous[d]) ]
+    todo = [ ZZ**T.nrows() ]
     decomp = []
     while todo:
-        T = self.Tp_d(p, d)
+        T = self.Tp_d(p, d).change_ring(ZZ)
         todo_next = []
         for V in todo:
             for A, is_irred in T.decomposition_of_subspace(V):
@@ -1140,7 +1567,7 @@ class quinary_module():
             break
     return decomp + todo
 
-  def eigentraces(self, space, bound, d = 1):
+  def eigentraces(self, space, bound, d = 1, k = 1):
 
     """
     Returns the traces for Tp1 of the Hecke-eigenspaces.
@@ -1155,3 +1582,15 @@ class quinary_module():
     """
 
     return [ self.Tp2_d(p, d).restrict(space).trace() for p in primes(bound) if self.disc()%p != 0]
+
+
+#q = quins[61][0]
+#qmod = quinary_module(q)
+#T2 = qmod.Tp_d(2, 1).transpose()
+#T3 = qmod.Tp_d(3, 1).transpose()
+#VM = T2.decomposition()[-1][0].matrix().transpose()
+#v = vector((1, 0, 0, 0, 0, 0, 0, 0))
+#MM = Matrix([T2^i*v*VM for i in range(6)]).transpose()
+#f = T2.fcp()[-1][0]
+#L.<a> = NumberField(f)
+#MM.inverse()*((T3*v)*VM)*vector((1, a, a^2, a^3, a^4, a^5))