text: """
  Local Density Congruence
  """
  
  
  ##########################################################################
  ## Methods which compute the local densities for representing a number
  ## by a quadratic form at a prime (possibly subject to additional
  ## congruence conditions).
  ##########################################################################
  
  from copy import deepcopy
  
  from sage.sets.set import Set
  from sage.rings.rational_field import QQ
  from sage.rings.arith import valuation
  from sage.misc.misc import verbose
  
  from sage.quadratic_forms.count_local_2 import count_modp__by_gauss_sum
  
  
  
  
  def count_modp_solutions__by_Gauss_sum(self, p, m):
      """
      Returns the number of solutions of `Q(x) = m (mod p)` of a
      non-degenerate quadratic form over the finite field `Z/pZ`,
      where `p` is a prime number > 2.
  
      Note: We adopt the useful convention that a zero-dimensional
      quadratic form has exactly one solution always (i.e. the empty
      vector).
  
      These are defined in Table 1 on p363 of Hanke's "Local
      Densities..." paper.
  
      INPUT:
  
          `p` -- a prime number > 2
  
          `m` -- an integer
  
      OUTPUT:
  
          an integer >= 0
  
      EXAMPLES::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1])
          sage: [Q.count_modp_solutions__by_Gauss_sum(3, m)  for m in range(3)]
          [9, 6, 12]
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,1,2])
          sage: [Q.count_modp_solutions__by_Gauss_sum(3, m)  for m in range(3)]
          [9, 12, 6]
  
      """
      if self.dim() == 0:
          return 1
      else:
          return count_modp__by_gauss_sum(self.dim(), p, m, self.Gram_det())
  
  
  
  
  
  def local_good_density_congruence_odd(self, p, m, Zvec, NZvec):
      """
      Finds the Good-type local density of Q representing `m` at `p`.
      (Assuming that `p` > 2 and Q is given in local diagonal form.)
  
      The additional congruence condition arguments Zvec and NZvec can
      be either a list of indices or None.  Zvec = [] is equivalent to
      Zvec = None which both impose no additional conditions, but NZvec
      = [] returns no solutions always while NZvec = None imposes no
      additional condition.
  
      TO DO: Add type checking for Zvec, NZvec, and that Q is in local
      normal form.
  
      INPUT:
  
          Q -- quadratic form assumed to be diagonal and p-integral
  
          `p` -- a prime number
  
          `m` -- an integer
  
          Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
  
      OUTPUT:
  
          a rational number
  
      EXAMPLES::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
          sage: Q.local_good_density_congruence_odd(3, 1, None, None)
          2/3
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
          sage: Q.local_good_density_congruence_odd(3, 1, None, None)
          8/9
  
      """
      n = self.dim()
  
      ## Put the Zvec congruence condition in a standard form
      if Zvec is None:
          Zvec = []
  
  
      ## Sanity Check on Zvec and NZvec:
      ## -------------------------------
      Sn = Set(range(n))
      if (Zvec is not None) and (len(Set(Zvec) + Sn) > n):
          raise RuntimeError("Zvec must be a subset of {0, ..., n-1}.")
      if (NZvec is not None) and (len(Set(NZvec) + Sn) > n):
          raise RuntimeError("NZvec must be a subset of {0, ..., n-1}.")
  
  
  
      ## Assuming Q is diagonal, find the indices of the p-unit (diagonal) entries
      UnitVec = [i  for i in range(n)  if (self[i,i] % p) != 0]
      NonUnitVec = list(Set(range(n)) - Set(UnitVec))
  
  
      ## Take cases on the existence of additional non-zero congruence conditions (mod p)
      UnitVec_minus_Zvec = list(Set(UnitVec) - Set(Zvec))
      NonUnitVec_minus_Zvec = list(Set(NonUnitVec) - Set(Zvec))
      Q_Unit_minus_Zvec = self.extract_variables(UnitVec_minus_Zvec)
  
      if (NZvec is None):
          if m % p != 0:
              total = Q_Unit_minus_Zvec.count_modp_solutions__by_Gauss_sum(p, m) * p**len(NonUnitVec_minus_Zvec)          ## m != 0 (mod p)
          else:
              total = (Q_Unit_minus_Zvec.count_modp_solutions__by_Gauss_sum(p, m) - 1) * p**len(NonUnitVec_minus_Zvec)    ## m == 0 (mod p)
  
      else:
          UnitVec_minus_ZNZvec = list(Set(UnitVec) - (Set(Zvec) + Set(NZvec)))
          NonUnitVec_minus_ZNZvec = list(Set(NonUnitVec) - (Set(Zvec) + Set(NZvec)))
          Q_Unit_minus_ZNZvec = self.extract_variables(UnitVec_minus_ZNZvec)
  
          if m % p != 0:         ## m != 0 (mod p)
              total = Q_Unit_minus_Zvec.count_modp_solutions__by_Gauss_sum(p, m) * p**len(NonUnitVec_minus_Zvec) \
                      - Q_Unit_minus_ZNZvec.count_modp_solutions__by_Gauss_sum(p, m) * p**len(NonUnitVec_minus_ZNZvec)
          else:                  ## m == 0 (mod p)
              total = (Q_Unit_minus_Zvec.count_modp_solutions__by_Gauss_sum(p, m) - 1) * p**len(NonUnitVec_minus_Zvec) \
                      - (Q_Unit_minus_ZNZvec.count_modp_solutions__by_Gauss_sum(p, m) - 1) * p**len(NonUnitVec_minus_ZNZvec)
  
      ## Return the Good-type representation density
      good_density = QQ(total) / p**(n-1)
      return good_density
  
  
  
  
  def local_good_density_congruence_even(self, m, Zvec, NZvec):
      """
      Finds the Good-type local density of Q representing `m` at `p=2`.
      (Assuming Q is given in local diagonal form.)
  
  
      The additional congruence condition arguments Zvec and NZvec can
      be either a list of indices or None.  Zvec = [] is equivalent to
      Zvec = None which both impose no additional conditions, but NZvec
      = [] returns no solutions always while NZvec = None imposes no
      additional condition.
  
      WARNING: Here the indices passed in Zvec and NZvec represent
      indices of the solution vector `x` of Q(`x`) = `m (mod p^k)`, and *not*
      the Jordan components of Q.  They therefore are required (and
      assumed) to include either all or none of the indices of a given
      Jordan component of Q.  This is only important when `p=2` since
      otherwise all Jordan blocks are 1x1, and so there the indices and
      Jordan blocks coincide.
  
      TO DO: Add type checking for Zvec, NZvec, and that Q is in local
      normal form.
  
  
      INPUT:
  
          Q -- quadratic form assumed to be block diagonal and 2-integral
  
          `p` -- a prime number
  
          `m` -- an integer
  
          Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
  
      OUTPUT:
  
          a rational number
  
      EXAMPLES::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
          sage: Q.local_good_density_congruence_even(1, None, None)
          1
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
          sage: Q.local_good_density_congruence_even(1, None, None)
          1
          sage: Q.local_good_density_congruence_even(2, None, None)
          3/2
          sage: Q.local_good_density_congruence_even(3, None, None)
          1
          sage: Q.local_good_density_congruence_even(4, None, None)
          1/2
  
      ::
  
          sage: Q = QuadraticForm(ZZ, 4, range(10))
          sage: Q[0,0] = 5
          sage: Q[1,1] = 10
          sage: Q[2,2] = 15
          sage: Q[3,3] = 20
          sage: Q
          Quadratic form in 4 variables over Integer Ring with coefficients:
          [ 5 1 2 3 ]
          [ * 10 5 6 ]
          [ * * 15 8 ]
          [ * * * 20 ]
          sage: Q.theta_series(20)
          1 + 2*q^5 + 2*q^10 + 2*q^14 + 2*q^15 + 2*q^16 + 2*q^18 + O(q^20)
          sage: Q.local_normal_form(2)
          Quadratic form in 4 variables over Integer Ring with coefficients:
          [ 0 1 0 0 ]
          [ * 0 0 0 ]
          [ * * 0 1 ]
          [ * * * 0 ]
          sage: Q.local_good_density_congruence_even(1, None, None)
          3/4
          sage: Q.local_good_density_congruence_even(2, None, None)
          1
          sage: Q.local_good_density_congruence_even(5, None, None)
          3/4
  
      """
      n = self.dim()
  
      ## Put the Zvec congruence condition in a standard form
      if Zvec is None:
          Zvec = []
  
  
      ## Sanity Check on Zvec and NZvec:
      ## -------------------------------
      Sn = Set(range(n))
      if (Zvec is not None) and (len(Set(Zvec) + Sn) > n):
          raise RuntimeError("Zvec must be a subset of {0, ..., n-1}.")
      if (NZvec is not None) and (len(Set(NZvec) + Sn) > n):
          raise RuntimeError("NZvec must be a subset of {0, ..., n-1}.")
  
  
  
      ## Find the indices of x for which the associated Jordan blocks are non-zero mod 8    TO DO: Move this to special Jordan block code separately!
      ## -------------------------------------------------------------------------------
      Not8vec = []
      for i in range(n):
  
          ## DIAGNOSTIC
          verbose(" i = " + str(i))
          verbose(" n = " + str(n))
          verbose(" Not8vec = " + str(Not8vec))
  
          nz_flag = False
  
          ## Check if the diagonal entry isn't divisible 8
          if  ((self[i,i] % 8) != 0):
              nz_flag = True
  
          ## Check appropriate off-diagonal entries aren't divisible by 8
          else:
  
              ## Special check for first off-diagonal entry
              if ((i == 0) and ((self[i,i+1] % 8) != 0)):
                  nz_flag = True
  
              ## Special check for last off-diagonal entry
              elif ((i == n-1) and ((self[i-1,i] % 8) != 0)):
                  nz_flag = True
  
              ## Check for the middle off-diagonal entries
              else:
                  if ( (i > 0)  and  (i < n-1)  and  (((self[i,i+1] % 8) != 0) or ((self[i-1,i] % 8) != 0)) ):
                      nz_flag = True
  
          ## Remember the (vector) index if it's not part of a Jordan block of norm divisible by 8
          if (nz_flag == True):
              Not8vec += [i]
  
  
  
      ## Compute the number of Good-type solutions mod 8:
      ## ------------------------------------------------
  
      ## Setup the indexing sets for additional zero congruence solutions
      Q_Not8 = self.extract_variables(Not8vec)
      Not8 = Set(Not8vec)
      Is8 = Set(range(n)) - Not8
  
      Z = Set(Zvec)
      Z_Not8 = Not8.intersection(Z)
      Z_Is8 = Is8.intersection(Z)
      Is8_minus_Z = Is8 - Z_Is8
  
  
      ## DIAGNOSTIC
      verbose("Z = " + str(Z))
      verbose("Z_Not8 = " + str(Z_Not8))
      verbose("Z_Is8 = " + str(Z_Is8))
      verbose("Is8_minus_Z = " + str(Is8_minus_Z))
  
  
      ## Take cases on the existence of additional non-zero congruence conditions (mod 2)
      if NZvec is None:
          total = (4 ** len(Z_Is8)) * (8 ** len(Is8_minus_Z)) \
              * Q_Not8.count_congruence_solutions__good_type(2, 3, m, list(Z_Not8), None)
      else:
          ZNZ = Z + Set(NZvec)
          ZNZ_Not8 = Not8.intersection(ZNZ)
          ZNZ_Is8 = Is8.intersection(ZNZ)
          Is8_minus_ZNZ = Is8 - ZNZ_Is8
  
          ## DIAGNOSTIC
          verbose("ZNZ = " + str(ZNZ))
          verbose("ZNZ_Not8 = " + str(ZNZ_Not8))
          verbose("ZNZ_Is8 = " + str(ZNZ_Is8))
          verbose("Is8_minus_ZNZ = " + str(Is8_minus_ZNZ))
  
          total = (4 ** len(Z_Is8)) * (8 ** len(Is8_minus_Z)) \
              * Q_Not8.count_congruence_solutions__good_type(2, 3, m, list(Z_Not8), None) \
              - (4 ** len(ZNZ_Is8)) * (8 ** len(Is8_minus_ZNZ)) \
              * Q_Not8.count_congruence_solutions__good_type(2, 3, m, list(ZNZ_Not8), None)
  
  
      ## DIAGNOSTIC
      verbose("total = " + str(total))
  
  
      ## Return the associated Good-type representation density
      good_density = QQ(total) / 8**(n-1)
      return good_density
  
  
  
  
  
  
  
  
  def local_good_density_congruence(self, p, m, Zvec=None, NZvec=None):
      """
      Finds the Good-type local density of Q representing `m` at `p`.
      (Front end routine for parity specific routines for p.)
  
      TO DO: Add Documentation about the additional congruence
      conditions Zvec and NZvec.
  
  
  
      INPUT:
  
          Q -- quadratic form assumed to be block diagonal and p-integral
  
          `p` -- a prime number
  
          `m` -- an integer
  
          Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
  
      OUTPUT:
  
          a rational number
  
      EXAMPLES::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
          sage: Q.local_good_density_congruence(2, 1, None, None)
          1
          sage: Q.local_good_density_congruence(3, 1, None, None)
          2/3
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
          sage: Q.local_good_density_congruence(2, 1, None, None)
          1
          sage: Q.local_good_density_congruence(3, 1, None, None)
          8/9
  
      """
      ## DIAGNOSTIC
      verbose(" In local_good_density_congruence with ")
      verbose(" Q is: \n" + str(self))
      verbose(" p = " + str(p))
      verbose(" m = " + str(m))
      verbose(" Zvec = " + str(Zvec))
      verbose(" NZvec = " + str(NZvec))
  
      ## Put the Zvec congruence condition in a standard form
      if Zvec is None:
          Zvec = []
  
  
      n = self.dim()
  
      ## Sanity Check on Zvec and NZvec:
      ## -------------------------------
      Sn = Set(range(n))
      if (Zvec is not None) and (len(Set(Zvec) + Sn) > n):
          raise RuntimeError("Zvec must be a subset of {0, ..., n-1}.")
      if (NZvec is not None) and (len(Set(NZvec) + Sn) > n):
          raise RuntimeError("NZvec must be a subset of {0, ..., n-1}.")
  
  
  
      ## Check that Q is in local normal form -- should replace this with a diagonalization check?
      ##   (it often may not be since the reduction procedure
      ##   often mixes up the order of the valuations...)
      #
      #if (self != self.local_normal_form(p))
      #    print "Warning in local_good_density_congruence: Q is not in local normal form! \n";
  
  
  
  
      ## Decide which routine to use to compute the Good-type density
      if (p > 2):
          return self.local_good_density_congruence_odd(p, m, Zvec, NZvec)
  
      if (p == 2):
          #print "\n Using the (p=2) Local_Good_Density_Even routine! \n"
          return self.local_good_density_congruence_even(m, Zvec, NZvec)
  
      raise RuntimeError("\n Error in Local_Good_Density: The 'prime' p = " + str(p) + " is < 2. \n")
  
  
  
  
  
  
  
  def local_zero_density_congruence(self, p, m, Zvec=None, NZvec=None):
      """
      Finds the Zero-type local density of Q representing `m` at `p`,
      allowing certain congruence conditions mod p.
  
  
      INPUT:
  
          Q -- quadratic form assumed to be block diagonal and `p`-integral
  
          `p` -- a prime number
  
          `m` -- an integer
  
          Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
  
      OUTPUT:
  
          a rational number
  
      EXAMPLES::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
          sage: Q.local_zero_density_congruence(2, 2, None, None)
          0
          sage: Q.local_zero_density_congruence(2, 4, None, None)
          1/2
          sage: Q.local_zero_density_congruence(3, 6, None, None)
          0
          sage: Q.local_zero_density_congruence(3, 9, None, None)
          2/9
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
          sage: Q.local_zero_density_congruence(2, 2, None, None)
          0
          sage: Q.local_zero_density_congruence(2, 4, None, None)
          1/4
          sage: Q.local_zero_density_congruence(3, 6, None, None)
          0
          sage: Q.local_zero_density_congruence(3, 9, None, None)
          8/81
  
      """
      ## DIAGNOSTIC
      verbose(" In local_zero_density_congruence with ")
      verbose(" Q is: \n" + str(self))
      verbose(" p = " + str(p))
      verbose(" m = " + str(m))
      verbose(" Zvec = " + str(Zvec))
      verbose(" NZvec = " + str(NZvec))
  
      ## Put the Zvec congruence condition in a standard form
      if Zvec is None:
          Zvec = []
  
  
      n = self.dim()
  
      ## Sanity Check on Zvec and NZvec:
      ## -------------------------------
      Sn = Set(range(n))
      if (Zvec is not None) and (len(Set(Zvec) + Sn) > n):
          raise RuntimeError("Zvec must be a subset of {0, ..., n-1}.")
      if (NZvec is not None) and (len(Set(NZvec) + Sn) > n):
          raise RuntimeError("NZvec must be a subset of {0, ..., n-1}.")
  
  
      p2 = p * p
  
      ## Check some conditions for no zero-type solutions to exist
      if ((m % (p2) != 0) or (NZvec is not None)):
          return 0
  
      ## Use the reduction procedure to return the result
      return self.local_density_congruence(p, m / p2, None, None) / p**(self.dim() - 2)
  
  
  
  
  
  
  
  def local_badI_density_congruence(self, p, m, Zvec=None, NZvec=None):
      """
      Finds the Bad-type I local density of Q representing `m` at `p`.
      (Assuming that p > 2 and Q is given in local diagonal form.)
  
  
      INPUT:
  
          Q -- quadratic form assumed to be block diagonal and `p`-integral
  
          `p` -- a prime number
  
          `m` -- an integer
  
          Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
  
      OUTPUT:
  
          a rational number
  
      EXAMPLES::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
          sage: Q.local_badI_density_congruence(2, 1, None, None)
          0
          sage: Q.local_badI_density_congruence(2, 2, None, None)
          1
          sage: Q.local_badI_density_congruence(2, 4, None, None)
          0
          sage: Q.local_badI_density_congruence(3, 1, None, None)
          0
          sage: Q.local_badI_density_congruence(3, 6, None, None)
          0
          sage: Q.local_badI_density_congruence(3, 9, None, None)
          0
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
          sage: Q.local_badI_density_congruence(2, 1, None, None)
          0
          sage: Q.local_badI_density_congruence(2, 2, None, None)
          0
          sage: Q.local_badI_density_congruence(2, 4, None, None)
          0
          sage: Q.local_badI_density_congruence(3, 2, None, None)
          0
          sage: Q.local_badI_density_congruence(3, 6, None, None)
          0
          sage: Q.local_badI_density_congruence(3, 9, None, None)
          0
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,3,3,9])
          sage: Q.local_badI_density_congruence(3, 1, None, None)
          0
          sage: Q.local_badI_density_congruence(3, 3, None, None)
          4/3
          sage: Q.local_badI_density_congruence(3, 6, None, None)
          4/3
          sage: Q.local_badI_density_congruence(3, 9, None, None)
          0
          sage: Q.local_badI_density_congruence(3, 18, None, None)
          0
  
  
      """
      ## DIAGNOSTIC
      verbose(" In local_badI_density_congruence with ")
      verbose(" Q is: \n" + str(self))
      verbose(" p = " + str(p))
      verbose(" m = " + str(m))
      verbose(" Zvec = " + str(Zvec))
      verbose(" NZvec = " + str(NZvec))
  
      ## Put the Zvec congruence condition in a standard form
      if Zvec is None:
          Zvec = []
  
  
      n = self.dim()
  
  
  
      ## Sanity Check on Zvec and NZvec:
      ## -------------------------------
      Sn = Set(range(n))
      if (Zvec is not None) and (len(Set(Zvec) + Sn) > n):
          raise RuntimeError("Zvec must be a subset of {0, ..., n-1}.")
      if (NZvec is not None) and (len(Set(NZvec) + Sn) > n):
          raise RuntimeError("NZvec must be a subset of {0, ..., n-1}.")
  
  
  
      ## Define the indexing set S_0, and determine if S_1 is empty:
      ## -----------------------------------------------------------
      S0 = []
      S1_empty_flag = True    ## This is used to check if we should be computing BI solutions at all!
                              ## (We should really to this earlier, but S1 must be non-zero to proceed.)
  
      ## Find the valuation of each variable (which will be the same over 2x2 blocks),
      ## remembering those of valuation 0 and if an entry of valuation 1 exists.
      for i in range(n):
  
          ## Compute the valuation of each index, allowing for off-diagonal terms
          if (self[i,i] == 0):
              if (i == 0):
                  val = valuation(self[i,i+1], p)    ## Look at the term to the right
              else:
                  if (i == n-1):
                      val = valuation(self[i-1,i], p)    ## Look at the term above
                  else:
                      val = valuation(self[i,i+1] + self[i-1,i], p)    ## Finds the valuation of the off-diagonal term since only one isn't zero
          else:
              val = valuation(self[i,i], p)
  
          if (val == 0):
              S0 += [i]
          elif (val == 1):
              S1_empty_flag = False    ## Need to have a non-empty S1 set to proceed with Bad-type I reduction...
  
  
  
  
  
      ## Check that S1 is non-empty and p|m to proceed, otherwise return no solutions.
      if (S1_empty_flag == True) or (m % p != 0):
          return 0
  
      ## Check some conditions for no bad-type I solutions to exist
      if (NZvec is not None) and (len(Set(S0).intersection(Set(NZvec))) != 0):
          return 0
  
  
  
      ## Check that the form is primitive...                     WHY DO WE NEED TO DO THIS?!?
      if (S0 == []):
          print " Using Q = " + str(self)
          print " and p = " + str(p)
          raise RuntimeError("Oops! The form is not primitive!")
  
  
  
      ## DIAGNOSTIC
      verbose(" m = " + str(m) + "   p = " + str(p))
      verbose(" S0 = " + str(S0))
      verbose(" len(S0) = " + str(len(S0)))
  
  
  
      ## Make the form Qnew for the reduction procedure:
      ## -----------------------------------------------
      Qnew = deepcopy(self)                                          ## TO DO:  DO THIS WITHOUT A copy(). =)
      for i in range(n):
          if i in S0:
              Qnew[i,i] = p * Qnew[i,i]
              if ((p == 2) and (i < n-1)):
                  Qnew[i,i+1] = p * Qnew[i,i+1]
          else:
              Qnew[i,i] = Qnew[i,i] / p
              if ((p == 2) and (i < n-1)):
                  Qnew[i,i+1] = Qnew[i,i+1] / p
  
  
  
      ## DIAGNOSTIC
      verbose("\n\n Check of Bad-type I reduction: \n")
      verbose(" Q is " + str(self))
      verbose(" Qnew is " + str(Qnew))
      verbose(" p = " + str(p))
      verbose(" m / p = " + str(m/p))
      verbose(" NZvec " + str(NZvec))
  
  
  
      ## Do the reduction
      Zvec_geq_1 = list(Set([i  for i in Zvec  if i not in S0]))
      if NZvec is None:
          NZvec_geq_1 = NZvec
      else:
          NZvec_geq_1 = list(Set([i  for i in NZvec  if i not in S0]))
  
      return QQ(p**(1 - len(S0))) * Qnew.local_good_density_congruence(p, m / p, Zvec_geq_1, NZvec_geq_1)
  
  
  
  
  
  
  
  def local_badII_density_congruence(self, p, m, Zvec=None, NZvec=None):
      """
      Finds the Bad-type II local density of Q representing `m` at `p`.
      (Assuming that `p` > 2 and Q is given in local diagonal form.)
  
  
       INPUT:
  
          Q -- quadratic form assumed to be block diagonal and p-integral
  
          `p` -- a prime number
  
          `m` -- an integer
  
          Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
  
      OUTPUT:
  
          a rational number
  
      EXAMPLES::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
          sage: Q.local_badII_density_congruence(2, 1, None, None)
          0
          sage: Q.local_badII_density_congruence(2, 2, None, None)
          0
          sage: Q.local_badII_density_congruence(2, 4, None, None)
          0
          sage: Q.local_badII_density_congruence(3, 1, None, None)
          0
          sage: Q.local_badII_density_congruence(3, 6, None, None)
          0
          sage: Q.local_badII_density_congruence(3, 9, None, None)
          0
          sage: Q.local_badII_density_congruence(3, 27, None, None)
          0
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,3,3,9,9])
          sage: Q.local_badII_density_congruence(3, 1, None, None)
          0
          sage: Q.local_badII_density_congruence(3, 3, None, None)
          0
          sage: Q.local_badII_density_congruence(3, 6, None, None)
          0
          sage: Q.local_badII_density_congruence(3, 9, None, None)
          4/27
          sage: Q.local_badII_density_congruence(3, 18, None, None)
          4/9
  
      """
      ## DIAGNOSTIC
      verbose(" In local_badII_density_congruence with ")
      verbose(" Q is: \n" + str(self))
      verbose(" p = " + str(p))
      verbose(" m = " + str(m))
      verbose(" Zvec = " + str(Zvec))
      verbose(" NZvec = " + str(NZvec))
  
      ## Put the Zvec congruence condition in a standard form
      if Zvec is None:
          Zvec = []
  
  
      n = self.dim()
  
  
      ## Sanity Check on Zvec and NZvec:
      ## -------------------------------
      Sn = Set(range(n))
      if (Zvec is not None) and (len(Set(Zvec) + Sn) > n):
          raise RuntimeError("Zvec must be a subset of {0, ..., n-1}.")
      if (NZvec is not None) and (len(Set(NZvec) + Sn) > n):
          raise RuntimeError("NZvec must be a subset of {0, ..., n-1}.")
  
  
      ## Define the indexing sets S_i:
      ## -----------------------------
      S0 = []
      S1 = []
      S2plus = []
  
      for i in range(n):
  
          ## Compute the valuation of each index, allowing for off-diagonal terms
          if (self[i,i] == 0):
              if (i == 0):
                  val = valuation(self[i,i+1], p)    ## Look at the term to the right
              elif (i == n-1):
                  val = valuation(self[i-1,i], p)    ## Look at the term above
              else:
                  val = valuation(self[i,i+1] + self[i-1,i], p)    ## Finds the valuation of the off-diagonal term since only one isn't zero
          else:
              val = valuation(self[i,i], p)
  
          ## Sort the indices into disjoint sets by their valuation
          if (val == 0):
              S0 += [i]
          elif (val == 1):
              S1 += [i]
          elif (val >= 2):
              S2plus += [i]
  
  
  
      ## Check that S2 is non-empty and p^2 divides m to proceed, otherwise return no solutions.
      p2 = p * p
      if (S2plus == []) or (m % p2 != 0):
          return 0
  
  
  
  
      ## Check some conditions for no bad-type II solutions to exist
      if (NZvec is not None) and (len(Set(S2plus).intersection(Set(NZvec))) == 0):
          return 0
  
  
  
      ## Check that the form is primitive...                     WHY IS THIS NECESSARY?
      if (S0 == []):
          print " Using Q = " + str(self)
          print " and p = " + str(p)
          raise RuntimeError("Oops! The form is not primitive!")
  
  
  
  
      ## DIAGNOSTIC
      verbose("\n Entering BII routine ")
      verbose(" S0 is " + str(S0))
      verbose(" S1 is " + str(S1))
      verbose(" S2plus is " + str(S2plus))
      verbose(" m = " + str(m) + "   p = " + str(p))
  
  
  
  
  
      ## Make the form Qnew for the reduction procedure:
      ## -----------------------------------------------
      Qnew = deepcopy(self)                                          ## TO DO:  DO THIS WITHOUT A copy(). =)
      for i in range(n):
          if i in S2plus:
              Qnew[i,i] = Qnew[i,i] / p2
              if (p == 2) and (i < n-1):
                  Qnew[i,i+1] = Qnew[i,i+1] / p2
  
      ## DIAGNOSTIC
      verbose("\n\n Check of Bad-type II reduction: \n")
      verbose(" Q is " + str(self))
      verbose(" Qnew is " + str(Qnew))
  
  
  
      ## Perform the reduction formula
      Zvec_geq_2 = list(Set([i  for i in Zvec  if i in S2plus]))
      if NZvec is None:
          NZvec_geq_2 = NZvec
      else:
          NZvec_geq_2 = list(Set([i  for i in NZvec  if i in S2plus]))
  
      return QQ(p**(len(S2plus) + 2 - n)) \
          * (Qnew.local_density_congruence(p, m / p2, Zvec_geq_2, NZvec_geq_2) \
          - Qnew.local_density_congruence(p, m / p2, S2plus , NZvec_geq_2))
  
  
  
  
  
  
  def local_bad_density_congruence(self, p, m, Zvec=None, NZvec=None):
      """
      Finds the Bad-type local density of Q representing
      `m` at `p`, allowing certain congruence conditions mod `p`.
  
      INPUT:
  
          Q -- quadratic form assumed to be block diagonal and p-integral
  
          `p` -- a prime number
  
          `m` -- an integer
  
          Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
  
      OUTPUT:
  
          a rational number
  
      EXAMPLES::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
          sage: Q.local_bad_density_congruence(2, 1, None, None)
          0
          sage: Q.local_bad_density_congruence(2, 2, None, None)
          1
          sage: Q.local_bad_density_congruence(2, 4, None, None)
          0
          sage: Q.local_bad_density_congruence(3, 1, None, None)
          0
          sage: Q.local_bad_density_congruence(3, 6, None, None)
          0
          sage: Q.local_bad_density_congruence(3, 9, None, None)
          0
          sage: Q.local_bad_density_congruence(3, 27, None, None)
          0
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,3,3,9,9])
          sage: Q.local_bad_density_congruence(3, 1, None, None)
          0
          sage: Q.local_bad_density_congruence(3, 3, None, None)
          4/3
          sage: Q.local_bad_density_congruence(3, 6, None, None)
          4/3
          sage: Q.local_bad_density_congruence(3, 9, None, None)
          4/27
          sage: Q.local_bad_density_congruence(3, 18, None, None)
          4/9
          sage: Q.local_bad_density_congruence(3, 27, None, None)
          8/27
  
      """
      return self.local_badI_density_congruence(p, m, Zvec, NZvec) + self.local_badII_density_congruence(p, m, Zvec, NZvec)
  
  
  
  #########################################################
  ## local_density and local_density_congruence routines ##
  #########################################################
  
  def local_density_congruence(self, p, m, Zvec=None, NZvec=None):
      """
      Finds the local density of Q representing `m` at `p`,
      allowing certain congruence conditions mod `p`.
  
      INPUT:
  
          Q -- quadratic form assumed to be block diagonal and p-integral
  
          `p` -- a prime number
  
          `m` -- an integer
  
          Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
  
      OUTPUT:
  
          a rational number
  
      EXAMPLES::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
          sage: Q.local_density_congruence(p=2, m=1, Zvec=None, NZvec=None)
          1
          sage: Q.local_density_congruence(p=3, m=1, Zvec=None, NZvec=None)
          8/9
          sage: Q.local_density_congruence(p=5, m=1, Zvec=None, NZvec=None)
          24/25
          sage: Q.local_density_congruence(p=7, m=1, Zvec=None, NZvec=None)
          48/49
          sage: Q.local_density_congruence(p=11, m=1, Zvec=None, NZvec=None)
          120/121
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
          sage: Q.local_density_congruence(2, 1, None, None)
          1
          sage: Q.local_density_congruence(2, 2, None, None)
          1
          sage: Q.local_density_congruence(2, 4, None, None)
          3/2
          sage: Q.local_density_congruence(3, 1, None, None)
          2/3
          sage: Q.local_density_congruence(3, 6, None, None)
          4/3
          sage: Q.local_density_congruence(3, 9, None, None)
          14/9
          sage: Q.local_density_congruence(3, 27, None, None)
          2
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,3,3,9,9])
          sage: Q.local_density_congruence(3, 1, None, None)
          2
          sage: Q.local_density_congruence(3, 3, None, None)
          4/3
          sage: Q.local_density_congruence(3, 6, None, None)
          4/3
          sage: Q.local_density_congruence(3, 9, None, None)
          2/9
          sage: Q.local_density_congruence(3, 18, None, None)
          4/9
  
      """
      return self.local_good_density_congruence(p, m, Zvec, NZvec) \
                  + self.local_zero_density_congruence(p, m, Zvec, NZvec) \
                  + self.local_bad_density_congruence(p, m, Zvec, NZvec)
  
  
  
  def local_primitive_density_congruence(self, p, m, Zvec=None, NZvec=None):
      """
      Finds the primitive local density of Q representing
      `m` at `p`, allowing certain congruence conditions mod `p`.
  
      Note: The following routine is not used internally, but is included for consistency.
  
      INPUT:
  
          Q -- quadratic form assumed to be block diagonal and p-integral
  
          `p` -- a prime number
  
          `m` -- an integer
  
          Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
  
      OUTPUT:
  
          a rational number
  
      EXAMPLES::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
          sage: Q.local_primitive_density_congruence(p=2, m=1, Zvec=None, NZvec=None)
          1
          sage: Q.local_primitive_density_congruence(p=3, m=1, Zvec=None, NZvec=None)
          8/9
          sage: Q.local_primitive_density_congruence(p=5, m=1, Zvec=None, NZvec=None)
          24/25
          sage: Q.local_primitive_density_congruence(p=7, m=1, Zvec=None, NZvec=None)
          48/49
          sage: Q.local_primitive_density_congruence(p=11, m=1, Zvec=None, NZvec=None)
          120/121
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
          sage: Q.local_primitive_density_congruence(2, 1, None, None)
          1
          sage: Q.local_primitive_density_congruence(2, 2, None, None)
          1
          sage: Q.local_primitive_density_congruence(2, 4, None, None)
          1
          sage: Q.local_primitive_density_congruence(3, 1, None, None)
          2/3
          sage: Q.local_primitive_density_congruence(3, 6, None, None)
          4/3
          sage: Q.local_primitive_density_congruence(3, 9, None, None)
          4/3
          sage: Q.local_primitive_density_congruence(3, 27, None, None)
          4/3
  
      ::
  
          sage: Q = DiagonalQuadraticForm(ZZ, [1,3,3,9,9])
          sage: Q.local_primitive_density_congruence(3, 1, None, None)
          2
          sage: Q.local_primitive_density_congruence(3, 3, None, None)
          4/3
          sage: Q.local_primitive_density_congruence(3, 6, None, None)
          4/3
          sage: Q.local_primitive_density_congruence(3, 9, None, None)
          4/27
          sage: Q.local_primitive_density_congruence(3, 18, None, None)
          4/9
          sage: Q.local_primitive_density_congruence(3, 27, None, None)
          8/27
          sage: Q.local_primitive_density_congruence(3, 81, None, None)
          8/27
          sage: Q.local_primitive_density_congruence(3, 243, None, None)
          8/27
  
      """
      return self.local_good_density_congruence(p, m, Zvec, NZvec) \
                  + self.local_bad_density_congruence(p, m, Zvec, NZvec)
  
  
  