def __init__(self, R, elements):
        """
        Initialize ``self``.

        EXAMPLES::

            sage: R.<x,y,z> = QQ[]
            sage: K = KoszulComplex(R, [x,y])
            sage: TestSuite(K).run()
        """
        # Generate the differentials
        self._elements = elements
        n = len(elements)
        I = range(n)
        diff = {}
        zero = R.zero()
        for i in I:
            M = matrix(R, binomial(n,i), binomial(n,i+1), zero)
            j = 0
            for comb in itertools.combinations(I, i+1):
                for k,val in enumerate(comb):
                    r = rank(comb[:k] + comb[k+1:], n, False)
                    M[r,j] = (-1)**k * elements[val]
                j += 1
            M.set_immutable()
            diff[i+1] = M
        diff[0] = matrix(R, 0, 1, zero)
        diff[0].set_immutable()
        diff[n+1] = matrix(R, 1, 0, zero)
        diff[n+1].set_immutable()
        ChainComplex_class.__init__(self, ZZ, ZZ(-1), R, diff)

    def cardinality(self):
        """
        EXAMPLES::

            sage: IntegerVectors(3,3, min_part=1).cardinality()
            1
            sage: IntegerVectors(5,3, min_part=1).cardinality()
            6
            sage: IntegerVectors(13, 4, min_part=2, max_part=4).cardinality()
            16
        """
        if not self._constraints:
            if self.k < 0:
                return +infinity
            if self.n >= 0:
                return binomial(self.n+self.k-1,self.n)
            else:
                return 0
        else:
            if len(self._constraints) == 1 and 'max_part' in self._constraints and self._constraints['max_part'] != infinity:
                m = self._constraints['max_part']
                if m >= self.n:
                    return binomial(self.n+self.k-1,self.n)
                else: #do by inclusion / exclusion on the number
                      #i of parts greater than m
                    return sum( [(-1)**i * binomial(self.n+self.k-1-i*(m+1), self.k-1)*binomial(self.k,i) for i in range(0, self.n/(m+1)+1)])
            else:
                return super(IntegerVectors_nkconstraints, self).cardinality()

    def cardinality(self):
        r"""
        Return the number of Baxter permutations of size ``self._n``.

        For any positive integer `n`, the number of Baxter
        permutations of size `n` equals

        .. MATH::

            \sum_{k=1}^n \dfrac
            {\binom{n+1}{k-1} \binom{n+1}{k} \binom{n+1}{k+1}}
            {\binom{n+1}{1} \binom{n+1}{2}} .

        This is :oeis:`A001181`.

        EXAMPLES::

            sage: [BaxterPermutations(n).cardinality() for n in range(13)]
            [1, 1, 2, 6, 22, 92, 422, 2074, 10754, 58202, 326240, 1882960, 11140560]

            sage: BaxterPermutations(3r).cardinality()
            6
            sage: parent(_)
            Integer Ring
        """
        if self._n == 0:
            return 1
        from sage.arith.all import binomial
        return sum((binomial(self._n + 1, k) *
                    binomial(self._n + 1, k + 1) *
                    binomial(self._n + 1, k + 2)) //
                   ((self._n + 1) * binomial(self._n + 1, 2))
                   for k in range(self._n))

def from_rank(r, n, k):
    r"""
    Returns the combination of rank ``r`` in the subsets of
    ``range(n)`` of size ``k`` when listed in lexicographic order.

    The algorithm used is based on combinadics and James McCaffrey's
    MSDN article. See: :wikipedia:`Combinadic`

    EXAMPLES::

        sage: import sage.combinat.combination as combination
        sage: combination.from_rank(0,3,0)
        ()
        sage: combination.from_rank(0,3,1)
        (0,)
        sage: combination.from_rank(1,3,1)
        (1,)
        sage: combination.from_rank(2,3,1)
        (2,)
        sage: combination.from_rank(0,3,2)
        (0, 1)
        sage: combination.from_rank(1,3,2)
        (0, 2)
        sage: combination.from_rank(2,3,2)
        (1, 2)
        sage: combination.from_rank(0,3,3)
        (0, 1, 2)
    """
    if k < 0:
        raise ValueError("k must be > 0")
    if k > n:
        raise ValueError("k must be <= n")

    a = n
    b = k
    x = binomial(n, k) - 1 - r  # x is the 'dual' of m
    comb = [None] * k

    for i in xrange(k):
        comb[i] = _comb_largest(a, b, x)
        x = x - binomial(comb[i], b)
        a = comb[i]
        b = b - 1

    for i in xrange(k):
        comb[i] = (n - 1) - comb[i]

    return tuple(comb)

    def unrank(self, r):
        """
        Return the subset of ``s`` that has rank ``k``.

        EXAMPLES::

            sage: Subsets(3).unrank(0)
            {}
            sage: Subsets([2,4,5]).unrank(1)
            {2}
            sage: Subsets([1,2,3]).unrank(257)
            Traceback (most recent call last):
            ...
            IndexError: index out of range

        """
        r = Integer(r)
        if r >= self.cardinality() or r < 0:
            raise IndexError("index out of range")
        else:
            k = ZZ_0
            n = self._s.cardinality()
            bin = Integer(1)
            while r >= bin:
                r -= bin
                k += 1
                bin = binomial(n,k)
            return self.element_class([self._s.unrank(i) for i in combination.from_rank(r, n, k)])

    def rank(self, sub):
        """
        Return the rank of ``sub`` as a subset of ``s``.

        EXAMPLES::

            sage: Subsets(3).rank([])
            0
            sage: Subsets(3).rank([1,2])
            4
            sage: Subsets(3).rank([1,2,3])
            7
            sage: Subsets(3).rank([2,3,4])
            Traceback (most recent call last):
            ...
            ValueError: {2, 3, 4} is not a subset of {1, 2, 3}
        """
        if sub not in Sets():
            ssub = Set(sub)
            if len(sub) != len(ssub):
                raise ValueError("repeated elements in {}".format(sub))
            sub = ssub

        try:
            index_list = sorted(self._s.rank(x) for x in sub)
        except (ValueError,IndexError):
            raise ValueError("{} is not a subset of {}".format(
                    Set(sub), self._s))

        n = self._s.cardinality()
        r = sum(binomial(n,i) for i in range(len(index_list)))
        return r + combination.rank(index_list,n)

    def cardinality(self):
        r"""
        Return the number of words in the shuffle product
        of ``w1`` and ``w2``.

        This is understood as a multiset cardinality, not as a
        set cardinality; it does not count the distinct words only.

        It is given by `\binom{l_1+l_2}{l_1}`, where `l_1` is the
        length of ``w1`` and where `l_2` is the length of ``w2``.

        EXAMPLES::

            sage: from sage.combinat.words.shuffle_product import ShuffleProduct_w1w2
            sage: w, u = map(Words("abcd"), ["ab", "cd"])
            sage: S = ShuffleProduct_w1w2(w,u)
            sage: S.cardinality()
            6

            sage: w, u = map(Words("ab"), ["ab", "ab"])
            sage: S = ShuffleProduct_w1w2(w,u)
            sage: S.cardinality()
            6
        """
        return binomial(self._w1.length()+self._w2.length(), self._w1.length())

    def rank(self, x):
        """
        Returns the position of a given element.

        INPUT:

        - ``x`` - a list with ``sum(x) == n`` and ``len(x) == k``

        TESTS::

            sage: IV = IntegerVectors(4,5) 
            sage: range(IV.cardinality()) == [IV.rank(x) for x in IV]
            True
        """

        if x not in self:
            raise ValueError("argument is not a member of IntegerVectors(%d,%d)" % (self.n, self.k))

        n = self.n
        k = self.k

        r = 0
        for i in range(k-1):
          k -= 1
          n -= x[i]
          r += binomial(k+n-1,k)

        return r

def upper_bound(min_length, max_length, floor, ceiling, min_slope, max_slope):
    """
    Compute a coarse upper bound on the size of a vector satisfying the
    constraints.

    TESTS::

        sage: import sage.combinat.integer_list_old as integer_list
        sage: f = lambda x: lambda i: x
        sage: integer_list.upper_bound(0,4,f(0), f(1),-infinity,infinity)
        4
        sage: integer_list.upper_bound(0, infinity, f(0), f(1), -infinity, infinity)
        inf
        sage: integer_list.upper_bound(0, infinity, f(0), f(1), -infinity, -1)
        1
        sage: integer_list.upper_bound(0, infinity, f(0), f(5), -infinity, -1)
        15
        sage: integer_list.upper_bound(0, infinity, f(0), f(5), -infinity, -2)
        9
    """
    from sage.functions.all import floor as flr
    if max_length < float('inf'):
        return sum( [ ceiling(j) for j in range(max_length)] )
    elif max_slope < 0 and ceiling(1) < float('inf'):
        maxl = flr(-ceiling(1)/max_slope)
        return ceiling(1)*(maxl+1) + binomial(maxl+1,2)*max_slope
    #FIXME: only checking the first 10000 values, but that should generally
    #be enough
    elif [ceiling(j) for j in range(10000)] == [0]*10000:
        return 0
    else:
        return float('inf')

def log_gamma_binomial(p, gamma, z, n, M):
    r"""
    Return the list of coefficients in the power series
    expansion (up to precision `M`) of `\binom{\log_p(z)/\log_p(\gamma)}{n}`

    INPUT:

    - ``p`` --  prime
    - ``gamma`` -- topological generator, e.g. `1+p`
    - ``z`` -- variable
    - ``n`` -- nonnegative integer
    - ``M`` -- precision

    OUTPUT:

    The list of coefficients in the power series expansion of
    `\binom{\log_p(z)/\log_p(\gamma)}{n}`

    EXAMPLES::

        sage: R.<z> = QQ['z']
        sage: from sage.modular.pollack_stevens.padic_lseries import log_gamma_binomial
        sage: log_gamma_binomial(5,1+5,z,2,4)
        [0, -3/205, 651/84050, -223/42025]
        sage: log_gamma_binomial(5,1+5,z,3,4)
        [0, 2/205, -223/42025, 95228/25845375]
    """
    L = sum([ZZ(-1) ** j / j * z ** j for j in range(1, M)])  # log_p(1+z)
    loggam = L / (L(gamma - 1))  # log_{gamma}(1+z)= log_p(1+z)/log_p(gamma)
    return z.parent()(binomial(loggam, n)).truncate(M).list()

    def _basic_integral(self, a, j):
        r"""
        Return `\int_{a+pZ_p} (z-{a})^j d\Phi(0-infty)`
        -- see formula [Pollack-Stevens, sec 9.2]

        INPUT:

        - ``a`` -- integer in range(p)
        - ``j`` -- integer in range(self.symbol().precision_relative())

        EXAMPLES::

            sage: from sage.modular.pollack_stevens.padic_lseries import pAdicLseries
            sage: E = EllipticCurve('11a3')
            sage: L = E.padic_lseries(5, implementation="pollackstevens", precision=4) #long time
            sage: L._basic_integral(1,2) # long time
            2*5^2 + 5^3 + O(5^4)
        """
        symb = self.symbol()
        M = symb.precision_relative()
        if j > M:
            raise PrecisionError("Too many moments requested")
        p = self.prime()
        ap = symb.Tq_eigenvalue(p)
        D = self._quadratic_twist
        ap = ap * kronecker(D, p)
        K = pAdicField(p, M)
        symb_twisted = symb.evaluate_twisted(a, D)
        return (
            sum(
                binomial(j, r) * ((a - ZZ(K.teichmuller(a))) ** (j - r)) * (p ** r) * symb_twisted.moment(r)
                for r in range(j + 1)
            )
            / ap
        )

    def _induced_flags(self, n, tg, type_edges):

        flag_counts = {}
        flags = []
        total = 0

        for p in Tuples([0, 1], binomial(n, 2) - binomial(tg.n, 2)):

            edges = list(type_edges)

            c = 0
            for i in range(tg.n + 1, n + 1):
                for j in range(1, i):
                    if p[c] == 0:
                        edges.append((i, j))
                    else:
                        edges.append((j, i))
                    c += 1

            ig = ThreeGraphFlag()
            ig.n = n
            ig.t = tg.n

            for s in Combinations(range(1, n + 1), 3):
                if self._variant:
                    if ((s[1], s[0]) in edges and (s[0], s[2]) in edges) or (
                            (s[2], s[0]) in edges and (s[0], s[1]) in edges):
                        ig.add_edge(s)
                else:
                    if ((s[0], s[1]) in edges and (s[1], s[2]) in edges and (s[2], s[0]) in edges) or (
                            (s[0], s[2]) in edges and (s[2], s[1]) in edges and (s[1], s[0]) in edges):
                        ig.add_edge(s)

            it = ig.induced_subgraph(range(1, tg.n + 1))
            if tg.is_labelled_isomorphic(it):
                ig.make_minimal_isomorph()

                ghash = hash(ig)
                if ghash in flag_counts:
                    flag_counts[ghash] += 1
                else:
                    flags.append(ig)
                    flag_counts[ghash] = 1

            total += 1

        return [(f, flag_counts[hash(f)] / Integer(total)) for f in flags]

    def unrank(self, r):
        """
        EXAMPLES::

            sage: c = Combinations([1,2,3])
            sage: c.list() == map(c.unrank, range(c.cardinality()))
            True
        """
        k = 0
        n = len(self.mset)
        b = binomial(n, k)
        while r >= b:
            r -= b
            k += 1
            b = binomial(n,k)

        return [self.mset[i] for i in from_rank(r, n, k)]

    def cardinality(self):
        r"""
        Returns the number of subwords of w of length k.

        EXAMPLES::

            sage: Subwords([1,2,3], 2).cardinality()
            3
        """
        return arith.binomial(Integer(len(self._w)), self._k)

    def cardinality(self):
        """
        Return the cardinality of ``self``.

        EXAMPLES::

            sage: IntegerVectors(3, 3, min_part=1).cardinality()
            1
            sage: IntegerVectors(5, 3, min_part=1).cardinality()
            6
            sage: IntegerVectors(13, 4, max_part=4).cardinality()
            20
            sage: IntegerVectors(k=4, max_part=3).cardinality()
            256
            sage: IntegerVectors(k=3, min_part=2, max_part=4).cardinality()
            27
            sage: IntegerVectors(13, 4, min_part=2, max_part=4).cardinality()
            16
        """
        if self.k is None:
            if self.n is None:
                return PlusInfinity()
            if ('max_length' not in self.constraints
                    and self.constraints.get('min_part', 0) <= 0):
                return PlusInfinity()
        elif ('max_part' in self.constraints
                and self.constraints['max_part'] != PlusInfinity()):
            if (self.n is None and len(self.constraints) == 2
                    and 'min_part' in self.constraints
                    and self.constraints['min_part'] >= 0):
                num = self.constraints['max_part'] - self.constraints['min_part'] + 1
                return Integer(num ** self.k)
            if len(self.constraints) == 1:
                m = self.constraints['max_part']
                if self.n is None:
                    return Integer((m + 1) ** self.k)
                if m >= self.n:
                    return Integer(binomial(self.n + self.k - 1, self.n))
                # do by inclusion / exclusion on the number
                # i of parts greater than m
                return Integer(sum( (-1)**i * binomial(self.n+self.k-1-i*(m+1), self.k-1) \
                    * binomial(self.k,i) for i in range(self.n/(m+1)+1) ))
        return ZZ.sum(ZZ.one() for x in self)

def tdesign_params(t, v, k, L):
    """
    Return the design's parameters: `(t, v, b, r , k, L)`. Note that `t` must be
    given.

    EXAMPLES::

        sage: BD = BlockDesign(7,[[0,1,2],[0,3,4],[0,5,6],[1,3,5],[1,4,6],[2,3,6],[2,4,5]])
        sage: from sage.combinat.designs.block_design import tdesign_params
        sage: tdesign_params(2,7,3,1)
        (2, 7, 7, 3, 3, 1)
    """
    x = binomial(v, t)
    y = binomial(k, t)
    b = divmod(L * x, y)[0]
    x = binomial(v-1, t-1)
    y = binomial(k-1, t-1)
    r = integer_floor(L * x/y)
    return (t, v, b, r, k, L)

    def cardinality(self):
        """
        Returns the number of multichoices of k things from a list of n
        things.

        EXAMPLES::

            sage: MultichooseNK(3,2).cardinality()
            6
        """
        n,k = self._n, self._k
        return binomial(n+k-1,k)

    def by_taylor_expansion(self, fs, k, is_integral=False) :
        r"""
        We combine the theta decomposition and the heat operator as in [Sko].
        This yields a bijections of Jacobi forms of weight `k` and
        `M_k \times S_{k+2} \times .. \times S_{k+2m}`.
        """
        ## we introduce an abbreviations
        if is_integral :
            PS = self.integral_power_series_ring()
        else :
            PS = self.power_series_ring()
            
        if not len(fs) == self.__precision.jacobi_index() + 1 :
            raise ValueError("fs must be a list of m + 1 elliptic modular forms or their fourier expansion")
        
        qexp_prec = self._qexp_precision()
        if qexp_prec is None : # there are no forms below the precision
            return dict()
        
        f_divs = dict()
        for (i, f) in enumerate(fs) :
            f_divs[(i, 0)] = PS(f(qexp_prec), qexp_prec)
        
        if self.__precision.jacobi_index() == 1 :
            return self._by_taylor_expansion_m1(f_divs, k, is_integral)
        
        for i in xrange(self.__precision.jacobi_index() + 1) :
            for j in xrange(1, self.__precision.jacobi_index() - i + 1) :
                f_divs[(i,j)] = f_divs[(i, j - 1)].derivative().shift(1)
            
        phi_divs = list()
        for i in xrange(self.__precision.jacobi_index() + 1) :
            ## This is the formula in Skoruppas thesis. He uses d/ d tau instead of d / dz which yields
            ## a factor 4 m
            phi_divs.append( sum( f_divs[(j, i - j)] * (4 * self.__precision.jacobi_index())**i
                                  * binomial(i,j) / 2**i#2**(self.__precision.jacobi_index() - i + 1)
                                  * prod(2*(i - l) + 1 for l in xrange(1, i))
                                  / factorial(i + k + j - 1)
                                  * factorial(2*self.__precision.jacobi_index() + k - 1) 
                                  for j in xrange(i + 1) ) )
            
        phi_coeffs = dict()
        for r in xrange(self.__precision.jacobi_index() + 1) :
            series = sum( map(operator.mul, self._theta_factors()[r], phi_divs) )
            series = self._eta_factor() * series

            for n in xrange(qexp_prec) :
                phi_coeffs[(n, r)] = series[n]

        return phi_coeffs

    def rank(self, x):
        """
        EXAMPLES::

            sage: c = Combinations([1,2,3])
            sage: range(c.cardinality()) == map(c.rank, c)
            True
        """
        x = [self.mset.index(_) for _ in x]
        r = 0
        n = len(self.mset)
        for i in range(len(x)):
            r += binomial(n, i)
        r += rank(x, n)
        return r

def s(f,k):
    """Given f monic of degree n with distinct roots, returns the monic
    polynomial of degree binomial(n,k) whose roots are all products of
    k distinct roots of f.
    """
    n = f.degree()
    x = f.variables()[0]
    if k==0:
        return x-1
    if k==1:
        return f
    c = (-1)**n * f(0)
    if k==n:
        return x-c
    if k>n-k: # use s(f,n-k)
        g = s(f,n-k)
        from sage.arith.all import binomial
        return ((-x)**binomial(n,k) * g(c/x) / c**binomial(n-1,k)).numerator()

    if k==2:
        return (star(f,f)//r(f,2)).nth_root(2)
    if k==3:
        f2 = s(f,2)
        return (star(f2,f)*r(f,3) // star(r(f,2),f)).nth_root(3)

    fkn = fkd = 1
    for j in range(1,k+1):
        g = star(r(f,j), s(f,k-j))
        if j%2:
            fkn *= g
        else:
            fkd *= g

    fk = fkn//fkd
    assert fk*fkd==fkn
    return fk.nth_root(k)