def mass_at_two_by_counting_mod_power(self, k):
    """
    Computes the local mass at `p=2` assuming that it's stable `(mod 2^k)`.

    Note: This is **way** too slow to be useful, even when k=1!!!

    TO DO: Remove this routine, or try to compile it!


    INPUT:
        k -- an integer >= 1

    OUTPUT:
        a rational number

    EXAMPLE::

        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1])
        sage: Q.mass_at_two_by_counting_mod_power(1)
        4

    """
    R = IntegerModRing(2**k)
    Q1 = self.base_change_to(R)
    n = self.dim()
    MS = MatrixSpace(R, n)

    ct = sum([1  for x in mrange([2**k] * (n**2))  if Q1(MS(x)) == Q1])   ## Count the solutions mod 2^k
    two_mass = ZZ(1)/2 * (ZZ(ct) / ZZ(2)**(k*n*(n-1)/2))
    return two_mass

 def __iter__(self):
     """
     Return an iterator over the elements of this group.
     
     EXAMPLES::
     
         sage: G = AbelianGroup([2,3], names = "ab")
         sage: [a for a in G]
         [1, b, b^2, a, a*b, a*b^2]
         sage: L = list(G); L
         [1, b, b^2, a, a*b, a*b^2]
     
     The returned list is a reference; mutating it does not allow the
     user to (accidentally?) alter the computed generators::
     
         sage: L[0] = 0
         sage: list(G)
         [1, b, b^2, a, a*b, a*b^2]
         sage: G = AbelianGroup([1], names="a")
         sage: list(G)
         [1]
         sage: G = AbelianGroup([])
         sage: G.list()
         [1]
         sage: list(G)
         [1]
     """
     invs = self.invariants()
     for t in mrange(invs):
         yield AbelianGroupElement(self, t)

 def dual_vectors( self):
     """
     Return a set of representatives
     for $L^#/L$.
     """
     D,U,V = self.gram_matrix().smith_form()
     # hence D = U * self.gram_matrix() * V
     if None == self.__dual_vectors:
         W = V*D**-1
         S = self.space()
         self.__dual_vectors = [ W*S(v) for v in mrange( D.diagonal())]
     return self.__dual_vectors

 def cs_range( f, subset = None):
     """
     For a symmetric semi-positive integral matrix $f$,
     return a list of all integral $n$-vectors $v$ such that
     $x^tfx - (v*x)^2 >= 0$ for all $x$.
     """
     n = f.dimensions()[0]
     b = vector( s.isqrt() for s in f.diagonal())
     zv = vector([0]*n)
     box = [b - vector(t) for t in mrange( (2*b + vector([1]*n)).list()) if b - vector(t) > zv]
     if subset:
         box = [v for v in box if v in subset]
     rge = [v for v in box if min( (f - matrix(v).transpose()*matrix(v)).eigenvalues()) >=0]
     return rge

    def list(self):
        """
        Return tuple of all elements of this group.

        EXAMPLES::

            sage: G = AbelianGroup([2,3], names="ab")
            sage: Gd = G.dual_group(names="AB")
            sage: Gd.list()
            (1, B, B^2, A, A*B, A*B^2)
        """
        if not(self.is_finite()):
           raise NotImplementedError("Group must be finite")
        invs = self.gens_orders()
        T = mrange(invs)
        n = self.order()
        L = tuple( self(t) for t in T )
        return L

    def __iter__(self):
        r"""
        Return an iterator of all ideal classes in this class group.

        EXAMPLES::

            sage: K.<a> = NumberField(x^4 + 23)
            sage: G = K.class_group()
            sage: G
            Class group of order 3 with structure C3 of Number Field
            in a with defining polynomial x^4 + 23
            sage: list(G)
            [Trivial principal fractional ideal class,
             Fractional ideal class (2, 1/4*a^3 - 1/4*a^2 + 1/4*a - 1/4),
             Fractional ideal class (2, 1/2*a^2 + 1/2)]
            sage: G.list()
            (Trivial principal fractional ideal class,
             Fractional ideal class (2, 1/4*a^3 - 1/4*a^2 + 1/4*a - 1/4),
             Fractional ideal class (2, 1/2*a^2 + 1/2))

        TESTS::

            sage: K.<a> = NumberField(x^2 + 1)
            sage: G = K.class_group()
            sage: G
            Class group of order 1 of Number Field in a with defining polynomial x^2 + 1
            sage: list(G)
            [Trivial principal fractional ideal class]
            sage: G.list()
            (Trivial principal fractional ideal class,)
        """
        from sage.misc.mrange import mrange
        orders = self.gens_orders()
        T = mrange(orders)
        g = self.gens()
        for t in T:
            I = self(1)
            for i, j in enumerate(t):
                I *= g[i]**j
            yield I
        if not T:
            yield self(1)

 def list(self):
     """
     Return list of all elements of this group.
     
     EXAMPLES::
     
         sage: G = AbelianGroup([2,3], names = "ab")
         sage: Gd = DualAbelianGroup(G, names = "AB")
         sage: Gd.list()
         [1, B, B^2, A, A*B, A*B^2]
     """
     try:
         return list(self.__list)
     except AttributeError:
         pass
     if not(self.is_finite()):
        raise NotImplementedError, "Group must be finite"
     invs = self.invariants()
     T = mrange(invs)
     n = self.order()
     L = [DualAbelianGroupElement(self, t) for t in T]
     self.__list = L
     return list(self.__list)

def Watson_mass_at_2(self):
    """
    Returns the local mass of the quadratic form when `p=2`, according
    to Watson's Theorem 1 of "The 2-adic density of a quadratic form"
    in Mathematika 23 (1976), pp 94--106.

    INPUT:
        none

    OUTPUT:
        a rational number

    EXAMPLES::

        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1])
        sage: Q.Watson_mass_at_2()               ## WARNING:  WE NEED TO CHECK THIS CAREFULLY!
        384

    """
    ## Make a 0-dim'l quadratic form (for initialization purposes)
    Null_Form = copy.deepcopy(self)
    Null_Form.__init__(ZZ, 0)

    ## Step 0: Compute Jordan blocks and bounds of the scales to keep track of
    Jordan_Blocks = self.jordan_blocks_by_scale_and_unimodular(2)
    scale_list = [B[0]  for B in Jordan_Blocks]
    s_min = min(scale_list)
    s_max = max(scale_list)

    ## Step 1: Compute dictionaries of the diagonal block and 2x2 block for each scale
    diag_dict = dict((i, Null_Form)  for i in range(s_min-2, s_max + 4))     ## Initialize with the zero form
    dim2_dict = dict((i, Null_Form)  for i in range(s_min, s_max + 4))       ## Initialize with the zero form
    for (s,L) in Jordan_Blocks:
        i = 0
        while (i < L.dim()-1) and (L[i,i+1] == 0):      ## Find where the 2x2 blocks start
            i = i + 1
        if i < (L.dim() - 1):
            diag_dict[s] = L.extract_variables(range(i))                ## Diagonal Form
            dim2_dict[s+1] = L.extract_variables(range(i, L.dim()))     ## Non-diagonal Form
        else:
            diag_dict[s] = L

    #print "diag_dict = ", diag_dict
    #print "dim2_dict = ", dim2_dict
    #print "Jordan_Blocks = ", Jordan_Blocks

    ## Step 2: Compute three dictionaries of invariants (for n_j, m_j, nu_j)
    n_dict = dict((j,0)  for j in range(s_min+1, s_max+2))
    m_dict = dict((j,0)  for j in range(s_min, s_max+4))
    for (s,L) in Jordan_Blocks:
        n_dict[s+1] = L.dim()
        if diag_dict[s].dim() == 0:
            m_dict[s+1] = ZZ(1)/ZZ(2) * L.dim()
        else:
            m_dict[s+1] = floor(ZZ(L.dim() - 1) / ZZ(2))
            #print " ==>", ZZ(L.dim() - 1) / ZZ(2), floor(ZZ(L.dim() - 1) / ZZ(2))

    nu_dict = dict((j,n_dict[j+1] - 2*m_dict[j+1])  for j in range(s_min, s_max+1))
    nu_dict[s_max+1] = 0

    #print "n_dict = ", n_dict
    #print "m_dict = ", m_dict
    #print "nu_dict = ", nu_dict

    ## Step 3: Compute the e_j dictionary
    eps_dict = {}
    for j in range(s_min, s_max+3):
        two_form = (diag_dict[j-2] + diag_dict[j] + dim2_dict[j]).scale_by_factor(2)
        j_form = (two_form + diag_dict[j-1]).base_change_to(IntegerModRing(4))

        if j_form.dim() == 0:
            eps_dict[j] = 1
        else:
            iter_vec = [4] * j_form.dim()
            alpha = sum([True  for x in mrange(iter_vec)  if j_form(x) == 0])
            beta = sum([True  for x in mrange(iter_vec)  if j_form(x) == 2])
            if alpha > beta:
                eps_dict[j] = 1
            elif alpha == beta:
                eps_dict[j] = 0
            else:
                eps_dict[j] = -1

    #print "eps_dict = ", eps_dict


    ## Step 4: Compute the quantities nu, q, P, E for the local mass at 2
    nu = sum([j * n_dict[j] * (ZZ(n_dict[j] + 1) / ZZ(2) + \
              sum([n_dict[r]  for r in range(j+1, s_max+2)]))  for j in range(s_min+1, s_max+2)])
    q = sum([sgn(nu_dict[j-1] * (n_dict[j] + sgn(nu_dict[j])))  for j in range(s_min+1, s_max+2)])
    P = prod([ prod([1 - QQ(4)**(-i)  for i in range(1, m_dict[j]+1)])  for j in range(s_min+1, s_max+2)])
    E = prod([ZZ(1)/ZZ(2) * (1 + eps_dict[j] * QQ(2)**(-m_dict[j]))  for j in range(s_min, s_max+3)])

    #print "\nFinal Summary:"
    #print "nu =", nu
    #print "q = ", q
    #print "P = ", P
    #print "E = ", E


#.........这里部分代码省略.........

def minkowski_reduction_for_4vars__SP(self):
    """
    Find a Minkowski-reduced form equivalent to the given one.
    This means that

        Q(`v_k`) <= Q(`s_1 * v_1 + ... + s_n * v_n`)

    for all `s_i` where GCD(`s_k, ... s_n`) = 1.

    Note: When Q has dim <= 4 we can take all `s_i` in {1, 0, -1}.

    References:
        Schulze-Pillot's paper on "An algorithm for computing genera
            of ternary and quaternary quadratic forms", p138.
        Donaldson's 1979 paper "Minkowski Reduction of Integral
            Matrices", p203.

    EXAMPLES::

        sage: Q = QuadraticForm(ZZ,4,[30,17,11,12,29,25,62,64,25,110])
        sage: Q
        Quadratic form in 4 variables over Integer Ring with coefficients:
        [ 30 17 11 12 ]
        [ * 29 25 62 ]
        [ * * 64 25 ]
        [ * * * 110 ]
        sage: Q.minkowski_reduction_for_4vars__SP()
        (
        Quadratic form in 4 variables over Integer Ring with coefficients:
        [ 29 -17 25 4 ]
        [ * 30 -11 5 ]
        [ * * 64 0 ]
        [ * * * 77 ]                                                       ,
        <BLANKLINE>
        [ 0  1  0  0]
        [ 1  0  0 -1]
        [ 0  0  1  0]
        [ 0  0  0  1]
        )
    """
    R = self.base_ring()
    n = self.dim()
    interior_reduced_flag = False
    Q = deepcopy(self)
    M = matrix(R, n, n)
    for i in range(n):
        M[i,i] = 1

    ## Only allow 4-variable forms
    if n != 4:
        raise TypeError("Oops!  The given quadratic form has " + str(n) +  \
                " != 4 variables. =|")


    ## Step 1: Begin the reduction
    done_flag = False
    while not done_flag:

        ## Loop through possible shorter vectors
        done_flag = True
        #print " j_range = ", range(n-1, -1, -1)
        for j in range(n-1, -1, -1):
            for a_first in mrange([2  for i in range(j)]):
                y = [x-1 for x in a_first] + [1] + [0 for k in range(n-1-j)]
                e_j = [0  for k in range(n)]
                e_j[j] = 1
                #print "j = ", j

                ## Reduce if a shorter vector is found
                #print "y = ", y, "   e_j = ", e_j, "\n"
                if Q(y) < Q(e_j):

                    ## Further n=4 computations
                    B_y_vec = Q.matrix() * vector(ZZ, y)
                        ## SP's B = our self.matrix()/2
                        ## SP's A = coeff matrix of his B
                        ## Here we compute the double of both and compare.
                    B_sum = sum([abs(B_y_vec[i])  for i in range(4)  if i != j])
                    A_sum = sum([abs(Q[i,j])  for i in range(4)  if i != j])
                    B_max = max([abs(B_y_vec[i])  for i in range(4)  if i != j])
                    A_max = max([abs(Q[i,j])  for i in range(4)  if i != j])

                    if (B_sum < A_sum) or ((B_sum == A_sum) and (B_max < A_max)):

                        ## Create the transformation matrix
                        M_new = matrix(R, n, n)
                        for k in range(n):
                            M_new[k,k] = 1
                        for k in range(n):
                            M_new[k,j] = y[k]

                        ## Perform the reduction and restart the loop
                        #print "Q_before = ", Q
                        Q = Q(M_new)
                        M = M * M_new
                        done_flag = False

                        ## DIAGNOSTIC
                        #print "Q(y) = ", Q(y)
                        #print "Q(e_j) = ", Q(e_j)
#.........这里部分代码省略.........

def minkowski_reduction(self):
    """
    Find a Minkowski-reduced form equivalent to the given one.
    This means that

    .. MATH::

            Q(v_k) <= Q(s_1 * v_1 + ... + s_n * v_n)

    for all `s_i` where GCD`(s_k, ... s_n) = 1`.

    Note: When Q has dim <= 4 we can take all `s_i` in {1, 0, -1}.

    References:
        Schulze-Pillot's paper on "An algorithm for computing genera
            of ternary and quaternary quadratic forms", p138.
        Donaldson's 1979 paper "Minkowski Reduction of Integral
            Matrices", p203.

    EXAMPLES::

        sage: Q = QuadraticForm(ZZ,4,[30,17,11,12,29,25,62,64,25,110])
        sage: Q
        Quadratic form in 4 variables over Integer Ring with coefficients:
        [ 30 17 11 12 ]
        [ * 29 25 62 ]
        [ * * 64 25 ]
        [ * * * 110 ]
        sage: Q.minkowski_reduction()
        (
        Quadratic form in 4 variables over Integer Ring with coefficients:
        [ 30 17 11 -5 ]
        [ * 29 25 4 ]
        [ * * 64 0 ]
        [ * * * 77 ]                                                       ,
        <BLANKLINE>
        [ 1  0  0  0]
        [ 0  1  0 -1]
        [ 0  0  1  0]
        [ 0  0  0  1]
        )
    """
    R = self.base_ring()
    n = self.dim()
    interior_reduced_flag = False
    Q = deepcopy(self)
    M = matrix(R, n, n)
    for i in range(n):
        M[i,i] = 1


    ## Begin the reduction
    done_flag = False
    while not done_flag:

        ## Loop through possible shorted vectors until
        done_flag = True
        #print " j_range = ", range(n-1, -1, -1)
        for j in range(n-1, -1, -1):
            for a_first in mrange([2  for i in range(j)]):
                y = [x-1 for x in a_first] + [1] + [0 for k in range(n-1-j)]
                e_j = [0  for k in range(n)]
                e_j[j] = 1
                #print "j = ", j

                ## Reduce if a shorter vector is found
                #print "y = ", y, "   e_j = ", e_j, "\n"
                if Q(y) < Q(e_j):

                    ## Create the transformation matrix
                    M_new = matrix(R, n, n)
                    for k in range(n):
                        M_new[k,k] = 1
                    for k in range(n):
                        M_new[k,j] = y[k]

                    ## Perform the reduction and restart the loop
                    #print "Q_before = ", Q
                    Q = Q(M_new)
                    M = M * M_new
                    done_flag = False

                    ## DIAGNOSTIC
                    #print "Q(y) = ", Q(y)
                    #print "Q(e_j) = ", Q(e_j)
                    #print "M_new = ", M_new
                    #print "Q_after = ", Q
                    #print

                if not done_flag:
                    break

            if not done_flag:
                break

    ## Return the results
    return Q, M

def is_globally_equivalent_to(self, other, return_matrix=False, check_theta_to_precision='sturm', check_local_equivalence=True):
    """
    Determines if the current quadratic form is equivalent to the
    given form over ZZ.  If return_matrix is True, then we also return
    the transformation matrix M so that self(M) == other.

    INPUT:
        a QuadraticForm

    OUTPUT:
        boolean, and optionally a matrix

    EXAMPLES::
    
        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
        sage: M = Matrix(ZZ, 4, 4, [1,2,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1])
        sage: Q1 = Q(M)
        sage: Q.(Q1)                                                    # optional -- souvigner
        True
        sage: MM = Q.is_globally_equivalent_to(Q1, return_matrix=True)  # optional -- souvigner
        sage: Q(MM) == Q1                                               # optional -- souvigner
        True
        
    ::
    
        sage: Q1 = QuadraticForm(ZZ, 3, [1, 0, -1, 2, -1, 5])
        sage: Q2 = QuadraticForm(ZZ, 3, [2, 1, 2, 2, 1, 3])
        sage: Q3 = QuadraticForm(ZZ, 3, [8, 6, 5, 3, 4, 2])
        sage: Q1.is_globally_equivalent_to(Q2)                          # optional -- souvigner
        False
        sage: Q1.is_globally_equivalent_to(Q3)                          # optional -- souvigner
        True
        sage: M = Q1.is_globally_equivalent_to(Q3, True) ; M            # optional -- souvigner
        [-1 -1  0]
        [ 1  1  1]
        [-1  0  0]
        sage: Q1(M) == Q3                                               # optional -- souvigner
        True
        
    ::
        
        sage: Q = DiagonalQuadraticForm(ZZ, [1, -1])
        sage: Q.is_globally_equivalent_to(Q)
        Traceback (most recent call last):
        ...
        ValueError: not a definite form in QuadraticForm.is_globally_equivalent_to()
    
    """
    ## only for definite forms
    if not self.is_definite():
        raise ValueError, "not a definite form in QuadraticForm.is_globally_equivalent_to()"

    ## Check that other is a QuadraticForm
    #if not isinstance(other, QuadraticForm):
    if not is_QuadraticForm(other): 
        raise TypeError, "Oops!  You must compare two quadratic forms, but the argument is not a quadratic form. =("


    ## Now use the Souvigner code by default! =)
    return other.is_globally_equivalent__souvigner(self, return_matrix)    ## Note: We switch this because the Souvigner code has the opposite mapping convention to us.  (It takes the second argument to the first!)



    ## ----------------------------------  Unused Code below  ---------------------------------------------------------

    ## Check if the forms are locally equivalent
    if (check_local_equivalence == True):
        if not self.is_locally_equivalent_to(other):
            return False

    ## Check that the forms have the same theta function up to the desired precision (this can be set so that it determines the cusp form)
    if check_theta_to_precision != None:
        if self.theta_series(check_theta_to_precision, var_str='', safe_flag=False) != other.theta_series(check_theta_to_precision, var_str='', safe_flag=False):
            return False


    ## Make all possible matrices which give an isomorphism -- can we do this more intelligently?
    ## ------------------------------------------------------------------------------------------

    ## Find a basis of short vectors for one form, and try to match them with vectors of that length in the other one.
    basis_for_self, self_lengths = self.basis_of_short_vectors(show_lengths=True)
    max_len = max(self_lengths)
    short_vectors_of_other = other.short_vector_list_up_to_length(max_len + 1)

    ## Make the matrix A:e_i |--> v_i to our new basis.
    A = Matrix(basis_for_self).transpose()
    Q2 = A.transpose() * self.matrix() * A       ## This is the matrix of 'self' in the new basis
    Q3 = other.matrix()
        
    ## Determine all automorphisms
    n = self.dim()
    Auto_list = []
    
    ## DIAGNOSTIC
    #print "n = " + str(n)
    #print "pivot_lengths = " + str(pivot_lengths)
    #print "vector_list_by_length = " + str(vector_list_by_length)
    #print "length of vector_list_by_length = " + str(len(vector_list_by_length))
    
    for index_vec in mrange([len(short_vectors_of_other[self_lengths[i]])  for i in range(n)]):
#.........这里部分代码省略.........