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_simplify_other():
    assert simplify(sin(x)**2 + cos(x)**2) == 1
    assert simplify(gamma(x + 1)/gamma(x)) == x
    assert simplify(sin(x)**2 + cos(x)**2 + factorial(x)/gamma(x)) == 1 + x
    assert simplify(Eq(sin(x)**2 + cos(x)**2, factorial(x)/gamma(x))) == Eq(1, x)
    nc = symbols('nc', commutative=False)
    assert simplify(x + x*nc) == x*(1 + nc)

def test_simplify_other():
    assert simplify(sin(x) ** 2 + cos(x) ** 2) == 1
    assert simplify(gamma(x + 1) / gamma(x)) == x
    assert simplify(sin(x) ** 2 + cos(x) ** 2 + factorial(x) / gamma(x)) == 1 + x
    assert simplify(Eq(sin(x) ** 2 + cos(x) ** 2, factorial(x) / gamma(x))) == Eq(x, 1)
    nc = symbols("nc", commutative=False)
    assert simplify(x + x * nc) == x * (1 + nc)
    # issue 6123
    # f = exp(-I*(k*sqrt(t) + x/(2*sqrt(t)))**2)
    # ans = integrate(f, (k, -oo, oo), conds='none')
    ans = I * (
        -pi
        * x
        * exp(-3 * I * pi / 4 + I * x ** 2 / (4 * t))
        * erf(x * exp(-3 * I * pi / 4) / (2 * sqrt(t)))
        / (2 * sqrt(t))
        + pi * x * exp(-3 * I * pi / 4 + I * x ** 2 / (4 * t)) / (2 * sqrt(t))
    ) * exp(-I * x ** 2 / (4 * t)) / (sqrt(pi) * x) - I * sqrt(pi) * (
        -erf(x * exp(I * pi / 4) / (2 * sqrt(t))) + 1
    ) * exp(
        I * pi / 4
    ) / (
        2 * sqrt(t)
    )
    assert simplify(ans) == -(-1) ** (S(3) / 4) * sqrt(pi) / sqrt(t)
    # issue 6370
    assert simplify(2 ** (2 + x) / 4) == 2 ** x

def test_evalf_integer_parts():
    a = floor(log(8)/log(2) - exp(-1000), evaluate=False)
    b = floor(log(8)/log(2), evaluate=False)
    assert a.evalf() == 3
    assert b.evalf() == 3
    # equals, as a fallback, can still fail but it might succeed as here
    assert ceiling(10*(sin(1)**2 + cos(1)**2)) == 10

    assert int(floor(factorial(50)/E, evaluate=False).evalf(70)) == \
        long(11188719610782480504630258070757734324011354208865721592720336800)
    assert int(ceiling(factorial(50)/E, evaluate=False).evalf(70)) == \
        long(11188719610782480504630258070757734324011354208865721592720336801)
    assert int(floor((GoldenRatio**999 / sqrt(5) + Rational(1, 2)))
               .evalf(1000)) == fibonacci(999)
    assert int(floor((GoldenRatio**1000 / sqrt(5) + Rational(1, 2)))
               .evalf(1000)) == fibonacci(1000)

    assert ceiling(x).evalf(subs={x: 3}) == 3
    assert ceiling(x).evalf(subs={x: 3*I}) == 3*I
    assert ceiling(x).evalf(subs={x: 2 + 3*I}) == 2 + 3*I
    assert ceiling(x).evalf(subs={x: 3.}) == 3
    assert ceiling(x).evalf(subs={x: 3.*I}) == 3*I
    assert ceiling(x).evalf(subs={x: 2. + 3*I}) == 2 + 3*I

    assert float((floor(1.5, evaluate=False)+1/9).evalf()) == 1 + 1/9
    assert float((floor(0.5, evaluate=False)+20).evalf()) == 20

def test_jacobi():
    n = Symbol("n")
    a = Symbol("a")
    b = Symbol("b")

    assert jacobi(0, a, b, x) == 1
    assert jacobi(1, a, b, x) == a/2 - b/2 + x*(a/2 + b/2 + 1)

    assert jacobi(n, a, a, x) == RisingFactorial(a + 1, n)*gegenbauer(n, a + S(1)/2, x)/RisingFactorial(2*a + 1, n)
    assert jacobi(n, a, -a, x) == ((-1)**a*(-x + 1)**(-a/2)*(x + 1)**(a/2)*assoc_legendre(n, a, x)*
                                   factorial(-a + n)*gamma(a + n + 1)/(factorial(a + n)*gamma(n + 1)))
    assert jacobi(n, -b, b, x) == ((-x + 1)**(b/2)*(x + 1)**(-b/2)*assoc_legendre(n, b, x)*
                                   gamma(-b + n + 1)/gamma(n + 1))
    assert jacobi(n, 0, 0, x) == legendre(n, x)
    assert jacobi(n, S.Half, S.Half, x) == RisingFactorial(S(3)/2, n)*chebyshevu(n, x)/factorial(n + 1)
    assert jacobi(n, -S.Half, -S.Half, x) == RisingFactorial(S(1)/2, n)*chebyshevt(n, x)/factorial(n)

    X = jacobi(n, a, b, x)
    assert isinstance(X, jacobi)

    assert jacobi(n, a, b, -x) == (-1)**n*jacobi(n, b, a, x)
    assert jacobi(n, a, b, 0) == 2**(-n)*gamma(a + n + 1)*hyper((-b - n, -n), (a + 1,), -1)/(factorial(n)*gamma(a + 1))
    assert jacobi(n, a, b, 1) == RisingFactorial(a + 1, n)/factorial(n)

    m = Symbol("m", positive=True)
    assert jacobi(m, a, b, oo) == oo*RisingFactorial(a + b + m + 1, m)

    assert conjugate(jacobi(m, a, b, x)) == jacobi(m, conjugate(a), conjugate(b), conjugate(x))

    assert diff(jacobi(n,a,b,x), n) == Derivative(jacobi(n, a, b, x), n)
    assert diff(jacobi(n,a,b,x), x) == (a/2 + b/2 + n/2 + S(1)/2)*jacobi(n - 1, a + 1, b + 1, x)

def test_factorials():
    n = Symbol('n', integer=True)
    assert factorial(-2) == 0
    assert factorial(0) == 1
    assert factorial(7) == 5040
    assert factorial(n).func == Factorial
    assert factorial(2*n).func == Factorial

    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 cases(n):         
    """
    Procedure Name: cases
    Purpose: Generates all possible arrival/departure sequences for
                n customers in an M/M/1 queue initially empty and
                idle.
    Arguments:  1. n: the total number of customers in the system
    Output:     1. C: a list of sequences consisting of 1s and -1s,
                        where 1s represent an arrival and -1s
                        represent a departure
    """
    # Compute the nth Catalan number
    c=factorial(2*n)/factorial(n)/factorial(n+1)
    C=np.zeros((c,2*n))
    for i in range(c):
        # Initialize the matrix C
        if i==0:
            C[i]=ini(n)
        # Produce the successor the C[i]
        else:
            C[i]=swapa(n,C[i-1])
        # Check to see if the successor is legal
        #   If not, call swapb
        if okay(n,C[i])==False:
            C[i]=swapb(n,C[i-1])
    return C

def test_factorials():
    n = Symbol('n', integer=True)
    assert factorial(-2) == 0
    assert factorial(0) == 1
    assert factorial(7) == 5040
    assert isinstance(factorial(n), Factorial)
    assert isinstance(factorial(2*n), Factorial)

    def Theta_Equation(self, l="undef", m="undef", theta=theta):
        """ l - orbital quantum number
            m - magnetic quantum number """

        if l == "undef":
            l = self.l_val
        if m == "undef":
            m = self.m_val

        # Prevents integer division and float l,m
        l = int(l)
        m = int(m)

        mode = self.select_exec_mode(theta)

        gL = self.Generalized_Legendre(l, m, theta).subs(theta, sympy.cos(theta))

        if gL:
            rat_part = sympy.sqrt(
                Rational((2 * l + 1), 2) * Rational(factorial(l - sympy.abs(m)), factorial(l + sympy.abs(m)))
            )
            if mode == "numer":
                return (rat_part * gL).evalf()
            elif mode == "analit":
                return (rat_part * gL).expand()

        return False

def test_latex_functions():
    assert latex(exp(x)) == "$e^{x}$"
    assert latex(exp(1) + exp(2)) == "$e + e^{2}$"

    f = Function("f")
    assert latex(f(x)) == "$\\operatorname{f}\\left(x\\right)$"

    beta = Function("beta")

    assert latex(beta(x)) == r"$\operatorname{beta}\left(x\right)$"
    assert latex(sin(x)) == r"$\operatorname{sin}\left(x\right)$"
    assert latex(sin(x), fold_func_brackets=True) == r"$\operatorname{sin}x$"
    assert latex(sin(2 * x ** 2), fold_func_brackets=True) == r"$\operatorname{sin}2 x^{2}$"
    assert latex(sin(x ** 2), fold_func_brackets=True) == r"$\operatorname{sin}x^{2}$"

    assert latex(asin(x) ** 2) == r"$\operatorname{asin}^{2}\left(x\right)$"
    assert latex(asin(x) ** 2, inv_trig_style="full") == r"$\operatorname{arcsin}^{2}\left(x\right)$"
    assert latex(asin(x) ** 2, inv_trig_style="power") == r"$\operatorname{sin}^{-1}\left(x\right)^{2}$"
    assert latex(asin(x ** 2), inv_trig_style="power", fold_func_brackets=True) == r"$\operatorname{sin}^{-1}x^{2}$"

    assert latex(factorial(k)) == r"$k!$"
    assert latex(factorial(-k)) == r"$\left(- k\right)!$"

    assert latex(floor(x)) == r"$\lfloor{x}\rfloor$"
    assert latex(ceiling(x)) == r"$\lceil{x}\rceil$"
    assert latex(abs(x)) == r"$\lvert{x}\rvert$"
    assert latex(re(x)) == r"$\Re{x}$"
    assert latex(im(x)) == r"$\Im{x}$"
    assert latex(conjugate(x)) == r"$\overline{x}$"
    assert latex(gamma(x)) == r"$\operatorname{\Gamma}\left(x\right)$"
    assert latex(Order(x)) == r"$\operatorname{\mathcal{O}}\left(x\right)$"

def test_rsolve():
    f = y(n+2) - y(n+1) - y(n)
    h = sqrt(5)*(S.Half + S.Half*sqrt(5))**n \
      - sqrt(5)*(S.Half - S.Half*sqrt(5))**n

    assert rsolve(f, y(n)) in [
        C0*(S.Half - S.Half*sqrt(5))**n + C1*(S.Half + S.Half*sqrt(5))**n,
        C1*(S.Half - S.Half*sqrt(5))**n + C0*(S.Half + S.Half*sqrt(5))**n,
    ]

    assert rsolve(f, y(n), [      0,      5 ]) == h
    assert rsolve(f, y(n), {   0 :0,   1 :5 }) == h
    assert rsolve(f, y(n), { y(0):0, y(1):5 }) == h
    assert rsolve(y(n) - y(n-1) - y(n-2), y(n), [0, 5]) == h

    f = (n-1)*y(n+2) - (n**2+3*n-2)*y(n+1) + 2*n*(n+1)*y(n)
    g = C1*factorial(n) + C0*2**n
    h = -3*factorial(n) + 3*2**n

    assert rsolve(f, y(n)) == g
    assert rsolve(f, y(n), []) == g
    assert rsolve(f, y(n), {}) == g

    assert rsolve(f, y(n), [      0,      3 ]) == h
    assert rsolve(f, y(n), {   0 :0,   1 :3 }) == h
    assert rsolve(f, y(n), { y(0):0, y(1):3 }) == h

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)
    g1c = (m + 1)*(m*(m**2 - 7*m + 3)*sqrt(5) - (3*m**3 - 7*m**2 + 19*m - 6))/(2*m**3*sqrt(5) + m**4 + m**2 - 1)/6
    g1d = -S(2)/231 + 24**m*(m + 1)*(63*m**4 + 112*m**3 + 18*m**2 - 22*m + 3)/(693*binomial(2*m, m))
    g1e = -S(9)/2 + (81*m**2 + 261*m + 200)*factorial(3*m + 2)/(40*27**m*factorial(m)*factorial(m + 1)*factorial(m + 2))
    g1f = (2*m + 1)**2*binomial(2*m, m)**2/(4**(2*m)*(m + 1))
    g1g = -binomial(2*m, m)**2/4**(2*m)
    g1h = -(2*m + 1)**2*(3*m + 4)*factorial(m - S(1)/2)**2/factorial(m + 1)**2

    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 create_mgf_expr(self):
     """ Create sympy expression for the simple taylor expansion of the mgf. """
     indivs = range(self.n_indivs)
     branches = []
     for branch_size in range(1, self.n_indivs):
         branches += itertools.combinations(indivs, branch_size)
     mgf_base = sympy.Integer(1)
     num_loci = sympy.symbols("L", integer=True)
     for mom_order in range(1, self.moment_order + 1):
         m_max = self.moment_order - mom_order + 1
         branch_combos = itertools.combinations_with_replacement(branches, mom_order)
         for branch_combo in branch_combos:
             branch_mult = list(Counter(branch_combo).values())
             term_factor = sympy.Integer(1)
             for multi in branch_mult:
                 term_factor *= sympy.factorial(multi)
             bc_term = (make_subscr(branch_combo)*
                        sympy.symbols('theta')**mom_order*term_factor**-1)
             for branch in branch_combo:
                 branch_sum = sympy.Integer(0)
                 for indiv in branch:
                     branch_sum += sympy.symbols('k_' + str(indiv))
                 branch_term = sympy.Integer(0)
                 for mut_order in range(1, m_max + 1):
                     branch_term += (sympy.symbols('m_' + str(mut_order))/
                                     sympy.factorial(mut_order)*branch_sum**mut_order)
                 bc_term *= branch_term
             mgf_base += bc_term
     return mgf_base**num_loci

 def d(cls, j, m, mp, beta):
     """Wigner's lowercase d function."""
     # TODO: This does not do a good job of simplifying the trig functions.
     # The Jacobi Polynomial expansion is probably best as it is a simple
     # sum. But, this version does give correct answers and uses
     # Eq. 7 in Section 4.3.2 of Varshalovich.
     j = sympify(j)
     cosbeta = Symbol('cosbeta')
     x = 1-cosbeta
     y = 1+cosbeta
     jmmp = j-mp
     jpm = j+m
     jpmp = j+mp
     jmm = j-m
     mpmm2 = (mp-m)//2
     mpmp2 = (m+mp)//2
     result = (-1)**jmmp
     result *= 2**(-j)
     result *= sqrt(factorial(jpm)/(factorial(jmm)*factorial(jpmp)*factorial(jmmp)))
     result *= (x**mpmm2)*(y**(-mpmp2))
     result *= diff((x**jmmp)*(y**jpmp), cosbeta, jmm)
     if (2*j).is_Integer and not j.is_Integer:
         result = result.subs(1+cosbeta, 2*cos(beta/2)**2)
         result = result.subs(1-cosbeta, 2*sin(beta/2)**2)
     else:
         result = result.subs(cosbeta, cos(beta))
     return result

def test_latex_functions():
    assert latex(exp(x)) == "$e^{x}$"
    assert latex(exp(1)+exp(2)) == "$e + e^{2}$"

    f = Function('f')
    assert latex(f(x)) == '$\\operatorname{f}\\left(x\\right)$'

    beta = Function('beta')

    assert latex(beta(x)) == r"$\operatorname{beta}\left(x\right)$"
    assert latex(sin(x)) == r"$\operatorname{sin}\left(x\right)$"
    assert latex(sin(x), fold_func_brackets=True) == r"$\operatorname{sin}x$"
    assert latex(sin(2*x**2), fold_func_brackets=True) == \
    r"$\operatorname{sin}2 x^{2}$"

    assert latex(factorial(k)) == r"$k!$"
    assert latex(factorial(-k)) == r"$\left(- k\right)!$"

    assert latex(floor(x)) == r"$\lfloor{x}\rfloor$"
    assert latex(ceiling(x)) == r"$\lceil{x}\rceil$"
    assert latex(abs(x)) == r"$\lvert{x}\rvert$"
    assert latex(re(x)) == r"$\Re{x}$"
    assert latex(im(x)) == r"$\Im{x}$"
    assert latex(conjugate(x)) == r"$\overline{x}$"
    assert latex(gamma(x)) == r"$\operatorname{\Gamma}\left(x\right)$"
    assert latex(Order(x)) == r"$\operatorname{\mathcal{O}}\left(x\right)$"

def test_equals():
    assert (-3 - sqrt(5) + (-sqrt(10) / 2 - sqrt(2) / 2) ** 2).equals(0)
    assert (x ** 2 - 1).equals((x + 1) * (x - 1))
    assert (cos(x) ** 2 + sin(x) ** 2).equals(1)
    assert (a * cos(x) ** 2 + a * sin(x) ** 2).equals(a)
    r = sqrt(2)
    assert (-1 / (r + r * x) + 1 / r / (1 + x)).equals(0)
    assert factorial(x + 1).equals((x + 1) * factorial(x))
    assert sqrt(3).equals(2 * sqrt(3)) is False
    assert (sqrt(5) * sqrt(3)).equals(sqrt(3)) is False
    assert (sqrt(5) + sqrt(3)).equals(0) is False
    assert (sqrt(5) + pi).equals(0) is False
    assert meter.equals(0) is False
    assert (3 * meter ** 2).equals(0) is False

    # from integrate(x*sqrt(1+2*x), x);
    # diff is zero, but differentiation does not show it
    i = 2 * sqrt(2) * x ** (S(5) / 2) * (1 + 1 / (2 * x)) ** (S(5) / 2) / 5 + 2 * sqrt(2) * x ** (S(3) / 2) * (
        1 + 1 / (2 * x)
    ) ** (S(5) / 2) / (-6 - 3 / x)
    ans = sqrt(2 * x + 1) * (6 * x ** 2 + x - 1) / 15
    diff = i - ans
    assert diff.equals(0) is not False  # should be True, but now it's None
    # XXX TODO add a force=True option to equals to posify both
    # self and other before beginning comparisions
    p = Symbol("p", positive=True)
    assert diff.subs(x, p).equals(0) is True

def test_simple_products():
    assert Product(S.NaN, (x, 1, 3)) is S.NaN
    assert product(S.NaN, (x, 1, 3)) is S.NaN
    assert Product(x, (n, a, a)).doit() == x
    assert Product(x, (x, a, a)).doit() == a
    assert Product(x, (y, 1, a)).doit() == x**a
    lo, hi = 1, 2
    s1 = Product(n, (n, lo, hi))
    s2 = Product(n, (n, hi, lo))
    assert s1 != s2
    assert s1.doit() == s2.doit() == 2
    lo, hi = x, x + 1
    s1 = Product(n, (n, lo, hi))
    s2 = Product(n, (n, hi, lo))
    assert s1 != s2
    assert s1.doit() == s2.doit() == x*(x + 1)
    assert Product(Integral(2*x, (x, 1, y)) + 2*x, (x, 1, 2)).doit() == \
        (y**2 + 1)*(y**2 + 3)
    assert product(2, (n, a, b)) == 2**(b - a + 1)
    assert product(n, (n, 1, b)) == factorial(b)
    assert product(n**3, (n, 1, b)) == factorial(b)**3
    assert product(3**(2 + n), (n, a, b)) \
        == 3**(2*(1 - a + b) + b/2 + (b**2)/2 + a/2 - (a**2)/2)
    assert product(cos(n), (n, 3, 5)) == cos(3)*cos(4)*cos(5)
    assert product(cos(n), (n, x, x + 2)) == cos(x)*cos(x + 1)*cos(x + 2)
    assert isinstance(product(cos(n), (n, x, x + S.Half)), Product)
    # If Product managed to evaluate this one, it most likely got it wrong!
    assert isinstance(Product(n**n, (n, 1, b)), Product)

    def symmetry(self):
        r"""
        The symmetry $\\sigma(t)$ of a rooted tree is...

        **Examples**::

            >>> from nodepy import rooted_trees as rt
            >>> tree=rt.RootedTree('{T^2{T{T}}}')
            >>> tree.symmetry()
            2

        **Reference**:

            - [butcher2003]_ p. 127, eq. 301(b)
        """
        from nodepy.strmanip import getint
        if self=='T': return 1
        sigma=1
        if self[1]=='T':
            try: sigma=factorial(getint(self[3:]))
            except: pass
        nleaves,subtrees=self._parse_subtrees()
        while len(subtrees)>0:
            st=subtrees[0]
            nst=subtrees.count(st)
            sigma*=factorial(nst)*st.symmetry()**nst
            while st in subtrees: subtrees.remove(st)
        return sigma

def test_factorial_diff():
    n = Symbol('n', integer=True)

    assert factorial(n).diff(n) == \
        gamma(1 + n)*polygamma(0, 1 + n)
    assert factorial(n**2).diff(n) == \
        2*n*gamma(1 + n**2)*polygamma(0, 1 + n**2)