def taylor_term(n, x, *previous_terms):
     if n < 0:
         return S.Zero
     else:
         x = sympify(x)
         if len(previous_terms) > 1:
             p = previous_terms[-1]
             return (3**(S(1)/3)*x * Abs(sin(2*pi*(n + S.One)/S(3))) * C.factorial((n - S.One)/S(3)) /
                     ((n + S.One) * Abs(cos(2*pi*(n + S.Half)/S(3))) * C.factorial((n - 2)/S(3))) * p)
         else:
             return (S.One/(root(3, 6)*pi) * gamma((n + S.One)/S(3)) * Abs(sin(2*pi*(n + S.One)/S(3))) /
                     C.factorial(n) * (root(3, 3)*x)**n)

 def _calc_bernoulli(n):
     s = 0
     a = int(C.binomial(n+3, n-6))
     for j in xrange(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 / C.binomial(n+3, n)

def test_subfactorial():
    assert all(subfactorial(i) == ans for i, ans in enumerate(
        [1, 0, 1, 2, 9, 44, 265, 1854, 14833, 133496]))
    assert subfactorial(oo) == oo

    x = Symbol('x')
    assert subfactorial(x).rewrite(C.uppergamma) == \
        C.uppergamma(x + 1, -1)/S.Exp1

    tt = Symbol('tt', integer=True, nonnegative=True)
    tf = Symbol('tf', integer=True, nonnegative=False)
    tn = Symbol('tf', integer=True)
    ft = Symbol('ft', integer=False, nonnegative=True)
    ff = Symbol('ff', integer=False, nonnegative=False)
    fn = Symbol('ff', integer=False)
    nt = Symbol('nt', nonnegative=True)
    nf = Symbol('nf', nonnegative=False)
    nn = Symbol('nf')
    assert subfactorial(tt).is_integer
    assert subfactorial(tf).is_integer is None
    assert subfactorial(tn).is_integer is None
    assert subfactorial(ft).is_integer is None
    assert subfactorial(ff).is_integer is None
    assert subfactorial(fn).is_integer is None
    assert subfactorial(nt).is_integer is None
    assert subfactorial(nf).is_integer is None
    assert subfactorial(nn).is_integer is None

    def _eval_rewrite_as_Sum(self, arg):
        if arg.is_even:
            k = C.Dummy("k", integer=True)
            j = C.Dummy("j", integer=True)
            n = self.args[0] / 2
            Em = (S.ImaginaryUnit * C.Sum( C.Sum( C.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 Pow(expr, assumptions):
        """
        Real**Integer              -> Real
        Positive**Real             -> Real
        Real**(Integer/Even)       -> Real if base is nonnegative
        Real**(Integer/Odd)        -> Real
        Imaginary**(Integer/Even)  -> Real
        Imaginary**(Integer/Odd)   -> not Real
        Imaginary**Real            -> ? since Real could be 0 (giving real) or 1 (giving imaginary)
        b**Imaginary               -> Real if log(b) is imaginary and b != 0 and exponent != integer multiple of I*pi/log(b)
        Real**Real                 -> ? e.g. sqrt(-1) is imaginary and sqrt(2) is not
        """
        if expr.is_number:
            return AskRealHandler._number(expr, assumptions)

        if expr.base.func == C.exp:
            if ask(Q.imaginary(expr.base.args[0]), assumptions):
                if ask(Q.imaginary(expr.exp), assumptions):
                    return True
            # If the i = (exp's arg)/(I*pi) is an integer or half-integer
            # multiple of I*pi then 2*i will be an integer. In addition,
            # exp(i*I*pi) = (-1)**i so the overall realness of the expr
            # can be determined by replacing exp(i*I*pi) with (-1)**i.
            i = expr.base.args[0]/I/pi
            if ask(Q.integer(2*i), assumptions):
                return ask(Q.real(((-1)**i)**expr.exp), assumptions)
            return

        if ask(Q.imaginary(expr.base), assumptions):
            if ask(Q.integer(expr.exp), assumptions):
                odd = ask(Q.odd(expr.exp), assumptions)
                if odd is not None:
                    return not odd
                return

        if ask(Q.imaginary(expr.exp), assumptions):
            imlog = ask(Q.imaginary(C.log(expr.base)), assumptions)
            if imlog is not None:
                # I**i -> real, log(I) is imag;
                # (2*I)**i -> complex, log(2*I) is not imag
                return imlog

        if ask(Q.real(expr.base), assumptions):
            if ask(Q.real(expr.exp), assumptions):
                if expr.exp.is_Rational and \
                        ask(Q.even(expr.exp.q), assumptions):
                    return ask(Q.positive(expr.base), assumptions)
                elif ask(Q.integer(expr.exp), assumptions):
                    return True
                elif ask(Q.positive(expr.base), assumptions):
                    return True
                elif ask(Q.negative(expr.base), assumptions):
                    return False

def test_curve():
    s = Symbol("s")
    z = Symbol("z")

    # this curve is independent of the indicated parameter
    C = Curve([2 * s, s ** 2], (z, 0, 2))

    assert C.parameter == z
    assert C.functions == (2 * s, s ** 2)
    assert C.arbitrary_point() == Point(2 * s, s ** 2)
    assert C.arbitrary_point(z) == Point(2 * s, s ** 2)

    # this is how it is normally used
    C = Curve([2 * s, s ** 2], (s, 0, 2))

    assert C.parameter == s
    assert C.functions == (2 * s, s ** 2)
    t = Symbol("t")
    assert C.arbitrary_point() != Point(2 * t, t ** 2)  # the t returned as assumptions
    t = Symbol("t", real=True)  # now t has the same assumptions so the test passes
    assert C.arbitrary_point() == Point(2 * t, t ** 2)
    assert C.arbitrary_point(z) == Point(2 * z, z ** 2)
    assert C.arbitrary_point(C.parameter) == Point(2 * s, s ** 2)

    raises(ValueError, "Curve((s, s + t), (s, 1, 2)).arbitrary_point()")
    raises(ValueError, "Curve((s, s + t), (t, 1, 2)).arbitrary_point(s)")

 def eval(cls, n, m=None):
     if m is None:
         m = S.One
     if n == oo:
         return C.zeta(m)
     if n.is_Integer and n.is_nonnegative and m.is_Integer:
         if n == 0:
             return S.Zero
         if not m in cls._functions:
             @recurrence_memo([0])
             def f(n, prev):
                 return prev[-1] + S.One / n**m
             cls._functions[m] = f
         return cls._functions[m](int(n))

    def Pow(expr, assumptions):
        """
        Imaginary**integer/odd  -> Imaginary
        Imaginary**integer/even -> Real if integer % 2 == 0
        b**Imaginary            -> !Imaginary if exponent is an integer multiple of I*pi/log(b)
        Imaginary**Real         -> ?
        Negative**even root     -> Imaginary
        Negative**odd root      -> Real
        Negative**Real          -> Imaginary
        Real**Integer           -> Real
        Real**Positive          -> Real
        """
        if expr.is_number:
            return AskImaginaryHandler._number(expr, assumptions)

        if expr.base.func == C.exp:
            if ask(Q.imaginary(expr.base.args[0]), assumptions):
                if ask(Q.imaginary(expr.exp), assumptions):
                    return False
                i = expr.base.args[0]/I/pi
                if ask(Q.integer(2*i), assumptions):
                    return ask(Q.imaginary(((-1)**i)**expr.exp), assumptions)

        if ask(Q.imaginary(expr.base), assumptions):
            if ask(Q.integer(expr.exp), assumptions):
                odd = ask(Q.odd(expr.exp), assumptions)
                if odd is not None:
                    return odd
                return

        if ask(Q.imaginary(expr.exp), assumptions):
            imlog = ask(Q.imaginary(C.log(expr.base)), assumptions)
            if imlog is not None:
                return False  # I**i -> real; (2*I)**i -> complex ==> not imaginary

        if ask(Q.real(expr.base), assumptions):
            if ask(Q.real(expr.exp), assumptions):
                if ask(Q.rational(expr.exp) & Q.even(denom(expr.exp)), assumptions):
                    return ask(Q.negative(expr.base), assumptions)
                elif ask(Q.integer(expr.exp), assumptions):
                    return False
                elif ask(Q.positive(expr.base), assumptions):
                    return False
                elif ask(Q.negative(expr.base), assumptions):
                    return True

def Ylm(l, m, theta, phi):
    """
    Spherical harmonics Ylm.

    Examples:

        >>> from sympy import symbols, Ylm
        >>> theta, phi = symbols("theta phi")
        >>> Ylm(0, 0, theta, phi)
        1/(2*sqrt(pi))
        >>> Ylm(1, -1, theta, phi)
        sqrt(6)*exp(-I*phi)*sin(theta)/(4*sqrt(pi))
        >>> Ylm(1, 0, theta, phi)
        sqrt(3)*cos(theta)/(2*sqrt(pi))

    """
    l, m, theta, phi = [sympify(x) for x in (l, m, theta, phi)]
    factorial = C.factorial
    return sqrt((2*l+1)/(4*pi) * factorial(l-m)/factorial(l+m)) * \
            Plmcos(l, m, theta) * C.exp(I*m*phi)

 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 xrange(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 xrange(n + 1):
                     result.append(C.binomial(n, k)*cls(k)*sym**(n - k))
                 return Add(*result)
         else:
             raise ValueError("Bernoulli numbers are defined only"
                              " for nonnegative integer indices.")

def _minpoly_sin(ex, x):
    """
    Returns the minimal polynomial of ``sin(ex)``
    see http://mathworld.wolfram.com/TrigonometryAngles.html
    """
    from sympy.functions.combinatorial.factorials import binomial
    c, a = ex.args[0].as_coeff_Mul()
    if a is pi:
        if c.is_rational:
            n = c.q
            q = sympify(n)
            if q.is_prime:
                # for a = pi*p/q with q odd prime, using chebyshevt
                # write sin(q*a) = mp(sin(a))*sin(a);
                # the roots of mp(x) are sin(pi*p/q) for p = 1,..., q - 1
                a = dup_chebyshevt(n, ZZ)
                return Add(*[x**(n - i - 1)*a[i] for i in range(n)])
            if c.p == 1:
                if q == 9:
                    return 64*x**6 - 96*x**4 + 36*x**2 - 3

            if n % 2 == 1:
                # for a = pi*p/q with q odd, use
                # sin(q*a) = 0 to see that the minimal polynomial must be
                # a factor of dup_chebyshevt(n, ZZ)
                a = dup_chebyshevt(n, ZZ)
                a = [x**(n - i)*a[i] for i in range(n + 1)]
                r = Add(*a)
                _, factors = factor_list(r)
                res = _choose_factor(factors, x, ex)
                return res

            expr = ((1 - C.cos(2*c*pi))/2)**S.Half
            res = _minpoly_compose(expr, x, QQ)
            return res

    raise NotAlgebraic("%s doesn't seem to be an algebraic element" % ex)

def test_line():
    p1 = Point(0, 0)
    p2 = Point(1, 1)
    p3 = Point(x1, x1)
    p4 = Point(y1, y1)
    p5 = Point(x1, 1 + x1)
    p6 = Point(1, 0)
    p7 = Point(0, 1)
    p8 = Point(2, 0)
    p9 = Point(2, 1)

    l1 = Line(p1, p2)
    l2 = Line(p3, p4)
    l3 = Line(p3, p5)
    l4 = Line(p1, p6)
    l5 = Line(p1, p7)
    l6 = Line(p8, p9)
    l7 = Line(p2, p9)
    raises(ValueError, lambda: Line(Point(0, 0), Point(0, 0)))

    # Basic stuff
    assert Line((1, 1), slope=1) == Line((1, 1), (2, 2))
    assert Line((1, 1), slope=oo) == Line((1, 1), (1, 2))
    assert Line((1, 1), slope=-oo) == Line((1, 1), (1, 2))
    raises(ValueError, lambda: Line((1, 1), 1))
    assert Line(p1, p2) == Line(p2, p1)
    assert l1 == l2
    assert l1 != l3
    assert l1.slope == 1
    assert l1.length == oo
    assert l3.slope == oo
    assert l4.slope == 0
    assert l4.coefficients == (0, 1, 0)
    assert l4.equation(x=x, y=y) == y
    assert l5.slope == oo
    assert l5.coefficients == (1, 0, 0)
    assert l5.equation() == x
    assert l6.equation() == x - 2
    assert l7.equation() == y - 1
    assert p1 in l1  # is p1 on the line l1?
    assert p1 not in l3
    assert Line((-x, x), (-x + 1, x - 1)).coefficients == (1, 1, 0)

    assert simplify(l1.equation()) in (x - y, y - x)
    assert simplify(l3.equation()) in (x - x1, x1 - x)

    assert Line(p1, p2).scale(2, 1) == Line(p1, p9)

    assert l2.arbitrary_point() in l2
    for ind in xrange(0, 5):
        assert l3.random_point() in l3

    # Orthogonality
    p1_1 = Point(-x1, x1)
    l1_1 = Line(p1, p1_1)
    assert l1.perpendicular_line(p1) == l1_1
    assert Line.is_perpendicular(l1, l1_1)
    assert Line.is_perpendicular(l1, l2) == False
    p = l1.random_point()
    assert l1.perpendicular_segment(p) == p

    # Parallelity
    p2_1 = Point(-2 * x1, 0)
    l2_1 = Line(p3, p5)
    assert l2.parallel_line(p1_1) == Line(p2_1, p1_1)
    assert l2_1.parallel_line(p1) == Line(p1, Point(0, 2))
    assert Line.is_parallel(l1, l2)
    assert Line.is_parallel(l2, l3) == False
    assert Line.is_parallel(l2, l2.parallel_line(p1_1))
    assert Line.is_parallel(l2_1, l2_1.parallel_line(p1))

    # Intersection
    assert intersection(l1, p1) == [p1]
    assert intersection(l1, p5) == []
    assert intersection(l1, l2) in [[l1], [l2]]
    assert intersection(l1, l1.parallel_line(p5)) == []

    # Concurrency
    l3_1 = Line(Point(5, x1), Point(-Rational(3, 5), x1))
    assert Line.is_concurrent(l1) == False
    assert Line.is_concurrent(l1, l3)
    assert Line.is_concurrent(l1, l3, l3_1)
    assert Line.is_concurrent(l1, l1_1, l3) == False

    # Projection
    assert l2.projection(p4) == p4
    assert l1.projection(p1_1) == p1
    assert l3.projection(p2) == Point(x1, 1)
    raises(GeometryError, lambda: Line(Point(0, 0), Point(1, 0)).projection(Circle(Point(0, 0), 1)))

    # Finding angles
    l1_1 = Line(p1, Point(5, 0))
    assert feq(Line.angle_between(l1, l1_1).evalf(), pi.evalf() / 4)

    # Testing Rays and Segments (very similar to Lines)
    assert Ray((1, 1), angle=pi / 4) == Ray((1, 1), (2, 2))
    assert Ray((1, 1), angle=pi / 2) == Ray((1, 1), (1, 2))
    assert Ray((1, 1), angle=-pi / 2) == Ray((1, 1), (1, 0))
    assert Ray((1, 1), angle=-3 * pi / 2) == Ray((1, 1), (1, 2))
    assert Ray((1, 1), angle=5 * pi / 2) == Ray((1, 1), (1, 2))
#.........这里部分代码省略.........

 def _eval_rewrite_as_hyper(self,n):
     return C.hyper([1-n,-n],[2],1)

 def _eval_rewrite_as_gamma(self,n):
     # The gamma function allows to generalize Catalan numbers to complex n
     return 4**n*C.gamma(n + S.Half)/(C.gamma(S.Half)*C.gamma(n+2))

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

 def fdiff(self, argindex=1):
     n = self.args[0]
     return catalan(n)*(C.polygamma(0,n+Rational(1,2))-C.polygamma(0,n+2)+C.log(4))

def Ylm(l, m, theta, phi):
    l, m, theta, phi = [sympify(x) for x in (l, m, theta, phi)]
    factorial = C.Factorial
    return sqrt((2 * l + 1) / (4 * pi) * factorial(l - m) / factorial(l + m)) * Plmcos(l, m, theta) * C.exp(I * m * phi)

def test_line():
    p1 = Point(0, 0)
    p2 = Point(1, 1)
    p3 = Point(x1, x1)
    p4 = Point(y1, y1)
    p5 = Point(x1, 1 + x1)
    p6 = Point(1, 0)
    p7 = Point(0, 1)
    p8 = Point(2, 0)
    p9 = Point(2, 1)

    l1 = Line(p1, p2)
    l2 = Line(p3, p4)
    l3 = Line(p3, p5)
    l4 = Line(p1, p6)
    l5 = Line(p1, p7)
    l6 = Line(p8, p9)
    l7 = Line(p2, p9)

    # Basic stuff
    assert Line((1, 1), slope=1) == Line((1, 1), (2, 2))
    assert Line((1, 1), slope=oo) == Line((1, 1), (1, 2))
    assert Line((1, 1), slope=-oo) == Line((1, 1), (1, 2))
    raises(ValueError, "Line((1, 1), 1)")
    assert Line(p1, p2) == Line(p2, p1)
    assert l1 == l2
    assert l1 != l3
    assert l1.slope == 1
    assert l3.slope == oo
    assert l4.slope == 0
    assert l4.coefficients == (0, 1, 0)
    assert l4.equation(x=x, y=y) == y
    assert l5.slope == oo
    assert l5.coefficients == (1, 0, 0)
    assert l5.equation() == x
    assert l6.equation() == x - 2
    assert l7.equation() == y - 1
    assert p1 in l1  # is p1 on the line l1?
    assert p1 not in l3

    assert simplify(l1.equation()) in (x - y, y - x)
    assert simplify(l3.equation()) in (x - x1, x1 - x)

    assert l2.arbitrary_point() in l2
    for ind in xrange(0, 5):
        assert l3.random_point() in l3

    # Orthogonality
    p1_1 = Point(-x1, x1)
    l1_1 = Line(p1, p1_1)
    assert l1.perpendicular_line(p1) == l1_1
    assert Line.is_perpendicular(l1, l1_1)
    assert Line.is_perpendicular(l1, l2) == False

    # Parallelity
    p2_1 = Point(-2 * x1, 0)
    l2_1 = Line(p3, p5)
    assert l2.parallel_line(p1_1) == Line(p2_1, p1_1)
    assert l2_1.parallel_line(p1) == Line(p1, Point(0, 2))
    assert Line.is_parallel(l1, l2)
    assert Line.is_parallel(l2, l3) == False
    assert Line.is_parallel(l2, l2.parallel_line(p1_1))
    assert Line.is_parallel(l2_1, l2_1.parallel_line(p1))

    # Intersection
    assert intersection(l1, p1) == [p1]
    assert intersection(l1, p5) == []
    assert intersection(l1, l2) in [[l1], [l2]]
    assert intersection(l1, l1.parallel_line(p5)) == []

    # Concurrency
    l3_1 = Line(Point(5, x1), Point(-Rational(3, 5), x1))
    assert Line.is_concurrent(l1, l3)
    assert Line.is_concurrent(l1, l3, l3_1)
    assert Line.is_concurrent(l1, l1_1, l3) == False

    # Projection
    assert l2.projection(p4) == p4
    assert l1.projection(p1_1) == p1
    assert l3.projection(p2) == Point(x1, 1)

    # Finding angles
    l1_1 = Line(p1, Point(5, 0))
    assert feq(Line.angle_between(l1, l1_1).evalf(), pi.evalf() / 4)

    # Testing Rays and Segments (very similar to Lines)
    assert Ray((1, 1), angle=pi / 4) == Ray((1, 1), (2, 2))
    assert Ray((1, 1), angle=pi / 2) == Ray((1, 1), (1, 2))
    assert Ray((1, 1), angle=-pi / 2) == Ray((1, 1), (1, 0))
    assert Ray((1, 1), angle=-3 * pi / 2) == Ray((1, 1), (1, 2))
    assert Ray((1, 1), angle=5 * pi / 2) == Ray((1, 1), (1, 2))
    assert Ray((1, 1), angle=5.0 * pi / 2) == Ray((1, 1), (1, 2))
    assert Ray((1, 1), angle=pi) == Ray((1, 1), (0, 1))
    assert Ray((1, 1), angle=3.0 * pi) == Ray((1, 1), (0, 1))
    assert Ray((1, 1), angle=4.0 * pi) == Ray((1, 1), (2, 1))
    assert Ray((1, 1), angle=0) == Ray((1, 1), (2, 1))
    # XXX don't know why this fails without str
    assert str(Ray((1, 1), angle=4.2 * pi)) == str(Ray(Point(1, 1), Point(2, 1 + C.tan(0.2 * pi))))
    assert Ray((1, 1), angle=5) == Ray((1, 1), (2, 1 + C.tan(5)))
    raises(ValueError, "Ray((1, 1), 1)")
#.........这里部分代码省略.........