def test_sympy_parser():
    x = Symbol('x')
    inputs = {
        '2*x': 2 * x,
        '3.00': Float(3),
        '22/7': Rational(22, 7),
        '2+3j': 2 + 3*I,
        'exp(x)': exp(x),
        'x!': factorial(x),
        'x!!': factorial2(x),
        '(x + 1)! - 1': factorial(x + 1) - 1,
        '3.[3]': Rational(10, 3),
        '.0[3]': Rational(1, 30),
        '3.2[3]': Rational(97, 30),
        '1.3[12]': Rational(433, 330),
        '1 + 3.[3]': Rational(13, 3),
        '1 + .0[3]': Rational(31, 30),
        '1 + 3.2[3]': Rational(127, 30),
        '.[0011]': Rational(1, 909),
        '0.1[00102] + 1': Rational(366697, 333330),
        '1.[0191]': Rational(10190, 9999),
        '10!': 3628800,
        '-(2)': -Integer(2),
        '[-1, -2, 3]': [Integer(-1), Integer(-2), Integer(3)],
        'Symbol("x").free_symbols': x.free_symbols,
        "S('S(3).n(n=3)')": 3.00,
        'factorint(12, visual=True)': Mul(
            Pow(2, 2, evaluate=False),
            Pow(3, 1, evaluate=False),
            evaluate=False),
        'Limit(sin(x), x, 0, dir="-")': Limit(sin(x), x, 0, dir='-'),

    }
    for text, result in inputs.items():
        assert parse_expr(text) == result

def test_catalan():
    n = Symbol('n', integer=True)
    m = Symbol('n', integer=True, positive=True)

    catalans = [1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786]
    for i, c in enumerate(catalans):
        assert catalan(i) == c
        assert catalan(n).rewrite(factorial).subs(n, i) == c
        assert catalan(n).rewrite(Product).subs(n, i).doit() == c

    assert catalan(x) == catalan(x)
    assert catalan(2*x).rewrite(binomial) == binomial(4*x, 2*x)/(2*x + 1)
    assert catalan(Rational(1, 2)).rewrite(gamma) == 8/(3*pi)
    assert catalan(Rational(1, 2)).rewrite(factorial).rewrite(gamma) ==\
        8 / (3 * pi)
    assert catalan(3*x).rewrite(gamma) == 4**(
        3*x)*gamma(3*x + Rational(1, 2))/(sqrt(pi)*gamma(3*x + 2))
    assert catalan(x).rewrite(hyper) == hyper((-x + 1, -x), (2,), 1)

    assert catalan(n).rewrite(factorial) == factorial(2*n) / (factorial(n + 1)
                                                              * factorial(n))
    assert isinstance(catalan(n).rewrite(Product), catalan)
    assert isinstance(catalan(m).rewrite(Product), Product)

    assert diff(catalan(x), x) == (polygamma(
        0, x + Rational(1, 2)) - polygamma(0, x + 2) + log(4))*catalan(x)

    assert catalan(x).evalf() == catalan(x)
    c = catalan(S.Half).evalf()
    assert str(c) == '0.848826363156775'
    c = catalan(I).evalf(3)
    assert str((re(c), im(c))) == '(0.398, -0.0209)'

def psi_n(n, x, m, omega):
    """
    Returns the wavefunction psi_{n} for the One-dimensional harmonic oscillator.

    ``n``
        the "nodal" quantum number.  Corresponds to the number of nodes in the
        wavefunction.  n >= 0
    ``x``
        x coordinate
    ``m``
        mass of the particle
    ``omega``
        angular frequency of the oscillator

    :Examples
    ========

    >>> from sympy.physics.qho_1d import psi_n
    >>> from sympy import var
    >>> var("x m omega")
    (x, m, omega)
    >>> psi_n(0, x, m, omega)
    (m*omega)**(1/4)*exp(-m*omega*x**2/(2*hbar))/(hbar**(1/4)*pi**(1/4))

    """

    # sympify arguments
    n, x, m, omega = map(S, [n, x, m, omega])
    nu = m * omega / hbar
    # normalization coefficient
    C =  (nu/pi)**(S(1)/4) * sqrt(1/(2**n*factorial(n)))

    return C * exp(-nu* x**2 /2) * hermite(n, sqrt(nu)*x)

    def rank(self):
        """
        Returns the lexicographic rank of the permutation.

        Examples:
        >>> from sympy.combinatorics.permutations import Permutation
        >>> p = Permutation([0,1,2,3])
        >>> p.rank
        0
        >>> p = Permutation([3,2,1,0])
        >>> p.rank
        23
        """
        rank = 0
        rho = self.array_form[:]
        n = self.size - 1
        psize = int(factorial(n))
        for j in xrange(self.size - 1):
            rank += rho[j]*psize
            for i in xrange(j + 1, self.size):
                if rho[i] > rho[j]:
                    rho[i] -= 1
            psize //= n
            n -= 1
        return rank

    def unrank_nonlex(self, n, r):
        """
        This is a linear time unranking algorithm that does not
        respect lexicographic order [3].

        [3] See below

        Examples:
        >>> from sympy.combinatorics.permutations import Permutation
        >>> Permutation.unrank_nonlex(4, 5)
        Permutation([2, 0, 3, 1])
        >>> Permutation.unrank_nonlex(4, -1)
        Permutation([0, 1, 2, 3])

        >>> Permutation.unrank_nonlex(6, 2)
        Permutation([1, 5, 3, 4, 0, 2])
        """
        def unrank1(n, r, a):
            if n > 0:
                a[n - 1], a[r % n] = a[r % n], a[n - 1]
                unrank1(n - 1, r//n, a)

        id_perm = range(n)
        n = int(n)
        r = int(r % factorial(n))
        unrank1(n, r, id_perm)
        return Permutation(id_perm)

    def _eval_rewrite_as_Sum(self, ap, bq, z):
        from sympy.functions import factorial, RisingFactorial, Piecewise

        n = C.Dummy("n", integer=True)
        rfap = Tuple(*[RisingFactorial(a, n) for a in ap])
        rfbq = Tuple(*[RisingFactorial(b, n) for b in bq])
        coeff = Mul(*rfap) / Mul(*rfbq)
        return Piecewise((C.Sum(coeff * z ** n / factorial(n), (n, 0, oo)), self.convergence_statement), (self, True))

def R_nl(n, l, nu, r):
    """
    Returns the radial wavefunction R_{nl} for a 3d isotropic harmonic
    oscillator.

    ``n``
        the "nodal" quantum number.  Corresponds to the number of nodes in
        the wavefunction.  n >= 0
    ``l``
        the quantum number for orbital angular momentum
    ``nu``
        mass-scaled frequency: nu = m*omega/(2*hbar) where `m` is the mass
        and `omega` the frequency of the oscillator.
        (in atomic units nu == omega/2)
    ``r``
        Radial coordinate

    Examples
    ========

    >>> from sympy.physics.sho import R_nl
    >>> from sympy import var
    >>> var("r nu l")
    (r, nu, l)
    >>> R_nl(0, 0, 1, r)
    2*2**(3/4)*exp(-r**2)/pi**(1/4)
    >>> R_nl(1, 0, 1, r)
    4*2**(1/4)*sqrt(3)*(-2*r**2 + 3/2)*exp(-r**2)/(3*pi**(1/4))

    l, nu and r may be symbolic:

    >>> R_nl(0, 0, nu, r)
    2*2**(3/4)*sqrt(nu**(3/2))*exp(-nu*r**2)/pi**(1/4)
    >>> R_nl(0, l, 1, r)
    r**l*sqrt(2**(l + 3/2)*2**(l + 2)/factorial2(2*l + 1))*exp(-r**2)/pi**(1/4)

    The normalization of the radial wavefunction is:

    >>> from sympy import Integral, oo
    >>> Integral(R_nl(0, 0, 1, r)**2 * r**2, (r, 0, oo)).n()
    1.00000000000000
    >>> Integral(R_nl(1, 0, 1, r)**2 * r**2, (r, 0, oo)).n()
    1.00000000000000
    >>> Integral(R_nl(1, 1, 1, r)**2 * r**2, (r, 0, oo)).n()
    1.00000000000000

    """
    n, l, nu, r = map(S, [n, l, nu, r])

    # formula uses n >= 1 (instead of nodal n >= 0)
    n = n + 1
    C = sqrt(
        ((2 * nu) ** (l + Rational(3, 2)) * 2 ** (n + l + 1) * factorial(n - 1))
        / (sqrt(pi) * (factorial2(2 * n + 2 * l - 1)))
    )
    return C * r ** (l) * exp(-nu * r ** 2) * assoc_laguerre(n - 1, l + S(1) / 2, 2 * nu * r ** 2)

def monomial_count(V, N):
    """Computes the number of monomials of degree N in #V variables.

       The number of monomials is given as (#V + N)! / (#V! N!), eg:

       >>> from sympy import *
       >>> x,y = symbols('xy')

       >>> monomial_count(2, 2)
       6

       >>> M = monomials([x, y], 2)

       >>> print M
       [1, x, y, x**2, y**2, x*y]
       >>> len(M)
       6

    """
    return factorial(V + N) / factorial(V) / factorial(N)

def monomial_count(V, N):
    """Computes the number of monomials of degree `N` in `#V` variables.

       The number of monomials is given as `(#V + N)! / (#V! N!)`, e.g.::

           >>> from sympy import monomials, monomial_count
           >>> from sympy.abc import x, y

           >>> monomial_count(2, 2)
           6

           >>> M = monomials([x, y], 2)

           >>> sorted(M)
           [1, x, y, x**2, y**2, x*y]
           >>> len(M)
           6

    """
    return factorial(V + N) / factorial(V) / factorial(N)

def monomial_count(V, N):
    """Computes the number of monomials of degree N in #V variables.

       The number of monomials is given as (#V + N)! / (#V! N!), e.g.:

       >>> from sympy.polys.monomial import monomial_count, monomials
       >>> from sympy.abc import x, y

       >>> monomial_count(2, 2)
       6

       >>> M = monomials([x, y], 2)

       >>> print M
       set([1, x, y, x*y, x**2, y**2])
       >>> len(M)
       6

    """
    return factorial(V + N) / factorial(V) / factorial(N)

def coherent_state(n, alpha):
    """
    Returns <n|alpha> for the coherent states of 1D harmonic oscillator.
    See http://en.wikipedia.org/wiki/Coherent_states

    ``n``
        the "nodal" quantum number
    ``alpha``
        the eigen value of annihilation operator
    """

    return exp(- Abs(alpha)**2/2)*(alpha**n)/sqrt(factorial(n))

    def unrank_lex(self, size, rank):
        """
        Lexicographic permutation unranking.

        Examples:
        >>> from sympy.combinatorics.permutations import Permutation
        >>> a = Permutation.unrank_lex(5,10)
        >>> a.rank
        10
        >>> a
        Permutation([0, 2, 4, 1, 3])
        """
        perm_array = [0] * size
        for i in xrange(size):
            d = (rank % int(factorial(i + 1))) / int(factorial(i))
            rank = rank - d*int(factorial(i))
            perm_array[size - i - 1] = d
            for j in xrange(size - i, size):
                if perm_array[j] > d-1:
                    perm_array[j] += 1
        return Permutation(perm_array)

def test_Function():
    assert mcode(sin(x) ** cos(x)) == "sin(x).^cos(x)"
    assert mcode(sign(x)) == "sign(x)"
    assert mcode(exp(x)) == "exp(x)"
    assert mcode(log(x)) == "log(x)"
    assert mcode(factorial(x)) == "factorial(x)"
    assert mcode(floor(x)) == "floor(x)"
    assert mcode(atan2(y, x)) == "atan2(y, x)"
    assert mcode(beta(x, y)) == 'beta(x, y)'
    assert mcode(polylog(x, y)) == 'polylog(x, y)'
    assert mcode(harmonic(x)) == 'harmonic(x)'
    assert mcode(bernoulli(x)) == "bernoulli(x)"
    assert mcode(bernoulli(x, y)) == "bernoulli(x, y)"

    def euler_maclaurin(self, n=0):
        """
        Return n-th order Euler-Maclaurin approximation of self.

        The 0-th order approximation is simply the corresponding
        integral
        """
        f, i, a, b = self.f, self.i, self.a, self.b
        x = Symbol('x', dummy=True)
        s = Basic.Integral(f.subs(i, x), (x, a, b)).doit()
        if n > 0:
            s += (f.subs(i, a) + f.subs(i, b))/2
        for k in range(1, n):
            g = f.diff(i, 2*k-1)
            B = Basic.bernoulli
            s += B(2*k)/factorial(2*k)*(g.subs(i,b)-g.subs(i,a))
        return s

def test_sympy_parser():
    x = Symbol('x')
    inputs = {
        '2*x': 2 * x,
        '3.00': Float(3),
        '22/7': Rational(22, 7),
        '2+3j': 2 + 3*I,
        'exp(x)': exp(x),
        'x!': factorial(x),
        '3.[3]': Rational(10, 3),
        '10!': 3628800,
        '(_kern)': Symbol('_kern'),
        '-(2)': -Integer(2),
        '[-1, -2, 3]': [Integer(-1), Integer(-2), Integer(3)]
    }
    for text, result in inputs.items():
        assert parse_expr(text) == result

def test_sympy_parser():
    x = Symbol("x")
    inputs = {
        "2*x": 2 * x,
        "3.00": Float(3),
        "22/7": Rational(22, 7),
        "2+3j": 2 + 3 * I,
        "exp(x)": exp(x),
        "x!": factorial(x),
        "3.[3]": Rational(10, 3),
        "10!": 3628800,
        "-(2)": -Integer(2),
        "[-1, -2, 3]": [Integer(-1), Integer(-2), Integer(3)],
        'Symbol("x").free_symbols': x.free_symbols,
        "S('S(3).n(n=3)')": 3.00,
        "factorint(12, visual=True)": Mul(Pow(2, 2, evaluate=False), Pow(3, 1, evaluate=False), evaluate=False),
        'Limit(sin(x), x, 0, dir="-")': Limit(sin(x), x, 0, dir="-"),
    }
    for text, result in inputs.items():
        assert parse_expr(text) == result

def test_sympy_parser():
    x = Symbol('x')
    inputs = {
        '2*x': 2 * x,
        '3.00': Float(3),
        '22/7': Rational(22, 7),
        '2+3j': 2 + 3*I,
        'exp(x)': exp(x),
        'x!': factorial(x),
        '3.[3]': Rational(10, 3),
        '10!': 3628800,
        '-(2)': -Integer(2),
        '[-1, -2, 3]': [Integer(-1), Integer(-2), Integer(3)],
        'Symbol("x").free_symbols': x.free_symbols,
        "S('S(3).n(n=3)')": 3.00,
        'factorint(12, visual=True)': Mul(
            Pow(2, 2, evaluate=False),
            Pow(3, 1, evaluate=False),
            evaluate=False),
        'Limit(sin(x), x, 0, dir="-")': Limit(sin(x), x, 0, dir='-'),

    }
    for text, result in inputs.items():
        assert parse_expr(text) == result

    def euler_maclaurin(self, m=0, n=0, eps=0, eval_integral=True):
        """
        Return an Euler-Maclaurin approximation of self, where m is the
        number of leading terms to sum directly and n is the number of
        terms in the tail.

        With m = n = 0, this is simply the corresponding integral
        plus a first-order endpoint correction.

        Returns (s, e) where s is the Euler-Maclaurin approximation
        and e is the estimated error (taken to be the magnitude of
        the first omitted term in the tail):

            >>> from sympy.abc import k, a, b
            >>> from sympy import Sum
            >>> Sum(1/k, (k, 2, 5)).doit().evalf()
            1.28333333333333
            >>> s, e = Sum(1/k, (k, 2, 5)).euler_maclaurin()
            >>> s
            -log(2) + 7/20 + log(5)
            >>> from sympy import sstr
            >>> print(sstr((s.evalf(), e.evalf()), full_prec=True))
            (1.26629073187415, 0.0175000000000000)

        The endpoints may be symbolic:

            >>> s, e = Sum(1/k, (k, a, b)).euler_maclaurin()
            >>> s
            -log(a) + log(b) + 1/(2*b) + 1/(2*a)
            >>> e
            Abs(1/(12*b**2) - 1/(12*a**2))

        If the function is a polynomial of degree at most 2n+1, the
        Euler-Maclaurin formula becomes exact (and e = 0 is returned):

            >>> Sum(k, (k, 2, b)).euler_maclaurin()
            (b**2/2 + b/2 - 1, 0)
            >>> Sum(k, (k, 2, b)).doit()
            b**2/2 + b/2 - 1

        With a nonzero eps specified, the summation is ended
        as soon as the remainder term is less than the epsilon.
        """
        from sympy.functions import bernoulli, factorial
        from sympy.integrals import Integral

        m = int(m)
        n = int(n)
        f = self.function
        if len(self.limits) != 1:
            raise ValueError("More than 1 limit")
        i, a, b = self.limits[0]
        if (a > b) == True:
            if a - b == 1:
                return S.Zero, S.Zero
            a, b = b + 1, a - 1
            f = -f
        s = S.Zero
        if m:
            if b.is_Integer and a.is_Integer:
                m = min(m, b - a + 1)
            if not eps or f.is_polynomial(i):
                for k in range(m):
                    s += f.subs(i, a + k)
            else:
                term = f.subs(i, a)
                if term:
                    test = abs(term.evalf(3)) < eps
                    if test == True:
                        return s, abs(term)
                    elif not (test == False):
                        # a symbolic Relational class, can't go further
                        return term, S.Zero
                s += term
                for k in range(1, m):
                    term = f.subs(i, a + k)
                    if abs(term.evalf(3)) < eps and term != 0:
                        return s, abs(term)
                    s += term
            if b - a + 1 == m:
                return s, S.Zero
            a += m
        x = Dummy('x')
        I = Integral(f.subs(i, x), (x, a, b))
        if eval_integral:
            I = I.doit()
        s += I

        def fpoint(expr):
            if b is S.Infinity:
                return expr.subs(i, a), 0
            return expr.subs(i, a), expr.subs(i, b)
        fa, fb = fpoint(f)
        iterm = (fa + fb)/2
        g = f.diff(i)
        for k in range(1, n + 2):
            ga, gb = fpoint(g)
            term = bernoulli(2*k)/factorial(2*k)*(gb - ga)
            if (eps and term and abs(term.evalf(3)) < eps) or (k > n):
                break
#.........这里部分代码省略.........

def test_rcode_functions():
    assert rcode(sin(x) ** cos(x)) == "sin(x)^cos(x)"
    assert rcode(factorial(x) + gamma(y)) == "factorial(x) + gamma(y)"
    assert rcode(beta(Min(x, y), Max(x, y))) == "beta(min(x, y), max(x, y))"

def test_harmonic_rewrite_polygamma():
    n = Symbol("n")
    m = Symbol("m")

    assert harmonic(n).rewrite(digamma) == polygamma(0, n + 1) + EulerGamma
    assert harmonic(n).rewrite(trigamma) ==  polygamma(0, n + 1) + EulerGamma
    assert harmonic(n).rewrite(polygamma) ==  polygamma(0, n + 1) + EulerGamma

    assert harmonic(n,3).rewrite(polygamma) == polygamma(2, n + 1)/2 - polygamma(2, 1)/2
    assert harmonic(n,m).rewrite(polygamma) == (-1)**m*(polygamma(m - 1, 1) - polygamma(m - 1, n + 1))/factorial(m - 1)

    assert expand_func(harmonic(n+4)) == harmonic(n) + 1/(n + 4) + 1/(n + 3) + 1/(n + 2) + 1/(n + 1)
    assert expand_func(harmonic(n-4)) == harmonic(n) - 1/(n - 1) - 1/(n - 2) - 1/(n - 3) - 1/n

    assert harmonic(n, m).rewrite("tractable") == harmonic(n, m).rewrite(polygamma)