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'))

def test_evalf_fast_series():
    # Euler transformed series for sqrt(1+x)
    assert NS(Sum(
        fac(2*n + 1)/fac(n)**2/2**(3*n + 1), (n, 0, oo)), 100) == NS(sqrt(2), 100)

    # Some series for exp(1)
    estr = NS(E, 100)
    assert NS(Sum(1/fac(n), (n, 0, oo)), 100) == estr
    assert NS(1/Sum((1 - 2*n)/fac(2*n), (n, 0, oo)), 100) == estr
    assert NS(Sum((2*n + 1)/fac(2*n), (n, 0, oo)), 100) == estr
    assert NS(Sum((4*n + 3)/2**(2*n + 1)/fac(2*n + 1), (n, 0, oo))**2, 100) == estr

    pistr = NS(pi, 100)
    # Ramanujan series for pi
    assert NS(9801/sqrt(8)/Sum(fac(
        4*n)*(1103 + 26390*n)/fac(n)**4/396**(4*n), (n, 0, oo)), 100) == pistr
    assert NS(1/Sum(
        binomial(2*n, n)**3 * (42*n + 5)/2**(12*n + 4), (n, 0, oo)), 100) == pistr
    # Machin's formula for pi
    assert NS(16*Sum((-1)**n/(2*n + 1)/5**(2*n + 1), (n, 0, oo)) -
        4*Sum((-1)**n/(2*n + 1)/239**(2*n + 1), (n, 0, oo)), 100) == pistr

    # Apery's constant
    astr = NS(zeta(3), 100)
    P = 126392*n**5 + 412708*n**4 + 531578*n**3 + 336367*n**2 + 104000* \
        n + 12463
    assert NS(Sum((-1)**n * P / 24 * (fac(2*n + 1)*fac(2*n)*fac(
        n))**3 / fac(3*n + 2) / fac(4*n + 3)**3, (n, 0, oo)), 100) == astr
    assert NS(Sum((-1)**n * (205*n**2 + 250*n + 77)/64 * fac(n)**10 /
              fac(2*n + 1)**5, (n, 0, oo)), 100) == astr

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

    X = UniformSum('x', n)
    assert density(X)(x) == (Sum((-1)**_k*(-_k + x)**(n - 1)
                        *binomial(n, _k), (_k, 0, floor(x)))/factorial(n - 1))

def _estimate_gradients_2d_global():

    #
    # Compute
    #
    #

    f1, f2, df1, df2, x = symbols(['f1', 'f2', 'df1', 'df2', 'x'])
    c = [f1, (df1 + 3*f1)/3, (df2 + 3*f2)/3, f2]

    w = 0
    for k in range(4):
        w += binomial(3, k) * c[k] * x**k*(1-x)**(3-k)

    wpp = w.diff(x, 2).expand()
    intwpp2 = (wpp**2).integrate((x, 0, 1)).expand()

    A = Matrix([[intwpp2.coeff(df1**2), intwpp2.coeff(df1*df2)/2],
                [intwpp2.coeff(df1*df2)/2, intwpp2.coeff(df2**2)]])

    B = Matrix([[intwpp2.coeff(df1).subs(df2, 0)],
                [intwpp2.coeff(df2).subs(df1, 0)]]) / 2

    print("A")
    print(A)
    print("B")
    print(B)
    print("solution")
    print(A.inv() * B)

def test_assoc_laguerre():
    n = Symbol("n")
    m = Symbol("m")
    alpha = Symbol("alpha")

    # generalized Laguerre polynomials:
    assert assoc_laguerre(0, alpha, x) == 1
    assert assoc_laguerre(1, alpha, x) == -x + alpha + 1
    assert assoc_laguerre(2, alpha, x).expand() == \
        (x**2/2 - (alpha + 2)*x + (alpha + 2)*(alpha + 1)/2).expand()
    assert assoc_laguerre(3, alpha, x).expand() == \
        (-x**3/6 + (alpha + 3)*x**2/2 - (alpha + 2)*(alpha + 3)*x/2 +
        (alpha + 1)*(alpha + 2)*(alpha + 3)/6).expand()

    # Test the lowest 10 polynomials with laguerre_poly, to make sure it works:
    for i in range(10):
        assert assoc_laguerre(i, 0, x).expand() == laguerre_poly(i, x)

    X = assoc_laguerre(n, m, x)
    assert isinstance(X, assoc_laguerre)

    assert assoc_laguerre(n, 0, x) == laguerre(n, x)
    assert assoc_laguerre(n, alpha, 0) == binomial(alpha + n, alpha)

    assert diff(assoc_laguerre(n, alpha, x), x) == \
        -assoc_laguerre(n - 1, alpha + 1, x)

    assert conjugate(assoc_laguerre(n, alpha, x)) == \
        assoc_laguerre(n, conjugate(alpha), conjugate(x))

    raises(ValueError, lambda: assoc_laguerre(-2.1, alpha, x))

def test_rf_eval_apply():
    x, y = symbols('x,y')
    n, k = symbols('n k', integer=True)
    m = Symbol('m', integer=True, nonnegative=True)

    assert rf(nan, y) == nan
    assert rf(x, nan) == nan

    assert rf(x, y) == rf(x, y)

    assert rf(oo, 0) == 1
    assert rf(-oo, 0) == 1

    assert rf(oo, 6) == oo
    assert rf(-oo, 7) == -oo

    assert rf(oo, -6) == oo
    assert rf(-oo, -7) == oo

    assert rf(x, 0) == 1
    assert rf(x, 1) == x
    assert rf(x, 2) == x*(x + 1)
    assert rf(x, 3) == x*(x + 1)*(x + 2)
    assert rf(x, 5) == x*(x + 1)*(x + 2)*(x + 3)*(x + 4)

    assert rf(x, -1) == 1/(x - 1)
    assert rf(x, -2) == 1/((x - 1)*(x - 2))
    assert rf(x, -3) == 1/((x - 1)*(x - 2)*(x - 3))

    assert rf(1, 100) == factorial(100)

    assert rf(x**2 + 3*x, 2) == (x**2 + 3*x)*(x**2 + 3*x + 1)
    assert isinstance(rf(x**2 + 3*x, 2), Mul)
    assert rf(x**3 + x, -2) == 1/((x**3 + x - 1)*(x**3 + x - 2))

    assert rf(Poly(x**2 + 3*x, x), 2) == Poly(x**4 + 8*x**3 + 19*x**2 + 12*x, x)
    assert isinstance(rf(Poly(x**2 + 3*x, x), 2), Poly)
    raises(ValueError, lambda: rf(Poly(x**2 + 3*x, x, y), 2))
    assert rf(Poly(x**3 + x, x), -2) == 1/(x**6 - 9*x**5 + 35*x**4 - 75*x**3 + 94*x**2 - 66*x + 20)
    raises(ValueError, lambda: rf(Poly(x**3 + x, x, y), -2))

    assert rf(x, m).is_integer is None
    assert rf(n, k).is_integer is None
    assert rf(n, m).is_integer is True
    assert rf(n, k + pi).is_integer is False
    assert rf(n, m + pi).is_integer is False
    assert rf(pi, m).is_integer is False

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

    import random
    from mpmath import rf as mpmath_rf
    for i in range(100):
        x = -500 + 500 * random.random()
        k = -500 + 500 * random.random()
        assert (abs(mpmath_rf(x, k) - rf(x, k)) < 10**(-15))