Python sympy.binomial函数代码示例

def test_limit_seq():
    assert limit(Sum(1/x, (x, 1, y)) - log(y), y, oo) == EulerGamma
    assert limit(Sum(1/x, (x, 1, y)) - 1/y, y, oo) == S.Infinity
    assert (limit(binomial(2*x, x) / Sum(binomial(2*y, y), (y, 1, x)), x, oo) ==
            S(3) / 4)
    assert (limit(Sum(y**2 * Sum(2**z/z, (z, 1, y)), (y, 1, x)) /
                  (2**x*x), x, oo) == 4)

def test_binomial_rewrite():
    n = Symbol("n", integer=True)
    k = Symbol("k", integer=True)

    assert binomial(n, k).rewrite(factorial) == factorial(n) / (factorial(k) * factorial(n - k))
    assert binomial(n, k).rewrite(gamma) == gamma(n + 1) / (gamma(k + 1) * gamma(n - k + 1))
    assert binomial(n, k).rewrite(ff) == ff(n, k) / factorial(k)

def test_gosper_sum():
    assert gosper_sum(1, (k, 0, n)) == 1 + n
    assert gosper_sum(k, (k, 0, n)) == n*(1 + n)/2
    assert gosper_sum(k**2, (k, 0, n)) == n*(1 + n)*(1 + 2*n)/6
    assert gosper_sum(k**3, (k, 0, n)) == n**2*(1 + n)**2/4

    assert gosper_sum(2**k, (k, 0, n)) == 2*2**n - 1

    assert gosper_sum(factorial(k), (k, 0, n)) is None
    assert gosper_sum(binomial(n, k), (k, 0, n)) is None

    assert gosper_sum(factorial(k)/k**2, (k, 0, n)) is None
    assert gosper_sum((k - 3)*factorial(k), (k, 0, n)) is None

    assert gosper_sum(k*factorial(k), k) == factorial(k)
    assert gosper_sum(
        k*factorial(k), (k, 0, n)) == n*factorial(n) + factorial(n) - 1

    assert gosper_sum((-1)**k*binomial(n, k), (k, 0, n)) == 0
    assert gosper_sum((
        -1)**k*binomial(n, k), (k, 0, m)) == -(-1)**m*(m - n)*binomial(n, m)/n

    assert gosper_sum((4*k + 1)*factorial(k)/factorial(2*k + 1), (k, 0, n)) == \
        (2*factorial(2*n + 1) - factorial(n))/factorial(2*n + 1)

    # issue 6033:
    assert gosper_sum(
        n*(n + a + b)*a**n*b**n/(factorial(n + a)*factorial(n + b)), \
        (n, 0, m)) == -a*b*(exp(m*log(a))*exp(m*log(b))*factorial(a)* \
        factorial(b) - factorial(a + m)*factorial(b + m))/(factorial(a)* \
        factorial(b)*factorial(a + m)*factorial(b + m))

def test_gosper_sum_AeqB_part2():
    f2a = n**2*a**n
    f2b = (n - r/2)*binomial(r, n)
    f2c = factorial(n - 1)**2/(factorial(n - x)*factorial(n + x))

    g2a = -a*(a + 1)/(a - 1)**3 + a**(
        m + 1)*(a**2*m**2 - 2*a*m**2 + m**2 - 2*a*m + 2*m + a + 1)/(a - 1)**3
    g2b = (m - r)*binomial(r, m)/2
    ff = factorial(1 - x)*factorial(1 + x)
    g2c = 1/ff*(
        1 - 1/x**2) + factorial(m)**2/(x**2*factorial(m - x)*factorial(m + x))

    g = gosper_sum(f2a, (n, 0, m))
    assert g is not None and simplify(g - g2a) == 0
    g = gosper_sum(f2b, (n, 0, m))
    assert g is not None and simplify(g - g2b) == 0
    g = gosper_sum(f2c, (n, 1, m))
    assert g is not None and simplify(g - g2c) == 0

    # delete these lines and unXFAIL the nan test below when it passes
    f2d = n*(n + a + b)*a**n*b**n/(factorial(n + a)*factorial(n + b))
    g2d = 1/(factorial(a - 1)*factorial(
        b - 1)) - a**(m + 1)*b**(m + 1)/(factorial(a + m)*factorial(b + m))
    assert simplify(
        sum(f2d.subs(n, i) for i in range(3)) - g2d.subs(m, 2)) == 0

def test_gosper_sum_AeqB_part1():
    f1a = n**4
    f1b = n**3*2**n
    f1c = 1/(n**2 + sqrt(5)*n - 1)
    f1d = n**4*4**n/binomial(2*n, n)
    f1e = factorial(3*n)/(factorial(n)*factorial(n + 1)*factorial(n + 2)*27**n)
    f1f = binomial(2*n, n)**2/((n + 1)*4**(2*n))
    f1g = (4*n - 1)*binomial(2*n, n)**2/((2*n - 1)**2*4**(2*n))
    f1h = n*factorial(n - S(1)/2)**2/factorial(n + 1)**2

    g1a = m*(m + 1)*(2*m + 1)*(3*m**2 + 3*m - 1)/30
    g1b = 26 + 2**(m + 1)*(m**3 - 3*m**2 + 9*m - 13)

    g = gosper_sum(f1a, (n, 0, m))
    assert g is not None and simplify(g - g1a) == 0
    g = gosper_sum(f1b, (n, 0, m))
    assert g is not None and simplify(g - g1b) == 0
    g = gosper_sum(f1c, (n, 0, m))
    assert g is not None # and simplify(g - g1c) == 0
    g = gosper_sum(f1d, (n, 0, m))
    assert g is not None # and simplify(g - g1d) == 0
    g = gosper_sum(f1e, (n, 0, m))
    assert g is not None # and simplify(g - g1e) == 0
    g = gosper_sum(f1f, (n, 0, m))
    assert g is not None # and simplify(g - g1f) == 0
    g = gosper_sum(f1g, (n, 0, m))
    assert g is not None # and simplify(g - g1g) == 0
    g = gosper_sum(f1h, (n, 0, m))
    assert g is not None # and simplify(g - g1h) == 0

def test_issue_9699():
    n, k = symbols('n k', real=True)
    x, y = symbols('x, y')
    assert combsimp((n + 1)*factorial(n)) == factorial(n + 1)
    assert combsimp((x + 1)*factorial(x)/gamma(y)) == gamma(x + 2)/gamma(y)
    assert combsimp(factorial(n)/n) == factorial(n - 1)
    assert combsimp(rf(x + n, k)*binomial(n, k)) == binomial(n, k)*gamma(k + n + x)/gamma(n + x)

def test_limitseq_sum():
    from sympy.abc import x, y, z
    assert limit_seq(Sum(1/x, (x, 1, y)) - log(y), y) == S.EulerGamma
    assert limit_seq(Sum(1/x, (x, 1, y)) - 1/y, y) == S.Infinity
    assert (limit_seq(binomial(2*x, x) / Sum(binomial(2*y, y), (y, 1, x)), x) ==
            S(3) / 4)
    assert (limit_seq(Sum(y**2 * Sum(2**z/z, (z, 1, y)), (y, 1, x)) /
                  (2**x*x), x) == 4)

def pxbar2( i, nup, bp ):
    bp2 = float(nup)/float(numX)
    if (nup==0):    
        if (i==0):
            return 1
        else:
            return 0
    return sp.binomial(Ns,i) * bp2**i* (1-bp2)**(Ns-i)* tdown(i/float(Ns)) / sum(sp.binomial(Ns,j) * bp2**j* (1-bp2)**(Ns-j)* tdown(j/float(Ns)) for j in xrange(0,Ns+1) )

 def dict(self):
     N, m, n = self.N, self.m, self.n
     N, m, n = list(map(sympify, (N, m, n)))
     density = dict((sympify(k),
                     Rational(binomial(m, k) * binomial(N - m, n - k),
                              binomial(N, n)))
                     for k in range(max(0, n + m - N), min(m, n) + 1))
     return density

def test_issue_2787():
    n, k = symbols('n k', positive=True, integer=True)
    p = symbols('p', positive=True)
    binomial_dist = binomial(n, k)*p**k*(1 - p)**(n - k)
    s = Sum(binomial_dist*k, (k, 0, n))
    res = s.doit().simplify()
    assert res == Piecewise((n*p, And(Or(-n + 1 < 0, -n + 1 >= 0),
        Or(-n + 1 < 0, Ne(p/(p - 1), 1)), p*Abs(1/(p - 1)) <= 1)),
        (Sum(k*p**k*(-p + 1)**(-k)*(-p + 1)**n*binomial(n, k), (k, 0, n)), True))

def test_issue_2787():
    n, k = symbols('n k', positive=True, integer=True)
    p = symbols('p', positive=True)
    binomial_dist = binomial(n, k)*p**k*(1 - p)**(n - k)
    s = Sum(binomial_dist*k, (k, 0, n))
    res = s.doit().simplify()
    assert res == Piecewise(
        (n*p, p/Abs(p - 1) <= 1),
        ((-p + 1)**n*Sum(k*p**k*(-p + 1)**(-k)*binomial(n, k), (k, 0, n)),
        True))

def test_gosper_sum_iterated():
    f1 = binomial(2*k, k)/4**k
    f2 = (1 + 2*n)*binomial(2*n, n)/4**n
    f3 = (1 + 2*n)*(3 + 2*n)*binomial(2*n, n)/(3*4**n)
    f4 = (1 + 2*n)*(3 + 2*n)*(5 + 2*n)*binomial(2*n, n)/(15*4**n)
    f5 = (1 + 2*n)*(3 + 2*n)*(5 + 2*n)*(7 + 2*n)*binomial(2*n, n)/(105*4**n)

    assert gosper_sum(f1, (k, 0, n)) == f2
    assert gosper_sum(f2, (n, 0, n)) == f3
    assert gosper_sum(f3, (n, 0, n)) == f4
    assert gosper_sum(f4, (n, 0, n)) == f5

def vasicek_base(N, k, p, rho):
    """

    :param N:
    :param k:
    :param p:
    :param rho:
    :return:
    """

    zmin = - settings.SCALE
    zmax = settings.SCALE
    grid = settings.GRID_POINTS

    dz = float(zmax - zmin) / float(grid - 1)
    a = stats.norm.ppf(p, loc=0.0, scale=1.0)
    integral = 0
    for i in range(1, grid):
        z = zmin + dz * i
        arg = (a - rho * z) / math.sqrt(1 - rho * rho)
        phi_den = stats.norm.pdf(z, loc=0.0, scale=1.0)
        phi_cum = stats.norm.cdf(arg, loc=0.0, scale=1.0)
        integrant = phi_den * math.pow(phi_cum, k) * math.pow(1 - phi_cum, N - k) * binomial(N, k)
        integral = integral + integrant
    return dz * integral

    def translate_from(self, src_expansion, src_coeff_exprs, src_rscale,
            dvec, tgt_rscale):
        if not isinstance(src_expansion, type(self)):
            raise RuntimeError("do not know how to translate %s to "
                    "Taylor multipole expansion"
                               % type(src_expansion).__name__)

        if not self.use_rscale:
            src_rscale = 1
            tgt_rscale = 1

        logger.info("building translation operator: %s(%d) -> %s(%d): start"
                % (type(src_expansion).__name__,
                    src_expansion.order,
                    type(self).__name__,
                    self.order))

        from sumpy.tools import mi_factorial

        src_mi_to_index = dict((mi, i) for i, mi in enumerate(
            src_expansion.get_coefficient_identifiers()))

        for i, mi in enumerate(src_expansion.get_coefficient_identifiers()):
            src_coeff_exprs[i] *= mi_factorial(mi)

        result = [0] * len(self.get_full_coefficient_identifiers())
        from pytools import generate_nonnegative_integer_tuples_below as gnitb

        for i, tgt_mi in enumerate(
                self.get_full_coefficient_identifiers()):

            tgt_mi_plus_one = tuple(mi_i + 1 for mi_i in tgt_mi)

            for src_mi in gnitb(tgt_mi_plus_one):
                try:
                    src_index = src_mi_to_index[src_mi]
                except KeyError:
                    # Omitted coefficients: not life-threatening
                    continue

                contrib = src_coeff_exprs[src_index]

                for idim in range(self.dim):
                    n = tgt_mi[idim]
                    k = src_mi[idim]
                    assert n >= k
                    from sympy import binomial
                    contrib *= (binomial(n, k)
                            * sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k))

                result[i] += (
                        contrib
                        * sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(src_mi))

            result[i] /= mi_factorial(tgt_mi)

        logger.info("building translation operator: done")
        return (
            self.derivative_wrangler.get_stored_mpole_coefficients_from_full(
                result, tgt_rscale))

    def __init__(self, states, interval, differential_order):
        """
        :param states: tuple of states in beginning and end of interval
        :param interval: time interval (tuple)
        :param differential_order: grade of differential flatness :math:`\\gamma`
        """
        self.yd = states
        self.t0 = interval[0]
        self.t1 = interval[1]
        self.dt = interval[1] - interval[0]
        gamma = differential_order  # + 1 # TODO check this against notes

        # setup symbolic expressions
        tau, k = sp.symbols('tau, k')

        alpha = sp.factorial(2 * gamma + 1)

        f = sp.binomial(gamma, k) * (-1) ** k * tau ** (gamma + k + 1) / (gamma + k + 1)
        phi = alpha / sp.factorial(gamma) ** 2 * sp.summation(f, (k, 0, gamma))

        # differentiate phi(tau), index in list corresponds to order
        dphi_sym = [phi]  # init with phi(tau)
        for order in range(differential_order):
            dphi_sym.append(dphi_sym[-1].diff(tau))

        # lambdify
        self.dphi_num = []
        for der in dphi_sym:
            self.dphi_num.append(sp.lambdify(tau, der, 'numpy'))