def q_binomial(n, k, p=None):
    """
    Returns the `q`-binomial coefficient.

    If `p` is unspecified, then it defaults to using the generator `q` for
    a univariate polynomial ring over the integers.

    The q-binomials are computed as a product of cyclotomic polynomials (cf. [CH2006]_).

    REFERENCES:

    .. [CH2006] William Y.C. Chen and Qing-Hu Hou, "Factors of the Gaussian
       coefficients", Discrete Mathematics 306 (2006), 1446-1449.
       http://dx.doi.org/10.1016/j.disc.2006.03.031

    EXAMPLES::

        sage: import sage.combinat.q_analogues as q_analogues
        sage: q_analogues.q_binomial(4,2)
        q^4 + q^3 + 2*q^2 + q + 1
        sage: q_analogues.q_binomial(4,5)
        0
        sage: p = ZZ['p'].0
        sage: q_analogues.q_binomial(4,2,p)
        p^4 + p^3 + 2*p^2 + p + 1
        sage: q_analogues.q_binomial(2,3)
        0

    The `q`-analogue of ``binomial(n,k)`` is currently only defined for
    `n` a nonnegative integer, it is zero for negative k  (trac #11411)::

        sage: q_analogues.q_binomial(5, -1)
        0
    """
    if not (n in ZZ and k in ZZ):
        raise ValueError("Argument (%s, %s) must be integers." % (n, k))
    if n < 0:
        raise NotImplementedError
    if 0 <= k and k <= n:
        if p == None:
            R = ZZ["q"]
        else:
            R = ZZ[p]
        from sage.functions.all import floor

        return prod(
            R.cyclotomic_polynomial(d) for d in range(2, n + 1) if floor(n / d) != floor(k / d) + floor((n - k) / d)
        )
    else:
        return 0

def Pall_mass_density_at_odd_prime(self, p):
    """
    Returns the local representation density of a form (for
    representing itself) defined over `ZZ`, at some prime `p>2`.

    REFERENCES:
        Pall's article "The Weight of a Genus of Positive n-ary Quadratic Forms"
        appearing in Proc. Symp. Pure Math. VIII (1965), pp95--105.

    INPUT:
        `p` -- a prime number > 2.

    OUTPUT:
        a rational number.

    EXAMPLES::

        sage: Q = QuadraticForm(ZZ, 3, [1,0,0,1,0,1])
        sage: Q.Pall_mass_density_at_odd_prime(3)
        [(0, Quadratic form in 3 variables over Integer Ring with coefficients:
        [ 1 0 0 ]
        [ * 1 0 ]
        [ * * 1 ])] [(0, 3, 8)] [8/9] 8/9
        8/9
    """
    ## Check that p is a positive prime -- unnecessary since it's done implicitly in the next step. =)
    if p<=2:
        raise TypeError, "Oops!  We need p to be a prime > 2."

    ## Step 1: Obtain a p-adic (diagonal) local normal form, and
    ## compute the invariants for each Jordan block.
    jordan_list = self.jordan_blocks_by_scale_and_unimodular(p)
    modified_jordan_list = [(a, Q.dim(), Q.det())  for (a,Q) in jordan_list]     ## List of pairs (scale, det)
    #print jordan_list
    #print modified_jordan_list

    ## Step 2: Compute the list of local masses for each Jordan block
    jordan_mass_list = []
    for (s,n,d) in modified_jordan_list:
        generic_factor = prod([1 - p**(-2*j)  for j in range(1, floor((n-1)/2)+1)])
        #print "generic factor: ", generic_factor
        if (n % 2 == 0):
            m = n/2
            generic_factor *= (1 + legendre_symbol(((-1)**m) * d, p) * p**(-m))
        #print "jordan_mass: ", generic_factor
        jordan_mass_list = jordan_mass_list + [generic_factor]

    ## Step 3: Compute the local mass $\al_p$ at p.
        MJL = modified_jordan_list
    s = len(modified_jordan_list)
    M = [sum([MJL[j][1]  for j in range(i, s)])  for i in range(s-1)]    ## Note: It's s-1 since we don't need the last M.
    #print "M = ", M
    nu = sum([M[i] * MJL[i][0] * MJL[i][1]  for i in range(s-1)]) - ZZ(sum([J[0] * J[1] * (J[1]-1)  for J in MJL]))/ZZ(2)
    p_mass = prod(jordan_mass_list)
    p_mass *= 2**(s-1) * p**nu

    print jordan_list, MJL, jordan_mass_list, p_mass

    ## Return the result
    return p_mass

    def test(self,do_twoy=False,up_to_M=0):
        r""" Return the number of digits we believe are correct (at least) 

        EXAMPLES::

        
            sage: G=MySubgroup(Gamma0(1))
            sage: R=mpmath.mpf(9.53369526135355755434423523592877032382125639510725198237579046413534)
            sage: F=MaassWaveformElement(G,R)    
            sage: F.test()
            7



        """
        # If we have a Gamma_0(N) we can use Hecke operators
        if(not do_twoy and (isinstance(self._space,sage.modular.arithgroup.congroup_gamma0.Gamma0_class) or \
            self._space._G.is_congruence)):
            # isinstance(self._space,sage.modular.arithgroup.congroup_sl2z.SL2Z_class_with_category))):
            if(self._space._verbose>1):
                print "Check Hecke relations!"
            er=test_Hecke_relations(2,3,self.coeffs)
            d=floor(-mpmath.log10(er))
            if(self._space._verbose>1):
                print "Hecke is ok up to ",d,"digits!"
            return d
        else:
            # Test two Y's
            nd=self._nd+5
            [M0,Y0]=find_Y_and_M(G,R,nd)
            Y1=Y0*0.95
            C1=Maassform_coeffs(self,R,Mset=M0,Yset=Y0 ,ndigs=nd )[0]
            C2=Maassform_coeffs(self,R,Mset=M0,Yset=Y1 ,ndigs=nd )[0]
            er=mpmath.mpf(0)
            for j in range(2,max(M0/2,up_to_M)):
                t=abs(C1[j]-C2[j])
                print "|C1-C2[",j,"]|=|",C1[j],'-',C2[j],'|=',t
                if(t>er):
                    er=t
            d=floor(-mpmath.log10(er))
            print "Hecke is ok up to ",d,"digits!"
            return d

def katz_expansions(k0,p,ellp,mdash,n):
    r"""
    Returns a list e of q-expansions, and the Eisenstein series `E_{p-1} = 1 +
    \dots`, all modulo `(p^\text{mdash},q^\text{ellp})`. The list e contains
    the elements `e_{i,s}` in the Katz expansions basis in Step 3 of Algorithm
    1 in [AGBL]_ when one takes as input to that algorithm p,m and k and define
    ``k0``, ``mdash``, n, ``ellp = ell*p`` as in Step 1.

    INPUT:

    - ``k0`` -- integer in range 0 to p-1.
    - ``p`` -- prime at least 5.
    - ``ellp,mdash,n`` -- positive integers.

    OUTPUT:

    - list of q-expansions and the Eisenstein series E_{p-1} modulo
      `(p^\text{mdash},q^\text{ellp})`.

    EXAMPLES::

        sage: from sage.modular.overconvergent.hecke_series import katz_expansions
        sage: katz_expansions(0,5,10,3,4)
        ([1 + O(q^10), q + 6*q^2 + 27*q^3 + 98*q^4 + 65*q^5 + 37*q^6 + 81*q^7 + 85*q^8 + 62*q^9 + O(q^10)],
        1 + 115*q + 35*q^2 + 95*q^3 + 20*q^4 + 115*q^5 + 105*q^6 + 60*q^7 + 25*q^8 + 55*q^9 + O(q^10))
    """
    S = Zmod(p**mdash)

    Ep1 = eisenstein_series_qexp(p-1, ellp, K=S, normalization="constant")
    E4 =  eisenstein_series_qexp(4,   ellp, K=S, normalization="constant")
    E6 =  eisenstein_series_qexp(6,   ellp, K=S, normalization="constant")

    delta = delta_qexp(ellp, K=S)
    h = delta / E6**2
    hj = delta.parent()(1)
    e = []

    # We compute negative powers of E_(p-1) successively (this saves a great
    # deal of time). The effect is that Ep1mi = Ep1 ** (-i).
    Ep1m1 = ~Ep1
    Ep1mi = 1
    for i in xrange(0,n+1):
        Wi,hj = compute_Wi(k0 + i*(p-1),p,h,hj,E4,E6)
        for bis in Wi:
            eis = p**floor(i/(p+1)) * Ep1mi * bis
            e.append(eis)
        Ep1mi = Ep1mi * Ep1m1

    return e,Ep1

def reduced_binary_form1(self):
    r"""
    Reduce the form `ax^2 + bxy+cy^2` to satisfy the reduced condition `|b| \le
    a \le c`, with `b \ge 0` if `a = c`. This reduction occurs within the
    proper class, so all transformations are taken to have determinant 1.

    EXAMPLES::

        sage: QuadraticForm(ZZ,2,[5,5,2]).reduced_binary_form1()
        (
        Quadratic form in 2 variables over Integer Ring with coefficients:
        [ 2 -1 ]
        [ * 2 ]                                                            ,
        <BLANKLINE>
        [ 0 -1]
        [ 1  1]
        )
    """
    if self.dim() != 2:
        raise TypeError("This must be a binary form for now...")

    R = self.base_ring()
    interior_reduced_flag = False
    Q = deepcopy(self)
    M = matrix(R, 2, 2, [1,0,0,1])

    while not interior_reduced_flag:
        interior_reduced_flag = True

        ## Arrange for a <= c
        if Q[0,0] > Q[1,1]:
            M_new = matrix(R,2,2,[0, -1, 1, 0])
            Q = Q(M_new)
            M = M * M_new
            interior_reduced_flag = False
            #print "A"

        ## Arrange for |b| <= a
        if abs(Q[0,1]) > Q[0,0]:
            r = R(floor(round(Q[0,1]/(2*Q[0,0]))))
            M_new = matrix(R,2,2,[1, -r, 0, 1])
            Q = Q(M_new)
            M = M * M_new
            interior_reduced_flag = False
            #print "B"

    return Q, M

def isqrt(x):
    """
    Returns an integer square root, i.e., the floor of a square root.

    EXAMPLES::

        sage: isqrt(10)
        3
        sage: isqrt(10r)
        3
    """
    try:
        return x.isqrt()
    except AttributeError:
        from sage.functions.all import floor
        n = sage.rings.integer.Integer(floor(x))
        return n.isqrt()

def hecke_series_degree_bound(p,N,k,m):
    r"""
    Returns the ``Wan bound`` on the degree of the characteristic series of the
    Atkin operator on p-adic overconvergent modular forms of level
    `\Gamma_0(N)` and weight k when reduced modulo `p^m`. This bound depends
    only upon p, `k \pmod{p-1}`, and N. It uses Lemma 3.1 in [DW]_.

    INPUT:

    - ``p`` -- prime at least 5.
    - ``N`` -- positive integer not divisible by p.
    - ``k`` -- even integer.
    - ``m`` -- positive integer.

    OUTPUT:

    A non-negative integer.

    EXAMPLES::

        sage: from sage.modular.overconvergent.hecke_series import hecke_series_degree_bound
        sage: hecke_series_degree_bound(13,11,100,5)
        39

    REFERENCES:

    .. [DW] Daqing Wan, "Dimension variation of classical and p-adic modular
       forms", Invent. Math. 133, (1998) 449-463.
    """
    k0 = k % (p-1)
    ds = [dimension_modular_forms(N, k0)]
    ms = [ds[0]]
    sum = 0
    u = 1

    ord = 0
    while ord < m:
        ds.append(dimension_modular_forms(N,k0 + u*(p-1)))
        ms.append(ds[u] - ds[u-1])
        sum = sum + u*ms[u]
        ord = floor(((p-1)/(p+1))*sum - ds[u])
        u = u + 1

    return (ds[u-1] - 1)

    def plot_lattice(self):
        r"""
        Plot the lattice (i.e. its points in the cut-off bounds of ``self``).

        OUTPUT:

        - a plot.

        EXAMPLES::

            sage: from sage.geometry.toric_plotter import ToricPlotter
            sage: tp = ToricPlotter(dict(), 2)
            sage: print tp.plot_lattice()
            Graphics object consisting of 1 graphics primitive
        """
        if not self.show_lattice:
            # Plot the origin anyway, otherwise rays/generators may look ugly.
            return self.plot_points([self.origin])
        d = self.dimension
        extra_options = self.extra_options
        if d == 1:
            points = ((x, 0) for x in range(ceil(self.xmin), floor(self.xmax) + 1))
        elif d == 2:
            points = (
                (x, y)
                for x in range(ceil(self.xmin), floor(self.xmax) + 1)
                for y in range(ceil(self.ymin), floor(self.ymax) + 1)
            )
        elif d == 3:
            points = (
                (x, y, z)
                for x in range(ceil(self.xmin), floor(self.xmax) + 1)
                for y in range(ceil(self.ymin), floor(self.ymax) + 1)
                for z in range(ceil(self.zmin), floor(self.zmax) + 1)
            )
        if self.mode == "round":
            r = 1.01 * self.radius  # To make sure integer values work OK.
            points = (pt for pt in points if vector(pt).norm() <= r)
        f = self.lattice_filter
        if f is not None:
            points = (pt for pt in points if f(pt))
        return self.plot_points(tuple(points))

def higher_level_UpGj(p,N,klist,m,modformsring,bound):
    r"""
    Returns a list ``[A_k]`` of square matrices over ``IntegerRing(p^m)``
    parameterised by the weights k in ``klist``. The matrix `A_k` is the finite
    square matrix which occurs on input p,k,N and m in Step 6 of Algorithm 2 in
    [AGBL]_. Notational change from paper: In Step 1 following Wan we defined
    j by `k = k_0 + j(p-1)` with `0 \le k_0 < p-1`. Here we replace j by
    ``kdiv`` so that we may use j as a column index for matrices.)

    INPUT:

    - ``p`` -- prime at least 5.
    - ``N`` -- integer at least 2 and not divisible by p (level).
    - ``klist`` -- list of integers all congruent modulo (p-1) (the weights).
    - ``m`` -- positive integer.
    - ``modformsring`` -- True or False.
    - ``bound`` -- (even) positive integer.

    OUTPUT:

    - list of square matrices.

    EXAMPLES::

        sage: from sage.modular.overconvergent.hecke_series import higher_level_UpGj
        sage: higher_level_UpGj(5,3,[4],2,true,6)
        [
        [ 1  0  0  0  0  0]
        [ 0  1  0  0  0  0]
        [ 0  7  0  0  0  0]
        [ 0  5 10 20  0  0]
        [ 0  7 20  0 20  0]
        [ 0  1 24  0 20  0]
        ]

    """
    t = cputime()
    # Step 1

    k0 = klist[0] % (p-1)
    n = floor(((p+1)/(p-1)) * (m+1))
    elldash = compute_elldash(p,N,k0,n)
    elldashp = elldash*p
    mdash = m + ceil(n/(p+1))

    verbose("done step 1",t)
    t = cputime()
    # Steps 2 and 3

    e,Ep1 = higher_level_katz_exp(p,N,k0,m,mdash,elldash,elldashp,modformsring,bound)
    ell = dimension(transpose(e)[0].parent())
    S = e[0,0].parent()

    verbose("done steps 2+3", t)
    t = cputime()
    # Step 4

    R = Ep1.parent()
    G = compute_G(p, Ep1)
    Alist = []

    verbose("done step 4a", t)
    t = cputime()
    for k in klist:
        k = ZZ(k) # convert to sage integer
        kdiv = k // (p-1)
        Gkdiv = G**kdiv

        T = matrix(S,ell,elldash)
        for i in xrange(ell):
            ei = R(e[i].list())
            Gkdivei = Gkdiv*ei; # act by G^kdiv
            for j in xrange(0, elldash):
                T[i,j] = Gkdivei[p*j]

        verbose("done steps 4b and 5", t)
        t = cputime()

        # Step 6: solve T = AE using fact E is upper triangular.
        # Warning: assumes that T = AE (rather than pT = AE) has
        # a solution over Z/(p^mdash). This has always been the case in
        # examples computed by the author, see Note 3.1.

        A = matrix(S,ell,ell)
        verbose("solving a square matrix problem of dimension %s" % ell)
        verbose("elldash is %s" % elldash)

        for i in xrange(0,ell):
            Ti = T[i]
            for j in xrange(0,ell):
                ej = Ti.parent()([e[j][l] for l in xrange(0,elldash)])
                ejleadpos = ej.nonzero_positions()[0]
                lj = ZZ(ej[ejleadpos])
                A[i,j] = S(ZZ(Ti[j])/lj)
                Ti = Ti - A[i,j]*ej

        Alist.append(MatrixSpace(Zmod(p**m),ell,ell)(A))
        verbose("done step 6", t)

    return Alist

def higher_level_katz_exp(p,N,k0,m,mdash,elldash,elldashp,modformsring,bound):
    r"""
    Returns a matrix `e` of size ``ell x elldashp`` over the integers modulo
    `p^\text{mdash}`, and the Eisenstein series `E_{p-1} = 1 + .\dots \bmod
    (p^\text{mdash},q^\text{elldashp})`. The matrix e contains the coefficients
    of the elements `e_{i,s}` in the Katz expansions basis in Step 3 of
    Algorithm 2 in [AGBL]_ when one takes as input to that algorithm
    `p`,`N`,`m` and `k` and define ``k0``, ``mdash``, ``n``, ``elldash``,
    ``elldashp = ell*dashp`` as in Step 1.

    INPUT:

    - ``p`` -- prime at least 5.
    - ``N`` -- positive integer at least 2 and not divisible by p (level).
    - ``k0`` -- integer in xrange 0 to p-1.
    - ``m,mdash,elldash,elldashp`` -- positive integers.
    - ``modformsring`` -- True or False.
    - ``bound`` -- positive (even) integer.

    OUTPUT:

    - matrix and q-expansion.

    EXAMPLES::

        sage: from sage.modular.overconvergent.hecke_series import higher_level_katz_exp
        sage: e,Ep1 = higher_level_katz_exp(5,2,0,1,2,4,20,true,6)
        sage: e
        [ 1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]
        [ 0  1 18 23 19  6  9  9 17  7  3 17 12  8 22  8 11 19  1  5]
        [ 0  0  1 11 20 16  0  8  4  0 18 15 24  6 15 23  5 18  7 15]
        [ 0  0  0  1  4 16 23 13  6  5 23  5  2 16  4 18 10 23  5 15]
        sage: Ep1
        1 + 15*q + 10*q^2 + 20*q^3 + 20*q^4 + 15*q^5 + 5*q^6 + 10*q^7 +
        5*q^9 + 10*q^10 + 5*q^11 + 10*q^12 + 20*q^13 + 15*q^14 + 20*q^15 + 15*q^16 +
        10*q^17 + 20*q^18 + O(q^20)
    """
    ordr = 1/(p+1)
    S = Zmod(p**mdash)
    Ep1 = eisenstein_series_qexp(p-1,prec=elldashp,K=S, normalization="constant")

    n = floor(((p+1)/(p-1))*(m+1))
    Wjs = complementary_spaces(N,p,k0,n,mdash,elldashp,elldash,modformsring,bound)

    Basis = []
    for j in xrange(n+1):
        Wj = Wjs[j]
        dimj = len(Wj)
        Ep1minusj = Ep1**(-j)
        for i in xrange(dimj):
            wji = Wj[i]
            b = p**floor(ordr*j) * wji * Ep1minusj
            Basis.append(b)

    # extract basis as a matrix

    ell = len(Basis)
    M = matrix(S,ell,elldashp)
    for i in xrange(ell):
        for j in xrange(elldashp):
            M[i,j] = Basis[i][j]

    ech_form(M,p) # put it into echelon form

    return M,Ep1

def theta_series_degree_2(Q, prec):
    r"""
    Compute the theta series of degree 2 for the quadratic form Q.

    INPUT:

    - ``prec`` -- an integer.

    OUTPUT:

    dictionary, where:

    - keys are `{\rm GL}_2(\ZZ)`-reduced binary quadratic forms (given as triples of
      coefficients)
    - values are coefficients

    EXAMPLES::

        sage: Q2 = QuadraticForm(ZZ, 4, [1,1,1,1, 1,0,0, 1,0, 1])
        sage: S = Q2.theta_series_degree_2(10)
        sage: S[(0,0,2)]
        24
        sage: S[(1,0,1)]
        144
        sage: S[(1,1,1)]
        192

    AUTHORS:

    - Gonzalo Tornaria (2010-03-23)

    REFERENCE:

    - Raum, Ryan, Skoruppa, Tornaria, 'On Formal Siegel Modular Forms'
      (preprint)
    """
    if Q.base_ring() != ZZ:
        raise TypeError("The quadratic form must be integral")
    if not Q.is_positive_definite():
        raise ValueError("The quadratic form must be positive definite")
    try:
        X = ZZ(prec-1)    # maximum discriminant
    except TypeError:
        raise TypeError("prec is not an integer")

    if X < -1:
        raise ValueError("prec must be >= 0")

    if X == -1:
        return {}

    V = ZZ ** Q.dim()
    H = Q.Hessian_matrix()

    t = cputime()
    max = int(floor((X+1)/4))
    v_list = (Q.vectors_by_length(max))        # assume a>0
    v_list = [[V(_) for _ in vs] for vs in v_list]  # coerce vectors into V
    verbose("Computed vectors_by_length" , t)

    # Deal with the singular part
    coeffs = {(0,0,0):ZZ(1)}
    for i in range(1,max+1):
        coeffs[(0,0,i)] = ZZ(2) * len(v_list[i])

    # Now deal with the non-singular part
    a_max = int(floor(sqrt(X/3)))
    for a in range(1, a_max + 1):
        t = cputime()
        c_max = int(floor((a*a + X)/(4*a)))
        for c in range(a, c_max + 1):
            for v1 in v_list[a]:
                v1_H = v1 * H
                def B_v1(v):
                    return v1_H * v2
                for v2 in v_list[c]:
                    b = abs(B_v1(v2))
                    if b <= a and 4*a*c-b*b <= X:
                        qf = (a,b,c)
                        count = ZZ(4) if b == 0 else ZZ(2)
                        coeffs[qf] = coeffs.get(qf, ZZ(0)) + count
        verbose("done a = %d" % a, t)

    return coeffs

 def rangea(t):
     tp0=t*wprime0
     tp1=t*wprime1
     return sorted(list(union(range(floor(xmin-_sage_const_1 -tp0-wprime0),ceil(xmax-tp0)+_sage_const_1 ),range(floor(ymin-_sage_const_1 -tp1-wprime1),ceil(ymax-tp1)+_sage_const_1 ))),key=abs)

        def coproduct_iterated(self, n=1):
            r"""
            Apply ``n`` coproducts to ``self``.

            .. TODO::

                Remove dependency on ``modules_with_basis`` methods.

            EXAMPLES::

                sage: Psi = NonCommutativeSymmetricFunctions(QQ).Psi()
                sage: Psi[2,2].coproduct_iterated(0)
                Psi[2, 2]
                sage: Psi[2,2].coproduct_iterated(2)
                Psi[] # Psi[] # Psi[2, 2] + 2*Psi[] # Psi[2] # Psi[2]
                 + Psi[] # Psi[2, 2] # Psi[] + 2*Psi[2] # Psi[] # Psi[2]
                 + 2*Psi[2] # Psi[2] # Psi[] + Psi[2, 2] # Psi[] # Psi[]

            TESTS::

                sage: p = SymmetricFunctions(QQ).p()
                sage: p[5,2,2].coproduct_iterated()
                p[] # p[5, 2, 2] + 2*p[2] # p[5, 2] + p[2, 2] # p[5]
                 + p[5] # p[2, 2] + 2*p[5, 2] # p[2] + p[5, 2, 2] # p[]
                sage: p([]).coproduct_iterated(3)
                p[] # p[] # p[] # p[]

            ::

                sage: Psi = NonCommutativeSymmetricFunctions(QQ).Psi()
                sage: Psi[2,2].coproduct_iterated(0)
                Psi[2, 2]
                sage: Psi[2,2].coproduct_iterated(3)
                Psi[] # Psi[] # Psi[] # Psi[2, 2] + 2*Psi[] # Psi[] # Psi[2] # Psi[2]
                 + Psi[] # Psi[] # Psi[2, 2] # Psi[] + 2*Psi[] # Psi[2] # Psi[] # Psi[2]
                 + 2*Psi[] # Psi[2] # Psi[2] # Psi[] + Psi[] # Psi[2, 2] # Psi[] # Psi[]
                 + 2*Psi[2] # Psi[] # Psi[] # Psi[2] + 2*Psi[2] # Psi[] # Psi[2] # Psi[]
                 + 2*Psi[2] # Psi[2] # Psi[] # Psi[] + Psi[2, 2] # Psi[] # Psi[] # Psi[]

            ::

                sage: m = SymmetricFunctionsNonCommutingVariables(QQ).m()
                sage: m[[1,3],[2]].coproduct_iterated(2)
                m{} # m{} # m{{1, 3}, {2}} + m{} # m{{1}} # m{{1, 2}}
                 + m{} # m{{1, 2}} # m{{1}} + m{} # m{{1, 3}, {2}} # m{}
                 + m{{1}} # m{} # m{{1, 2}} + m{{1}} # m{{1, 2}} # m{}
                 + m{{1, 2}} # m{} # m{{1}} + m{{1, 2}} # m{{1}} # m{}
                 + m{{1, 3}, {2}} # m{} # m{}
                sage: m[[]].coproduct_iterated(3), m[[1,3],[2]].coproduct_iterated(0)
                (m{} # m{} # m{} # m{}, m{{1, 3}, {2}})
            """
            if n < 0:
                raise ValueError("cannot take fewer than 0 coproduct iterations: %s < 0" % str(n))
            if n == 0:
                return self
            if n == 1:
                return self.coproduct()
            from sage.functions.all import floor, ceil
            from sage.rings.all import Integer

            # Use coassociativity of `\Delta` to perform many coproducts simultaneously.
            fn = floor(Integer(n-1)/2); cn = ceil(Integer(n-1)/2)
            split = lambda a,b: tensor([a.coproduct_iterated(fn), b.coproduct_iterated(cn)])
            return self.coproduct().apply_multilinear_morphism(split)

def berlekamp_welsh(deg, points):
    r"""
    Reconstruct polynomial with Berlekamp-Welsh algorithm.

    INPUT:

    - ``deg``    --  degree of polynomial to reconstruct.
    - ``points`` --  array of points (list of (x,y)-tuples).

    OUTPUT:

    Reconstructed polynomial.

    EXAMPLES::

        sage: from sage.crypto.smc.berlekamp_welsh import berlekamp_welsh
        sage: from sage.rings.finite_rings.constructor import FiniteField
        sage: from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

    Reconstruction with errors::

        sage: order = 2**8
        sage: F = FiniteField(order, 'a')
        sage: P = PolynomialRing(F, 'x')
        sage: n = 7
        sage: deg = 2
        sage: poly = F.fetch_int(42)
        sage: for i in range(1, deg+1): poly += F.random_element() * P.gen()**i

        sage: # evaluate polynomial at different points (shares)
        sage: points = [(F.fetch_int(i), poly(F.fetch_int(i))) for i in range(1, n+1)]
        sage: poly == berlekamp_welsh(deg, points)
        True

        sage: # introduce error
        sage: points[0] = (points[0][0], points[0][1] + F.fetch_int(9))
        sage: poly == berlekamp_welsh(deg, points)
        True

    """
    # check input vector
    F = points[0][0].parent()
    if not F.is_field():
        raise TypeError("points must be of field type.")
    for x, y in points:
        if x.parent() != F or y.parent() != F:
            raise TypeError("all points must be from same field.")
        
    # generate and solve system of linear equations
    from sage.functions.all import floor
    deg_E = floor((len(points) - (deg + 1)) / 2.)
    deg_Q = deg_E + deg
    from sage.matrix.all import Matrix
    from sage.all import vector
    sys_size = deg_Q + 1 + deg_E
    A = Matrix(F, sys_size)
    b = vector(F, sys_size)
    for n, (x, y) in enumerate(points):
        A[n] = ([x**i for i in range(deg_Q+1)]+[-y * x**i for i in range(deg_E)])
        b[n] = (y * x**deg_E)
    QE = A.solve_right(b)

    # reconstruct polynomial
    from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    P = PolynomialRing(F, 'x')
    Q = sum([coeff * P.gen()**i for i, coeff in enumerate(QE[:deg_Q+1])]);
    E = P.gen()**deg_E
    E += sum([coeff * P.gen()**i for i, coeff in enumerate(QE[deg_Q+1:])]);
    P = Q.quo_rem(E)[0]
    return 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 = floor((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

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

def level1_UpGj(p,klist,m):
    r"""
    Returns a list `[A_k]` of square matrices over ``IntegerRing(p^m)``
    parameterised by the weights k in ``klist``. The matrix `A_k` is the finite
    square matrix which occurs on input p,k and m in Step 6 of Algorithm 1 in
    [AGBL]_. Notational change from paper: In Step 1 following Wan we defined
    j by `k = k_0 + j(p-1)` with `0 \le k_0 < p-1`. Here we replace j by
    ``kdiv`` so that we may use j as a column index for matrices.

    INPUT:

    - ``p`` -- prime at least 5.
    - ``klist`` -- list of integers congruent modulo `(p-1)` (the weights).
    - ``m`` -- positive integer.

    OUTPUT:

    - list of square matrices.

    EXAMPLES::

        sage: from sage.modular.overconvergent.hecke_series import level1_UpGj
        sage: level1_UpGj(7,[100],5)
        [
        [    1   980  4802     0     0]
        [    0 13727 14406     0     0]
        [    0 13440  7203     0     0]
        [    0  1995  4802     0     0]
        [    0  9212 14406     0     0]
        ]
    """
    # Step 1
    t = cputime()

    k0 = klist[0] % (p-1)
    n = floor(((p+1)/(p-1)) * (m+1))
    ell = dimension_modular_forms(1, k0 + n*(p-1))
    ellp = ell*p
    mdash = m + ceil(n/(p+1))

    verbose("done step 1", t)
    t = cputime()
    # Steps 2 and 3

    e,Ep1 = katz_expansions(k0,p,ellp,mdash,n)

    verbose("done steps 2+3", t)
    t=cputime()
    # Step 4

    G = compute_G(p, Ep1)
    Alist = []

    verbose("done step 4a", t)
    t=cputime()
    for k in klist:
        k = ZZ(k) # convert to sage integer
        kdiv = k // (p-1)
        Gkdiv = G**kdiv
        u = []
        for i in xrange(0,ell):
            ei = e[i]
            ui = Gkdiv*ei
            u.append(ui)

        verbose("done step 4b", t)
        t = cputime()
        # Step 5 and computation of T in Step 6

        S = e[0][0].parent()
        T = matrix(S,ell,ell)

        for i in xrange(0,ell):
            for j in xrange(0,ell):
                T[i,j] = u[i][p*j]

        verbose("done step 5", t)
        t = cputime()
        # Step 6: solve T = AE using fact E is upper triangular.
        # Warning: assumes that T = AE (rather than pT = AE) has
        # a solution over Z/(p^mdash). This has always been the case in
        # examples computed by the author, see Note 3.1.

        A = matrix(S,ell,ell)
        verbose("solving a square matrix problem of dimension %s" % ell, t)

        for i in xrange(0,ell):
            Ti = T[i]
            for j in xrange(0,ell):
                ej = Ti.parent()([e[j][l] for l in xrange(0,ell)])
                lj = ZZ(ej[j])
                A[i,j] = S(ZZ(Ti[j])/lj)
                Ti = Ti - A[i,j]*ej

        Alist.append(MatrixSpace(Zmod(p**m),ell,ell)(A))
        verbose("done step 6", t)

    return Alist

def reduced_binary_form(self):
    """
    Find a form which is reduced in the sense that no further binary
    form reductions can be done to reduce the original form.

    EXAMPLES::

        sage: QuadraticForm(ZZ,2,[5,5,2]).reduced_binary_form()
        (
        Quadratic form in 2 variables over Integer Ring with coefficients:
        [ 2 -1 ]
        [ * 2 ]                                                            ,
        <BLANKLINE>
        [ 0 -1]
        [ 1  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


    while not interior_reduced_flag:
        interior_reduced_flag = True

        #print Q

        ## Arrange for (weakly) increasing diagonal entries
        for i in range(n):
            for j in range(i+1,n):
                if Q[i,i] > Q[j,j]:
                    M_new = matrix(R,n,n)
                    for k in range(n):
                        M_new[k,k] = 1
                    M_new[i,j] = -1
                    M_new[j,i] = 1
                    M_new[i,i] = 0
                    M_new[j,j] = 1

                    Q = Q(M_new)
                    M = M * M_new
                    interior_reduced_flag = False
                    #print "A"

                ## Arrange for |b| <= a
                if abs(Q[i,j]) > Q[i,i]:
                    r = R(floor(round(Q[i,j]/(2*Q[i,i]))))

                    M_new = matrix(R,n,n)
                    for k in range(n):
                        M_new[k,k] = 1
                    M_new[i,j] = -r

                    Q = Q(M_new)
                    M = M * M_new
                    interior_reduced_flag = False
                    #print "B"

    return Q, M

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


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