def _calc_bernoulli(n):
     s = 0
     a = int(binomial(n + 3, n - 6))
     for j in range(1, n//6 + 1):
         s += a * bernoulli(n - 6*j)
         # Avoid computing each binomial coefficient from scratch
         a *= _product(n - 6 - 6*j + 1, n - 6*j)
         a //= _product(6*j + 4, 6*j + 9)
     if n % 6 == 4:
         s = -Rational(n + 3, 6) - s
     else:
         s = Rational(n + 3, 3) - s
     return s / binomial(n + 3, n)

    def rule(n, k):
        coeff, rewrite = S.One, False

        cn, _n = n.as_coeff_Add()

        if _n and cn.is_Integer and cn:
            coeff *= _rf(_n + 1, cn)/_rf(_n - k + 1, cn)
            rewrite = True
            n = _n

        # this sort of binomial has already been removed by
        # rising factorials but is left here in case the order
        # of rule application is changed
        if k.is_Add:
            ck, _k = k.as_coeff_Add()
            if _k and ck.is_Integer and ck:
                coeff *= _rf(n - ck - _k + 1, ck)/_rf(_k + 1, ck)
                rewrite = True
                k = _k

        if count_ops(k) > count_ops(n - k):
            rewrite = True
            k = n - k

        if rewrite:
            return coeff*binomial(n, k)

    def _eval_rewrite_as_Sum(self, arg):
        if arg.is_even:
            k = Dummy("k", integer=True)
            j = Dummy("j", integer=True)
            n = self.args[0] / 2
            Em = (S.ImaginaryUnit * C.Sum(C.Sum(binomial(k, j) * ((-1)**j * (k - 2*j)**(2*n + 1)) /
                  (2**k*S.ImaginaryUnit**k * k), (j, 0, k)), (k, 1, 2*n + 1)))

            return Em

 def get_size(self):
     r"""
     Returns
     -------
     size: int
         The size of set T. Set T is the set of all possible
         monomials of the n variables for degree equal to the
         degree_m
     """
     return binomial(self.degree_m + self.n - 1, self.n - 1)

def _stirling1(n, k):
    if n == k == 0:
        return S.One
    if 0 in (n, k):
        return S.Zero
    n1 = n - 1

    # some special values
    if n == k:
        return S.One
    elif k == 1:
        return factorial(n1)
    elif k == n1:
        return binomial(n, 2)
    elif k == n - 2:
        return (3*n - 1)*binomial(n, 3)/4
    elif k == n - 3:
        return binomial(n, 2)*binomial(n, 4)

    # general recurrence
    return n1*_stirling1(n1, k) + _stirling1(n1, k - 1)

def _stirling2(n, k):
    if n == k == 0:
        return S.One
    if 0 in (n, k):
        return S.Zero
    n1 = n - 1

    # some special values
    if k == n1:
        return binomial(n, 2)
    elif k == 2:
        return 2**n1 - 1

    # general recurrence
    return k*_stirling2(n1, k) + _stirling2(n1, k - 1)

    def eval(cls, n, sym=None):
        if n.is_Number:
            if n.is_Integer and n.is_nonnegative:
                if n is S.Zero:
                    return S.One
                elif n is S.One:
                    if sym is None:
                        return -S.Half
                    else:
                        return sym - S.Half
                # Bernoulli numbers
                elif sym is None:
                    if n.is_odd:
                        return S.Zero
                    n = int(n)
                    # Use mpmath for enormous Bernoulli numbers
                    if n > 500:
                        p, q = bernfrac(n)
                        return Rational(int(p), int(q))
                    case = n % 6
                    highest_cached = cls._highest[case]
                    if n <= highest_cached:
                        return cls._cache[n]
                    # To avoid excessive recursion when, say, bernoulli(1000) is
                    # requested, calculate and cache the entire sequence ... B_988,
                    # B_994, B_1000 in increasing order
                    for i in range(highest_cached + 6, n + 6, 6):
                        b = cls._calc_bernoulli(i)
                        cls._cache[i] = b
                        cls._highest[case] = i
                    return b
                # Bernoulli polynomials
                else:
                    n, result = int(n), []
                    for k in range(n + 1):
                        result.append(binomial(n, k)*cls(k)*sym**(n - k))
                    return Add(*result)
            else:
                raise ValueError("Bernoulli numbers are defined only"
                                 " for nonnegative integer indices.")

        if sym is None:
            if n.is_odd and (n - 1).is_positive:
                return S.Zero

    def f(rv):
        if not rv.is_Mul:
            return rv
        rvd = rv.as_powers_dict()
        nd_fact_args = [[], []] # numerator, denominator

        for k in rvd:
            if isinstance(k, factorial) and rvd[k].is_Integer:
                if rvd[k].is_positive:
                    nd_fact_args[0].extend([k.args[0]]*rvd[k])
                else:
                    nd_fact_args[1].extend([k.args[0]]*-rvd[k])
                rvd[k] = 0
        if not nd_fact_args[0] or not nd_fact_args[1]:
            return rv

        hit = False
        for m in range(2):
            i = 0
            while i < len(nd_fact_args[m]):
                ai = nd_fact_args[m][i]
                for j in range(i + 1, len(nd_fact_args[m])):
                    aj = nd_fact_args[m][j]

                    sum = ai + aj
                    if sum in nd_fact_args[1 - m]:
                        hit = True

                        nd_fact_args[1 - m].remove(sum)
                        del nd_fact_args[m][j]
                        del nd_fact_args[m][i]

                        rvd[binomial(sum, ai if count_ops(ai) <
                                count_ops(aj) else aj)] += (
                                -1 if m == 0 else 1)
                        break
                else:
                    i += 1

        if hit:
            return Mul(*([k**rvd[k] for k in rvd] + [factorial(k)
                    for k in nd_fact_args[0]]))/Mul(*[factorial(k)
                    for k in nd_fact_args[1]])
        return rv

    def eval(cls, n, alpha, x):
        # L_{n}^{0}(x)  --->  L_{n}(x)
        if alpha == S.Zero:
            return laguerre(n, x)

        if not n.is_Number:
            # We can evaluate for some special values of x
            if x == S.Zero:
                return binomial(n + alpha, alpha)
            elif x == S.Infinity and n > S.Zero:
                return S.NegativeOne**n * S.Infinity
            elif x == S.NegativeInfinity and n > S.Zero:
                return S.Infinity
        else:
            # n is a given fixed integer, evaluate into polynomial
            if n.is_negative:
                raise ValueError(
                    "The index n must be nonnegative integer (got %r)" % n)
            else:
                return laguerre_poly(n, x, alpha)

def nC(n, k=None, replacement=False):
    """Return the number of combinations of ``n`` items taken ``k`` at a time.

    Possible values for ``n``::
        integer - set of length ``n``
        sequence - converted to a multiset internally
        multiset - {element: multiplicity}

    If ``k`` is None then the total of all combinations of length 0
    through the number of items represented in ``n`` will be returned.

    If ``replacement`` is True then a given item can appear more than once
    in the ``k`` items. (For example, for 'ab' sets of 2 would include 'aa',
    'ab', and 'bb'.) The multiplicity of elements in ``n`` is ignored when
    ``replacement`` is True but the total number of elements is considered
    since no element can appear more times than the number of elements in
    ``n``.

    Examples
    ========

    >>> from sympy.functions.combinatorial.numbers import nC
    >>> from sympy.utilities.iterables import multiset_combinations
    >>> nC(3, 2)
    3
    >>> nC('abc', 2)
    3
    >>> nC('aab', 2)
    2

    When ``replacement`` is True, each item can have multiplicity
    equal to the length represented by ``n``:

    >>> nC('aabc', replacement=True)
    35
    >>> [len(list(multiset_combinations('aaaabbbbcccc', i))) for i in range(5)]
    [1, 3, 6, 10, 15]
    >>> sum(_)
    35

    If there are ``k`` items with multiplicities ``m_1, m_2, ..., m_k``
    then the total of all combinations of length 0 hrough ``k`` is the
    product, ``(m_1 + 1)*(m_2 + 1)*...*(m_k + 1)``. When the multiplicity
    of each item is 1 (i.e., k unique items) then there are 2**k
    combinations. For example, if there are 4 unique items, the total number
    of combinations is 16:

    >>> sum(nC(4, i) for i in range(5))
    16

    References
    ==========

    .. [1] http://en.wikipedia.org/wiki/Combination
    .. [2] http://tinyurl.com/cep849r

    See Also
    ========
    sympy.utilities.iterables.multiset_combinations
    """
    from sympy.functions.combinatorial.factorials import binomial
    from sympy.core.mul import prod

    if isinstance(n, SYMPY_INTS):
        if k is None:
            if not replacement:
                return 2**n
            return sum(nC(n, i, replacement) for i in range(n + 1))
        if k < 0:
            raise ValueError("k cannot be negative")
        if replacement:
            return binomial(n + k - 1, k)
        return binomial(n, k)
    if isinstance(n, _MultisetHistogram):
        N = n[_N]
        if k is None:
            if not replacement:
                return prod(m + 1 for m in n[_M])
            return sum(nC(n, i, replacement) for i in range(N + 1))
        elif replacement:
            return nC(n[_ITEMS], k, replacement)
        # assert k >= 0
        elif k in (1, N - 1):
            return n[_ITEMS]
        elif k in (0, N):
            return 1
        return _AOP_product(tuple(n[_M]))[k]
    else:
        return nC(_multiset_histogram(n), k, replacement)

 def _eval_rewrite_as_polynomial(self, n, x):
     from sympy import Sum
     k = Dummy("k")
     kern = (-1)**k*binomial(n, k)**2*((1 + x)/2)**(n - k)*((1 - x)/2)**k
     return Sum(kern, (k, 0, n))

 def _eval_rewrite_as_binomial(self, n):
     return binomial(2*n, n)/(n + 1)

def combsimp(expr):
    r"""
    Simplify combinatorial expressions.

    This function takes as input an expression containing factorials,
    binomials, Pochhammer symbol and other "combinatorial" functions,
    and tries to minimize the number of those functions and reduce
    the size of their arguments. The result is be given in terms of
    binomials and factorials.

    The algorithm works by rewriting all combinatorial functions as
    expressions involving rising factorials (Pochhammer symbols) and
    applies recurrence relations and other transformations applicable
    to rising factorials, to reduce their arguments, possibly letting
    the resulting rising factorial to cancel. Rising factorials with
    the second argument being an integer are expanded into polynomial
    forms and finally all other rising factorial are rewritten in terms
    more familiar functions. If the initial expression contained any
    combinatorial functions, the result is expressed using binomial
    coefficients and gamma functions. If the initial expression consisted
    of gamma functions alone, the result is expressed in terms of gamma
    functions.

    If the result is expressed using gamma functions, the following three
    additional steps are performed:

    1. Reduce the number of gammas by applying the reflection theorem
       gamma(x)*gamma(1-x) == pi/sin(pi*x).
    2. Reduce the number of gammas by applying the multiplication theorem
       gamma(x)*gamma(x+1/n)*...*gamma(x+(n-1)/n) == C*gamma(n*x).
    3. Reduce the number of prefactors by absorbing them into gammas, where
       possible.

    All transformation rules can be found (or was derived from) here:

    1. http://functions.wolfram.com/GammaBetaErf/Pochhammer/17/01/02/
    2. http://functions.wolfram.com/GammaBetaErf/Pochhammer/27/01/0005/

    Examples
    ========

    >>> from sympy.simplify import combsimp
    >>> from sympy import factorial, binomial
    >>> from sympy.abc import n, k

    >>> combsimp(factorial(n)/factorial(n - 3))
    n*(n - 2)*(n - 1)
    >>> combsimp(binomial(n+1, k+1)/binomial(n, k))
    (n + 1)/(k + 1)

    """

    # as a rule of thumb, if the expression contained gammas initially, it
    # probably makes sense to retain them
    as_gamma = not expr.has(factorial, binomial)


    expr = expr.replace(binomial,
        lambda n, k: _rf((n - k + 1).expand(), k.expand())/_rf(1, k.expand()))
    expr = expr.replace(factorial,
        lambda n: _rf(1, n.expand()))
    expr = expr.rewrite(gamma)
    expr = expr.replace(gamma,
        lambda n: _rf(1, (n - 1).expand()))

    if as_gamma:
        expr = expr.replace(_rf,
            lambda a, b: gamma(a + b)/gamma(a))
    else:
        expr = expr.replace(_rf,
            lambda a, b: binomial(a + b - 1, b)*factorial(b))

    def rule(n, k):
        coeff, rewrite = S.One, False

        cn, _n = n.as_coeff_Add()

        if _n and cn.is_Integer and cn:
            coeff *= _rf(_n + 1, cn)/_rf(_n - k + 1, cn)
            rewrite = True
            n = _n

        # this sort of binomial has already been removed by
        # rising factorials but is left here in case the order
        # of rule application is changed
        if k.is_Add:
            ck, _k = k.as_coeff_Add()
            if _k and ck.is_Integer and ck:
                coeff *= _rf(n - ck - _k + 1, ck)/_rf(_k + 1, ck)
                rewrite = True
                k = _k

        if rewrite:
            return coeff*binomial(n, k)

    expr = expr.replace(binomial, rule)

    def rule_gamma(expr, level=0):
        """ Simplify products of gamma functions further. """

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

#.........这里部分代码省略.........
                    y1.append(j * other)
                return HolonomicFunction(ann_self, self.x, self.x0, y1)

        ann_other = other.annihilator
        list_self = []
        list_other = []
        a = ann_self.order
        b = ann_other.order

        R = ann_self.parent.base
        K = R.get_field()

        for j in ann_self.listofpoly:
            list_self.append(K.new(j.rep))

        for j in ann_other.listofpoly:
            list_other.append(K.new(j.rep))

        # will be used to reduce the degree
        self_red = [-list_self[i] / list_self[a] for i in range(a)]

        other_red = [-list_other[i] / list_other[b] for i in range(b)]

        # coeff_mull[i][j] is the coefficient of Dx^i(f).Dx^j(g)
        coeff_mul = [[S(0) for i in range(b + 1)] for j in range(a + 1)]
        coeff_mul[0][0] = S(1)

        # making the ansatz
        lin_sys = [[coeff_mul[i][j] for i in range(a) for j in range(b)]]

        homo_sys = [[S(0) for q in range(a * b)]]
        homo_sys = NewMatrix(homo_sys).transpose()

        sol = (NewMatrix(lin_sys).transpose()).gauss_jordan_solve(homo_sys)

        # until a non trivial solution is found
        while sol[0].is_zero:

            # updating the coefficents Dx^i(f).Dx^j(g) for next degree
            for i in range(a - 1, -1, -1):
                for j in range(b - 1, -1, -1):
                    coeff_mul[i][j + 1] += coeff_mul[i][j]
                    coeff_mul[i + 1][j] += coeff_mul[i][j]
                    if isinstance(coeff_mul[i][j], K.dtype):
                        coeff_mul[i][j] = DMFdiff(coeff_mul[i][j])
                    else:
                        coeff_mul[i][j] = coeff_mul[i][j].diff()

            # reduce the terms to lower power using annihilators of f, g
            for i in range(a + 1):
                if not coeff_mul[i][b] == S(0):
                    for j in range(b):
                        coeff_mul[i][j] += other_red[j] * \
                            coeff_mul[i][b]
                    coeff_mul[i][b] = S(0)

            # not d2 + 1, as that is already covered in previous loop
            for j in range(b):
                if not coeff_mul[a][j] == 0:
                    for i in range(a):
                        coeff_mul[i][j] += self_red[i] * \
                            coeff_mul[a][j]
                    coeff_mul[a][j] = S(0)

            lin_sys.append([coeff_mul[i][j] for i in range(a)
                            for j in range(b)])

            sol = (NewMatrix(lin_sys).transpose()).gauss_jordan_solve(homo_sys)


        sol_ann = _normalize(sol[0][0:], self.annihilator.parent, negative=False)

        if self._have_init_cond and other._have_init_cond:
            if self.x0 == other.x0:

                # try to find more inital conditions
                y0_self = _extend_y0(self, sol_ann.order)
                y0_other = _extend_y0(other, sol_ann.order)
                # h(x0) = f(x0) * g(x0)
                y0 = [y0_self[0] * y0_other[0]]

                # coefficient of Dx^j(f)*Dx^i(g) in Dx^i(fg)
                for i in range(1, min(len(y0_self), len(y0_other))):
                    coeff = [[0 for i in range(i + 1)] for j in range(i + 1)]
                    for j in range(i + 1):
                        for k in range(i + 1):
                            if j + k == i:
                                coeff[j][k] = binomial(i, j)

                    sol = 0
                    for j in range(i + 1):
                        for k in range(i + 1):
                            sol += coeff[j][k]* y0_self[j] * y0_other[k]

                    y0.append(sol)
                return HolonomicFunction(sol_ann, self.x, self.x0, y0)
            else:
                raise NotImplementedError

        return HolonomicFunction(sol_ann, self.x)

 def fit(d, r, p):
     return -np.log(np.float(c.binomial(d, r)) * (p ** r) * (1 - p) ** (d - r))

def test_sympy__functions__combinatorial__factorials__binomial():
    from sympy.functions.combinatorial.factorials import binomial
    assert _test_args(binomial(2, x))

#.........这里部分代码省略.........
    order : int

    Examples
    ========

    >>> from sympy import log, atan, I
    >>> from sympy.series.formal import rational_algorithm as ra
    >>> from sympy.abc import x, k

    >>> ra(1 / (1 - x), x, k)
    (1, 0, 0)
    >>> ra(log(1 + x), x, k)
    (-(-1)**(-k)/k, 0, 1)

    >>> ra(atan(x), x, k, full=True)
    ((-I*(-I)**(-k)/2 + I*I**(-k)/2)/k, 0, 1)

    Notes
    =====

    By setting ``full=True``, range of admissible functions to be solved using
    ``rational_algorithm`` can be increased. This option should be used
    carefully as it can signifcantly slow down the computation as ``doit`` is
    performed on the :class:`RootSum` object returned by the ``apart`` function.
    Use ``full=False`` whenever possible.

    See Also
    ========

    sympy.polys.partfrac.apart

    References
    ==========

    .. [1] Formal Power Series - Dominik Gruntz, Wolfram Koepf
    .. [2] Power Series in Computer Algebra - Wolfram Koepf
    """
    from sympy.polys import RootSum, apart
    from sympy.integrals import integrate

    diff = f
    ds = []  # list of diff

    for i in range(order + 1):
        if i:
            diff = diff.diff(x)

        if diff.is_rational_function(x):
            coeff, sep = S.Zero, S.Zero

            terms = apart(diff, x, full=full)
            if terms.has(RootSum):
                terms = terms.doit()

            for t in Add.make_args(terms):
                num, den = t.as_numer_denom()
                if not den.has(x):
                    sep += t
                else:
                    if isinstance(den, Mul):
                        # m*(n*x - a)**j -> (n*x - a)**j
                        ind = den.as_independent(x)
                        den = ind[1]
                        num /= ind[0]

                    # (n*x - a)**j -> (x - b)
                    den, j = den.as_base_exp()
                    a, xterm = den.as_coeff_add(x)

                    # term -> m/x**n
                    if not a:
                        sep += t
                        continue

                    xc = xterm[0].coeff(x)
                    a /= -xc
                    num /= xc**j

                    ak = ((-1)**j * num *
                          binomial(j + k - 1, k).rewrite(factorial) /
                          a**(j + k))
                    coeff += ak

            # Hacky, better way?
            if coeff is S.Zero:
                return None
            if (coeff.has(x) or coeff.has(zoo) or coeff.has(oo) or
                    coeff.has(nan)):
                return None

            for j in range(i):
                coeff = (coeff / (k + j + 1))
                sep = integrate(sep, x)
                sep += (ds.pop() - sep).limit(x, 0)  # constant of integration
            return (coeff.subs(k, k - i), sep, i)

        else:
            ds.append(diff)

    return None

    def __mul__(self, other):

        ann_self = self.annihilator
        ann_other = other.annihilator
        list_self = ann_self.listofpoly
        list_other = ann_other.listofpoly
        a = ann_self.order
        b = ann_other.order

        for i, j in enumerate(list_self):
            if isinstance(j, ann_self.parent.base.dtype):
                list_self[i] = ann_self.parent.base.to_sympy(j)

        for i, j in enumerate(list_other):
            if isinstance(j, ann_other.parent.base.dtype):
                list_other[i] = ann_other.parent.base.to_sympy(j)
        # will be used to reduce the degree
        self_red = [-list_self[i] / list_self[a]
                    for i in range(a)]

        other_red = [-list_other[i] / list_other[b]
                     for i in range(b)]
        # coeff_mull[i][j] is the coefficient of Dx^i(f).Dx^j(g)
        coeff_mul = [[S(0) for i in range(b + 1)]
                     for j in range(a + 1)]
        coeff_mul[0][0] = S(1)
        # making the ansatz
        lin_sys = [[coeff_mul[i][j]
                    for i in range(a) for j in range(b)]]

        homo_sys = [[0 for q in range(a * b)]]
        homo_sys = Matrix(homo_sys).transpose()

        sol = (Matrix(lin_sys).transpose()).gauss_jordan_solve(homo_sys)
        # until a non trivial solution is found
        while sol[0].is_zero:
            # updating the coefficents Dx^i(f).Dx^j(g) for next degree
            for i in range(a - 1, -1, -1):
                for j in range(b - 1, -1, -1):
                    coeff_mul[i][j + 1] += coeff_mul[i][j]
                    coeff_mul[i + 1][j] += coeff_mul[i][j]
                    coeff_mul[i][j] = coeff_mul[i][j].diff()
            # reduce the terms to lower power using annihilators of f, g
            for i in range(a + 1):
                if not coeff_mul[i][b] == S(0):
                    for j in range(b):
                        coeff_mul[i][j] += other_red[j] * \
                            coeff_mul[i][b]
                    coeff_mul[i][b] = S(0)

            # not d2 + 1, as that is already covered in previous loop
            for j in range(b):
                if not coeff_mul[a][j] == 0:
                    for i in range(a):
                        coeff_mul[i][j] += self_red[i] * \
                            coeff_mul[a][j]
                    coeff_mul[a][j] = S(0)

            lin_sys.append([coeff_mul[i][j] for i in range(a)
                            for j in range(b)])

            sol = (Matrix(lin_sys).transpose()).gauss_jordan_solve(homo_sys)

        sol = sol[0] / sol[1]

        sol_ann = DifferentialOperator(_normalize(
            sol[0:], self.x, negative=False), ann_self.parent)

        if self._have_init_cond and other._have_init_cond:
            if self.x0 == other.x0:

                # try to find more inital conditions
                y0_self = _extend_y0(self, sol_ann.order)
                y0_other = _extend_y0(other, sol_ann.order)
                # h(x0) = f(x0) * g(x0)
                y0 = [y0_self[0] * y0_other[0]]

                # coefficient of Dx^j(f)*Dx^i(g) in Dx^i(fg)
                for i in range(1, min(len(y0_self), len(y0_other))):
                    coeff = [[0 for i in range(i + 1)] for j in range(i + 1)]
                    for j in range(i + 1):
                        for k in range(i + 1):
                            if j + k == i:
                                coeff[j][k] = binomial(i, j)

                    sol = 0
                    for j in range(i + 1):
                        for k in range(i + 1):
                            sol += coeff[j][k]* y0_self[j] * y0_other[k]

                    y0.append(sol)
                return HolonomicFunction(sol_ann, self.x, self.x0, y0)
            else:
                raise NotImplementedError

        return HolonomicFunction(sol_ann, self.x)