def _mul_(self, other):
        """
        Multiply this Hecke operator by another element of the same algebra. If
        the other element is of the form `T_m` for some m, we check whether the
        product is equal to `T_{mn}` and return that; if the product is not
        (easily seen to be) of the form `T_{mn}`, then we calculate the product
        of the two matrices and return a Hecke algebra element defined by that.

        EXAMPLES: We create the space of modular symbols of level
        `11` and weight `2`, then compute `T_2`
        and `T_3` on it, along with their composition.

        ::

            sage: M = ModularSymbols(11)
            sage: t2 = M.hecke_operator(2); t3 = M.hecke_operator(3)
            sage: t2*t3 # indirect doctest
            Hecke operator T_6 on Modular Symbols space of dimension 3 for Gamma_0(11) of weight 2 with sign 0 over Rational Field
            sage: t3.matrix() * t2.matrix()
            [12  0 -2]
            [ 0  2  0]
            [ 0  0  2]
            sage: (t2*t3).matrix()
            [12  0 -2]
            [ 0  2  0]
            [ 0  0  2]

        When we compute `T_2^5` the result is not (easily seen to
        be) a Hecke operator of the form `T_n`, so it is returned
        as a Hecke module homomorphism defined as a matrix::

            sage: t2**5
            Hecke operator on Modular Symbols space of dimension 3 for Gamma_0(11) of weight 2 with sign 0 over Rational Field defined by:
            [243   0 -55]
            [  0 -32   0]
            [  0   0 -32]
        """
        if isinstance(other, HeckeOperator) and other.parent() == self.parent():
            n = None
            if arith.gcd(self.__n, other.__n) == 1:
                n = self.__n * other.__n
            else:
                P = set(arith.prime_divisors(self.domain().level()))
                if P.issubset(set(arith.prime_divisors(self.__n))) and \
                   P.issubset(set(arith.prime_divisors(other.__n))):
                    n = self.__n * other.__n
            if n:
                return HeckeOperator(self.parent(), n)
        # otherwise
        return self.matrix_form() * other

def CS_genus_symbol_list(self, force_recomputation=False):
    """
    Returns the list of Conway-Sloane genus symbols in increasing order of primes dividing 2*det.

    EXAMPLES::

        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3,4])
        sage: Q.CS_genus_symbol_list()
        [Genus symbol at 2 : [[1, 2, 3, 1, 4], [2, 1, 1, 1, 1], [3, 1, 1, 1, 1]],
        Genus symbol at 3 : [[0, 3, 1], [1, 1, -1]]]

    """
    ## Try to use the cached list
    if force_recomputation == False:
        try:
            return self.__CS_genus_symbol_list
        except AttributeError:
            pass

    ## Otherwise recompute and cache the list
    list_of_CS_genus_symbols = [ ]

    for p in prime_divisors(2 * self.det()):
        list_of_CS_genus_symbols.append(self.local_genus_symbol(p))

    self.__CS_genus_symbol_list = list_of_CS_genus_symbols
    return list_of_CS_genus_symbols

    def __unit_gens_primecase(self, p):
        """
        Assuming the modulus is prime, returns the smallest generator
        of the group of units.

        EXAMPLES::

            sage: Zmod(17)._IntegerModRing_generic__unit_gens_primecase(17)
            3
        """
        if p==2:
            return integer_mod.Mod(1,p)
        P = prime_divisors(p-1)
        ord = integer.Integer(p-1)
        one = integer_mod.Mod(1,p)
        x = 2
        while x < p:
            generator = True
            z = integer_mod.Mod(x,p)
            for q in P:
                if z**(ord//q) == one:
                    generator = False
                    break
            if generator:
                return z
            x += 1
        #end for
        assert False, "didn't find primitive root for p=%s"%p

def is_locally_equivalent_to(self, other, check_primes_only=False, force_jordan_equivalence_test=False):
    """
    Determines if the current quadratic form (defined over ZZ) is
    locally equivalent to the given form over the real numbers and the
    `p`-adic integers for every prime p.

    This works by comparing the local Jordan decompositions at every
    prime, and the dimension and signature at the real place.

    INPUT:

    a QuadraticForm

    OUTPUT:

    boolean

    EXAMPLES::

        sage: Q1 = QuadraticForm(ZZ, 3, [1, 0, -1, 2, -1, 5])
        sage: Q2 = QuadraticForm(ZZ, 3, [2, 1, 2, 2, 1, 3])
        sage: Q1.is_globally_equivalent_to(Q2)
        False
        sage: Q1.is_locally_equivalent_to(Q2)
        True

    """
    ## TO IMPLEMENT:
    if self.det() == 0:
        raise NotImplementedError("OOps!  We need to think about whether this still works for degenerate forms...  especially check the signature.")

    ## Check that both forms have the same dimension and base ring
    if (self.dim() != other.dim()) or (self.base_ring() != other.base_ring()):
        return False

    ## Check that the determinant and level agree
    if (self.det() != other.det()) or (self.level() != other.level()):
        return False

    ## -----------------------------------------------------

    ## Test equivalence over the real numbers
    if self.signature() != other.signature():
        return False

    ## Test equivalence over Z_p for all primes
    if (self.base_ring() == ZZ) and (force_jordan_equivalence_test == False):

        ## Test equivalence with Conway-Sloane genus symbols (default over ZZ)
        if self.CS_genus_symbol_list() != other.CS_genus_symbol_list():
            return False
    else:
        ## Test equivalence via the O'Meara criterion.
        for p in prime_divisors(ZZ(2) * self.det()):
            #print "checking the prime p = ", p
            if not self.has_equivalent_Jordan_decomposition_at_prime(other, p):
                return False

    ## All tests have passed!
    return True

def anisotropic_primes(self):
    """
    Returns a list with all of the anisotropic primes of the quadratic form.


    INPUT:

        None

    OUTPUT:

        Returns a list of prime numbers >0.

    EXAMPLES::

        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1])
        sage: Q.anisotropic_primes()
        [2]

    ::

        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
        sage: Q.anisotropic_primes()
        [2]

    ::

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

    """

    ## Look at all prime divisors of 2 * Det(Q) to find the anisotropic primes...
    possible_primes = prime_divisors(2 * self.det())
    AnisoPrimes = []

    ## DIAGNSOTIC
    #print " Possible anisotropic primes are: " + str(possible_primes)

    for p in possible_primes:
        if (self.is_anisotropic(p)):
            AnisoPrimes += [p]

    ## DIAGNSOTIC
    #print " leaving anisotropic_primes..."

    return AnisoPrimes

 def __unit_gens_primecase(self, p):
     if p==2:
         return integer_mod.Mod(1,p)
     P = prime_divisors(p-1)
     ord = integer.Integer(p-1)
     one = integer_mod.Mod(1,p)
     x = 2
     while x < p:
         generator = True
         z = integer_mod.Mod(x,p)
         for q in P:
             if z**(ord//q) == one:
                 generator = False
                 break
         if generator:
             return z
         x += 1
     #end for
     assert False, "didn't find primitive root for p=%s"%p

def mass__by_Siegel_densities(self, odd_algorithm="Pall", even_algorithm="Watson"):
    """
    Gives the mass of transformations (det 1 and -1).

    WARNING: THIS IS BROKEN RIGHT NOW... =(

    Optional Arguments:

    - When p > 2  --  odd_algorithm = "Pall" (only one choice for now)
    - When p = 2  --  even_algorithm = "Kitaoka" or "Watson"

    REFERENCES:

    - Nipp's Book "Tables of Quaternary Quadratic Forms".
    - Papers of Pall (only for p>2) and Watson (for `p=2` -- tricky!).
    - Siegel, Milnor-Hussemoller, Conway-Sloane Paper IV, Kitoaka (all of which
      have problems...)

    EXAMPLES::

        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
        sage: Q.mass__by_Siegel_densities()
        1/384
        sage: Q.mass__by_Siegel_densities() - (2^Q.dim() * factorial(Q.dim()))^(-1)
        0

    ::

        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1])
        sage: Q.mass__by_Siegel_densities()
        1/48
        sage: Q.mass__by_Siegel_densities() - (2^Q.dim() * factorial(Q.dim()))^(-1)
        0

    """
    ## Setup
    n = self.dim()
    s = (n-1) // 2
    if n % 2 != 0:
        char_d = squarefree_part(2*self.det())   ## Accounts for the det as a QF
    else:
        char_d = squarefree_part(self.det())

    ## Form the generic zeta product
    generic_prod = ZZ(2) * (pi)**(-ZZ(n) * (n+1) / 4)
    ##########################################
    generic_prod *= (self.det())**(ZZ(n+1)/2)  ## ***** This uses the Hessian Determinant ********
    ##########################################
    #print "gp1 = ", generic_prod
    generic_prod *= prod([gamma__exact(ZZ(j)/2)  for j in range(1,n+1)])
    #print "\n---", [(ZZ(j)/2, gamma__exact(ZZ(j)/2))  for j in range(1,n+1)]
    #print "\n---", prod([gamma__exact(ZZ(j)/2)  for j in range(1,n+1)])
    #print "gp2 = ", generic_prod
    generic_prod *= prod([zeta__exact(ZZ(j))  for j in range(2, 2*s+1, 2)])
    #print "\n---", [zeta__exact(ZZ(j))  for j in range(2, 2*s+1, 2)]
    #print "\n---", prod([zeta__exact(ZZ(j))  for j in range(2, 2*s+1, 2)])
    #print "gp3 = ", generic_prod
    if (n % 2 == 0):
        generic_prod *= ZZ(1) * quadratic_L_function__exact(n/2, (-1)**(n/2) * char_d)
        #print " NEW = ", ZZ(1) * quadratic_L_function__exact(n/2, (-1)**(n/2) * char_d)
        #print
    #print "gp4 = ", generic_prod

    #print "generic_prod =", generic_prod

    ## Determine the adjustment factors
    adj_prod = 1
    for p in prime_divisors(2 * self.det()):
        ## Cancel out the generic factors
        p_adjustment = prod([1 - ZZ(p)**(-j)  for j in range(2, 2*s+1, 2)])
        if (n % 2 == 0):
            p_adjustment *= ZZ(1) * (1 - kronecker((-1)**(n/2) * char_d, p) * ZZ(p)**(-n/2))
            #print " EXTRA = ", ZZ(1) * (1 - kronecker((-1)**(n/2) * char_d, p) * ZZ(p)**(-n/2))
        #print "Factor to cancel the generic one:", p_adjustment

        ## Insert the new mass factors
        if p == 2:
            if even_algorithm == "Kitaoka":
                p_adjustment = p_adjustment / self.Kitaoka_mass_at_2()
            elif even_algorithm == "Watson":
                p_adjustment = p_adjustment / self.Watson_mass_at_2()
            else:
                raise TypeError("There is a problem -- your even_algorithm argument is invalid.  Try again. =(")
        else:
            if odd_algorithm == "Pall":
                p_adjustment = p_adjustment / self.Pall_mass_density_at_odd_prime(p)
            else:
                raise TypeError("There is a problem -- your optional arguments are invalid.  Try again. =(")

        #print "p_adjustment for p =", p, "is", p_adjustment

        ## Put them together (cumulatively)
        adj_prod *= p_adjustment

    #print "Cumulative adj_prod =", adj_prod

        ## Extra adjustment for the case of a 2-dimensional form.
    #if (n == 2):
    #    generic_prod *= 2

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

#.........这里部分代码省略.........
        for D in BSD.curve.heegner_discriminants_list(10):
            max_height = max(13,BSD.curve.quadratic_twist(D).CPS_height_bound())
            heegner_primes = -1
            while heegner_primes == -1:
                if max_height > 21: break
                heegner_primes, _, exact = BSD.curve.heegner_index_bound(D, max_height=max_height)
                max_height += 1
            if isinstance(heegner_primes, list):
                break
        if not isinstance(heegner_primes, list):
            raise RuntimeError("Tried 10 Heegner discriminants, and heegner_index_bound failed each time.")
        if exact is not False:
            heegner_index = exact
            BSD.heegner_indexes[D] = exact
        else:
            BSD.heegner_index_upper_bound[D] = max(heegner_primes+[1])
        if 2 in heegner_primes:
            heegner_primes.remove(2)
    else: # rank 1
        for D in BSD.curve.heegner_discriminants_list(10):
            I = BSD.curve.heegner_index(D)
            J = I.is_int()
            if J[0] and J[1]>0:
                I = J[1]
            else:
                J = (2*I).is_int()
                if J[0] and J[1]>0:
                    I = J[1]
                else:
                    continue
            heegner_index = I
            BSD.heegner_indexes[D] = I
            break
        heegner_primes = [p for p in arith.prime_divisors(heegner_index) if p!=2]
    
    assert BSD.sha_an in ZZ and BSD.sha_an > 0
    if BSD.curve.has_cm():
        if BSD.curve.analytic_rank() == 0:
            if verbosity > 0:
                print '    p >= 5: true by Rubin'
            BSD.primes.append(3)
        else:
            K = rings.QuadraticField(BSD.curve.cm_discriminant(), 'a')
            D_K = K.disc()
            D_E = BSD.curve.discriminant()
            if len(K.factor(3)) == 1: # 3 does not split in K
                BSD.primes.append(3)
            for p in arith.prime_divisors(D_K):
                if p >= 5:
                    BSD.primes.append(p)
            for p in arith.prime_divisors(D_E):
                if p >= 5 and D_K%p and len(K.factor(p)) == 1: # p is inert in K
                    BSD.primes.append(p)
            for p in heegner_primes:
                if p >= 5 and D_E%p != 0 and D_K%p != 0 and len(K.factor(p)) == 1: # p is good for E and inert in K
                    kolyvagin_primes.append(p)
            for p in arith.prime_divisors(BSD.sha_an):
                if p >= 5 and D_K%p != 0 and len(K.factor(p)) == 1:
                    if BSD.curve.is_good(p):
                        if verbosity > 2 and p in heegner_primes and heegner_index is None:
                            print 'ALERT: Prime p (%d) >= 5 dividing sha_an, good for E, inert in K, in heegner_primes, should not divide the actual Heegner index'
                        # Note that the following check is not entirely
                        # exhaustive, in case there is a p not dividing 
                        # the Heegner index in heegner_primes,
                        # for which only an outer bound was computed
                        if p not in heegner_primes:

#.........这里部分代码省略.........
            1
            sage: m = EllipticCurve('37a1').modular_symbol(use_eclib=True)
            sage: m._scaling
            -1
            sage: m = EllipticCurve('389a1').modular_symbol(use_eclib=True)
            sage: m._scaling
            -1/2
            sage: m = EllipticCurve('389a1').modular_symbol(use_eclib=False)
            sage: m._scaling
            2
            sage: m = EllipticCurve('196a1').modular_symbol(use_eclib=False)
            sage: m._scaling  
            1/2  
            
        Some harder cases fail::
        
            sage: m = EllipticCurve('121b1').modular_symbol(use_eclib=False)
            Warning : Could not normalize the modular symbols, maybe all further results will be multiplied by -1, 2 or -2.
            sage: m._scaling  
            1           
            
        TESTS::

            sage: rk0 = ['11a1', '11a2', '15a1', '27a1', '37b1']
            sage: for la in rk0:  # long time (3s on sage.math, 2011)
            ...          E = EllipticCurve(la)
            ...          me = E.modular_symbol(use_eclib = True)
            ...          ms = E.modular_symbol(use_eclib = False)
            ...          print E.lseries().L_ratio()*E.real_components(), me(0), ms(0)
            1/5 1/5 1/5
            1 1 1
            1/4 1/4 1/4
            1/3 1/3 1/3
            2/3 2/3 2/3

            sage: rk1 = ['37a1','43a1','53a1', '91b1','91b2','91b3']
            sage: [EllipticCurve(la).modular_symbol(use_eclib=True)(0) for la in rk1]  # long time (1s on sage.math, 2011)
            [0, 0, 0, 0, 0, 0]
            sage: for la in rk1:  # long time (8s on sage.math, 2011)
            ...       E = EllipticCurve(la)
            ...       m = E.modular_symbol(use_eclib = True)
            ...       lp = E.padic_lseries(5)
            ...       for D in [5,17,12,8]:
            ...           ED = E.quadratic_twist(D)
            ...           md = sum([kronecker(D,u)*m(ZZ(u)/D) for u in range(D)])
            ...           etaa = lp._quotient_of_periods_to_twist(D)
            ...           assert ED.lseries().L_ratio()*ED.real_components()*etaa == md

        """
        E = self._E
        self._scaling = 1 # by now.
        self._failed_to_scale = False

        if self._sign == 1 :
            at0 = self(0)
            # print 'modular symbol evaluates to ',at0,' at 0'
            if at0 != 0 :
                l1 = self.__lalg__(1)
                if at0 != l1:
                    verbose('scale modular symbols by %s'%(l1/at0))
                    self._scaling = l1/at0
            else :
                # if [0] = 0, we can still hope to scale it correctly by considering twists of E
                Dlist = [5,8,12,13,17,21,24,28,29, 33, 37, 40, 41, 44, 53, 56, 57, 60, 61, 65, 69, 73, 76, 77, 85, 88, 89, 92, 93, 97]  # a list of positive fundamental discriminants
                j = 0
                at0 = 0
                # computes [0]+ for the twist of E by D until one value is non-zero
                while j < 30 and at0 == 0 :
                    D = Dlist[j]
                    # the following line checks if the twist of the newform of E by D is a newform
                    # this is to avoid that we 'twist back'
                    if all( valuation(E.conductor(),ell)<= valuation(D,ell) for ell in prime_divisors(D) ) :
                        at0 = sum([kronecker_symbol(D,u) * self(ZZ(u)/D) for u in range(1,abs(D))])
                    j += 1
                if j == 30 and at0 == 0: # curves like "121b1", "225a1", "225e1", "256a1", "256b1", "289a1", "361a1", "400a1", "400c1", "400h1", "441b1", "441c1", "441d1", "441f1 .. will arrive here
                    self.__scale_by_periods_only__()
                else :
                    l1 = self.__lalg__(D)
                    if at0 != l1:
                        verbose('scale modular symbols by %s found at D=%s '%(l1/at0,D), level=2)
                        self._scaling = l1/at0

        else : # that is when sign = -1
            Dlist = [-3,-4,-7,-8,-11,-15,-19,-20,-23,-24, -31, -35, -39, -40, -43, -47, -51, -52, -55, -56, -59, -67, -68, -71, -79, -83, -84, -87, -88, -91]  # a list of negative fundamental discriminants
            j = 0
            at0 = 0
            while j < 30 and at0 == 0 :
                # computes [0]+ for the twist of E by D until one value is non-zero
                D = Dlist[j]
                if all( valuation(E.conductor(),ell)<= valuation(D,ell) for ell in prime_divisors(D) ) :
                    at0 = - sum([kronecker_symbol(D,u) * self(ZZ(u)/D) for u in range(1,abs(D))])
                j += 1
            if j == 30 and at0 == 0: # no more hope for a normalization
                # we do at least a scaling with the quotient of the periods
                self.__scale_by_periods_only__()
            else :
                l1 = self.__lalg__(D)
                if at0 != l1:
                    verbose('scale modular symbols by %s'%(l1/at0))
                    self._scaling = l1/at0

#.........这里部分代码省略.........
            * QQ(sqrt((2 ** n) *  f) / (u * d)) \
            * abs(QuadraticBernoulliNumber(m, d1) / bernoulli(2*m))



    ## DIAGNOSTIC
    verbose("siegel_product Break 2. \n")


    ## Make the even generic factor
    if ((n % 2) == 0):
        m = n / 2
        d1 = fundamental_discriminant(((-1)**m) * d)  
        f = abs(d1)                                

        ## DIAGNOSTIC
        #cout << " mpz_class(-1)^m = " << (mpz_class(-1)^m) << " and d = " << d << endl;
        #cout << " f = " << f << " and d1 = " << d1 << endl;


        genericfactor = m / QQ(sqrt(f*d)) \
            * ((u/2) ** (m-1)) * (f ** m) \
            / abs(QuadraticBernoulliNumber(m, d1)) \
            * (2 ** m)                                               ## This last factor compensates for using the matrix of 2*Q


    ##return genericfactor
  
  
    ## Omit the generic factors in S and compute them separately
    omit = 1
    include = 1
  
    S_divisors = prime_divisors(S)

    ## DIAGNOSTIC
    #cout << "\n S is " << S << endl;
    #cout << " The Prime divisors of S are :";
    #PrintV(S_divisors);


    for p in S_divisors:    
        Q_normal = self.local_normal_form(p)
    

        ## DIAGNOSTIC
        verbose(" p = " + str(p) + " and its Kronecker symbol (d1/p) = (" + str(d1) + "/" + str(p) + ") is " + str(kronecker_symbol(d1, p)) + "\n")

        omit *= 1 / (1 - (kronecker_symbol(d1, p) / (p**m))) 


        ## DIAGNOSTIC
        verbose(" omit = " + str(omit) + "\n")
        verbose(" Q_normal is \n" + str(Q_normal) + "\n")
        verbose(" Q_normal = \n" + str(Q_normal))
        verbose(" p = " + str(p) + "\n")
        verbose(" u = " +str(u) + "\n")
        verbose(" include = " + str(include) + "\n")


        include *= Q_normal.local_density(p, u)


        ## DIAGNOSTIC
        #cout << " Including the p = " << p << " factor: " << local_density(Q_normal, p, u) << endl; 

    def old_submodule(self, p=None):
        """
        Returns the old or p-old submodule of self, i.e. the sum of the images
        of the degeneracy maps from level `N/p` (for the given prime `p`, or
        for all primes `p` dividing `N` if `p` is not given).

        INPUT:

        - ``p`` - (default: None); if not None, return only the p-old
          submodule.

        OUTPUT: the old or p-old submodule of self

        EXAMPLES::

            sage: m = ModularSymbols(33); m.rank()
            9
            sage: m.old_submodule().rank()
            7
            sage: m.old_submodule(3).rank()
            6
            sage: m.new_submodule(11).rank()
            8

        ::

            sage: e = DirichletGroup(16)([-1, 1])
            sage: M = ModularSymbols(e, 3, sign=1); M
            Modular Symbols space of dimension 4 and level 16, weight 3, character [-1, 1], sign 1, over Rational Field
            sage: M.old_submodule()
            Modular Symbols subspace of dimension 3 of Modular Symbols space of dimension 4 and level 16, weight 3, character [-1, 1], sign 1, over Rational Field

        Illustrate that trac 10664 is fixed::

            sage: ModularSymbols(DirichletGroup(42)[7], 6, sign=1).old_subspace(3)
            Modular Symbols subspace of dimension 0 of Modular Symbols space of dimension 40 and level 42, weight 6, character [-1, -1], sign 1, over Rational Field

        """
        try:
            if self.__is_old[p]:
                return self
        except AttributeError:
            self.__is_old = {}
        except KeyError:
            pass

        if self.rank() == 0:
            self.__is_old[p] = True
            return self
        try:
            return self.__old_submodule[p]
        except AttributeError:
            self.__old_submodule = {}
        except KeyError:
            pass

        # Construct the degeneracy map d.
        N = self.level()
        d = None

        eps = self.character()
        if eps is None:
            f = 1
        else:
            f = eps.conductor()

        if p is None:
            D = arith.prime_divisors(N)
        else:
            if N % p != 0:
                raise ValueError, "p must divide the level."
            D = [p]

        for q in D:
            NN = N//q
            if NN % f == 0:
                M = self.hecke_module_of_level(NN)

                # Here it is vital to pass self as an argument to
                # degeneracy_map, because M and the level N don't uniquely
                # determine self (e.g. the degeneracy map from level 1 to level
                # N could go to Gamma0(N), Gamma1(N) or anything in between)
                d1 = M.degeneracy_map(self, 1).matrix()

                if d is None:
                    d = d1
                else:
                    d = d.stack(d1)
                d = d.stack(M.degeneracy_map(self, q).matrix())
            #end if
        #end for
        if d is None:
            os = self.zero_submodule()
        else:
            os = self.submodule(d.image(), check=False)

        self.__is_old[p] = (os == self)

        os.__is_old = {p:True}
        os._is_full_hecke_module = True
#.........这里部分代码省略.........

    def new_submodule(self, p=None):
        """
        Returns the new or p-new submodule of self.

        INPUT:

        -  ``p`` - (default: None); if not None, return only
           the p-new submodule.

        OUTPUT: the new or p-new submodule of self, i.e. the intersection of
        the kernel of the degeneracy lowering maps to level `N/p` (for the
        given prime `p`, or for all prime divisors of `N` if `p` is not given).

        If self is cuspidal this is a Hecke-invariant complement of the
        corresponding old submodule, but this may break down on Eisenstein
        subspaces (see the amusing example in William Stein's book of a form
        which is new and old at the same time).

        EXAMPLES::

            sage: m = ModularSymbols(33); m.rank()
            9
            sage: m.new_submodule().rank()
            3
            sage: m.new_submodule(3).rank()
            4
            sage: m.new_submodule(11).rank()
            8
        """
        try:
            if self.__is_new[p]:
                return self
        except AttributeError:
            self.__is_new = {}
        except KeyError:
            pass

        if self.rank() == 0:
            self.__is_new[p] = True
            return self
        try:
            return self.__new_submodule[p]
        except AttributeError:
            self.__new_submodule = {}
        except KeyError:
            pass

        # Construct the degeneracy map d.
        N = self.level()
        d = None
        eps = self.character()
        if eps == None:
            f = 1
        else:
            f = eps.conductor()
        if p == None:
            D = arith.prime_divisors(N)
        else:
            if N % p != 0:
                raise ValueError, "p must divide the level."
            D = [p]
        for q in D:
            # Here we are only using degeneracy *lowering* maps, so it is fine
            # to be careless and pass an integer for the level. One needs to be
            # a bit more careful with degeneracy *raising* maps for the Gamma1
            # and GammaH cases.
            if ((N//q) % f) == 0:
                NN = N//q
                d1 = self.degeneracy_map(NN,1).matrix()
                if d is None:
                    d = d1
                else:
                    d = d.augment(d1)
                d = d.augment(self.degeneracy_map(NN,q).matrix())
            #end if
        #end for
        if d is None or d == 0:
            self.__is_new[p] = True
            return self
        else:
            self.__is_new[p] = False
        ns = self.submodule(d.kernel(), check=False)
        ns.__is_new = {p:True}
        ns._is_full_hecke_module = True
        self.__new_submodule[p] = ns
        return ns

    def basis(self, reduce=True):
        r"""
        Produce a basis for the free abelian group of eta-products of level
        N (under multiplication), attempting to find basis vectors of the
        smallest possible degree.

        INPUT:


        -  ``reduce`` - a boolean (default True) indicating
           whether or not to apply LLL-reduction to the calculated basis


        EXAMPLE::

            sage: EtaGroup(5).basis()
            [Eta product of level 5 : (eta_1)^6 (eta_5)^-6]
            sage: EtaGroup(12).basis()
            [Eta product of level 12 : (eta_1)^2 (eta_2)^1 (eta_3)^2 (eta_4)^-1 (eta_6)^-7 (eta_12)^3,
            Eta product of level 12 : (eta_1)^-4 (eta_2)^2 (eta_3)^4 (eta_6)^-2,
            Eta product of level 12 : (eta_1)^-1 (eta_2)^3 (eta_3)^3 (eta_4)^-2 (eta_6)^-9 (eta_12)^6,
            Eta product of level 12 : (eta_1)^1 (eta_2)^-1 (eta_3)^-3 (eta_4)^-2 (eta_6)^7 (eta_12)^-2,
            Eta product of level 12 : (eta_1)^-6 (eta_2)^9 (eta_3)^2 (eta_4)^-3 (eta_6)^-3 (eta_12)^1]
            sage: EtaGroup(12).basis(reduce=False) # much bigger coefficients
            [Eta product of level 12 : (eta_2)^24 (eta_12)^-24,
            Eta product of level 12 : (eta_1)^-336 (eta_2)^576 (eta_3)^696 (eta_4)^-216 (eta_6)^-576 (eta_12)^-144,
            Eta product of level 12 : (eta_1)^-8 (eta_2)^-2 (eta_6)^2 (eta_12)^8,
            Eta product of level 12 : (eta_1)^1 (eta_2)^9 (eta_3)^13 (eta_4)^-4 (eta_6)^-15 (eta_12)^-4,
            Eta product of level 12 : (eta_1)^15 (eta_2)^-24 (eta_3)^-29 (eta_4)^9 (eta_6)^24 (eta_12)^5]

        ALGORITHM: An eta product of level `N` is uniquely
        determined by the integers `r_d` for `d | N` with
        `d < N`, since `\sum_{d | N} r_d = 0`. The valid
        `r_d` are those that satisfy two congruences modulo 24,
        and one congruence modulo 2 for every prime divisor of N. We beef
        up the congruences modulo 2 to congruences modulo 24 by multiplying
        by 12. To calculate the kernel of the ensuing map
        `\ZZ^m \to (\ZZ/24\ZZ)^n`
        we lift it arbitrarily to an integer matrix and calculate its Smith
        normal form. This gives a basis for the lattice.

        This lattice typically contains "large" elements, so by default we
        pass it to the reduce_basis() function which performs
        LLL-reduction to give a more manageable basis.
        """

        N = self.level()
        divs = divisors(N)[:-1]
        s = len(divs)
        primedivs = prime_divisors(N)

        rows = []
        for i in xrange(s):
            # generate a row of relation matrix
            row = [ Mod(divs[i], 24) - Mod(N, 24), Mod(N/divs[i], 24) - Mod(1, 24)]
            for p in primedivs:
                row.append( Mod(12*(N/divs[i]).valuation(p), 24))
            rows.append(row)
        M = matrix(IntegerModRing(24), rows)
        Mlift = M.change_ring(ZZ)
        # now we compute elementary factors of Mlift
        S,U,V = Mlift.smith_form()
        good_vects = []
        for i in xrange(U.nrows()):
            vect = U.row(i)
            nf = (i < S.ncols() and S[i,i]) or 0
            good_vects.append((vect * 24/gcd(nf, 24)).list())
        for v in good_vects:
            v.append(-sum([r for r in v]))
        dicts = []
        for v in good_vects:
            dicts.append({})
            for i in xrange(s):
                dicts[-1][divs[i]] = v[i]
            dicts[-1][N] = v[-1]
        if reduce:
            return self.reduce_basis([ self(d) for d in dicts])
        else:
            return [self(d) for d in dicts]

    def series(self, n=2, quadratic_twist=+1, prec=5):
        r"""
        Returns the `n`-th approximation to the `p`-adic L-series as a
        power series in `T` (corresponding to `\gamma-1` with
        `\gamma=1+p` as a generator of `1+p\ZZ_p`).  Each coefficient
        is a `p`-adic number whose precision is provably correct.
        
        Here the normalization of the `p`-adic L-series is chosen such
        that `L_p(J,1) = (1-1/\alpha)^2 L(J,1)/\Omega_J` where
        `\alpha` is the unit root

        INPUT:
        
            - ``n`` - (default: 2) a positive integer
            - ``quadratic_twist`` - (default: +1) a fundamental
              discriminant of a quadratic field, coprime to the
              conductor of the curve
            - ``prec`` - (default: 5) maximal number of terms of the
              series to compute; to compute as many as possible just
              give a very large number for ``prec``; the result will
              still be correct.

        ALIAS: power_series is identical to series.

        EXAMPLES:

	    sage: J = J0(188)[0]
	    sage: p = 7
	    sage: L = J.padic_lseries(p)
	    sage: L.is_ordinary()
	    True
	    sage: f = L.series(2)
	    sage: f[0]
	    O(7^20)
	    sage: f[1].norm()
	    3 + 4*7 + 3*7^2 + 6*7^3 + 5*7^4 + 5*7^5 + 6*7^6 + 4*7^7 + 5*7^8 + 7^10 + 5*7^11 + 4*7^13 + 4*7^14 + 5*7^15 + 2*7^16 + 5*7^17 + 7^18 + 7^19 + O(7^20)

        """
        n = ZZ(n)
        if n < 1:
            raise ValueError, "n (=%s) must be a positive integer"%n
        if not self.is_ordinary():
            raise ValueError, "p (=%s) must be an ordinary prime"%p
        # check if the conditions on quadratic_twist are satisfied
        D = ZZ(quadratic_twist)
        if D != 1:
            if D % 4 == 0:
                d = D//4
                if not d.is_squarefree() or d % 4 == 1:
                    raise ValueError, "quadratic_twist (=%s) must be a fundamental discriminant of a quadratic field"%D
            else:
                if not D.is_squarefree() or D % 4 != 1:
                    raise ValueError, "quadratic_twist (=%s) must be a fundamental discriminant of a quadratic field"%D
            if gcd(D,self._p) != 1:
                raise ValueError, "quadratic twist (=%s) must be coprime to p (=%s) "%(D,self._p)
            if gcd(D,self._E.conductor())!= 1:
                for ell in prime_divisors(D):
                    if valuation(self._E.conductor(),ell) > valuation(D,ell) :
                        raise ValueError, "can not twist a curve of conductor (=%s) by the quadratic twist (=%s)."%(self._E.conductor(),D)
                    
            
        p = self._p
        if p == 2 and self._normalize :
            print 'Warning : For p=2 the normalization might not be correct !'
        #verbose("computing L-series for p=%s, n=%s, and prec=%s"%(p,n,prec))
        
#        bounds = self._prec_bounds(n,prec)
#        padic_prec = max(bounds[1:]) + 5
        padic_prec = 10
#        verbose("using p-adic precision of %s"%padic_prec)
        
        res_series_prec = min(p**(n-1), prec)
        verbose("using series precision of %s"%res_series_prec)
        
        ans = self._get_series_from_cache(n, res_series_prec,D)
        if not ans is None:
            verbose("found series in cache")
            return ans
 
        K = QQ
        gamma = K(1 + p)
        R = PowerSeriesRing(K,'T',res_series_prec)
        T = R(R.gen(),res_series_prec )
        #L = R(0) 
        one_plus_T_factor = R(1) 
        gamma_power = K(1)
        teich = self.teichmuller(padic_prec)
        p_power = p**(n-1)
#        F = Qp(p,padic_prec)

        verbose("Now iterating over %s summands"%((p-1)*p_power))
        verbose_level = get_verbose()
        count_verb = 0
        alphas = self.alpha()
        #print len(alphas)
        Lprod = []
        self._emb = 0
        if len(alphas) == 2:
            split = True
        else:
#.........这里部分代码省略.........