def _eval_rewrite_as_polynomial(self, n, a, b, x):
     from sympy import Sum
     # TODO: Make sure n \in N
     k = Dummy("k")
     kern = (RisingFactorial(-n, k) * RisingFactorial(a + b + n + 1, k) * RisingFactorial(a + k + 1, n - k) /
             factorial(k) * ((1 - x)/2)**k)
     return 1 / factorial(n) * Sum(kern, (k, 0, n))

 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) ** (-n)
                 * (3 ** (S(1) / 3) * x) ** (n + 1)
                 * sin(pi * (2 * n / 3 + S(4) / 3))
                 * factorial(n)
                 * gamma(n / 3 + S(2) / 3)
                 / (sin(pi * (2 * n / 3 + S(2) / 3)) * factorial(n + 1) * gamma(n / 3 + S(1) / 3))
                 * p
             )
         else:
             return (
                 S.One
                 / (3 ** (S(2) / 3) * pi)
                 * gamma((n + S.One) / S(3))
                 * sin(2 * pi * (n + S.One) / S(3))
                 / factorial(n)
                 * (root(3, 3) * x) ** n
             )

def wedge(first_tensor, second_tensor):
    """
    Finds outer product (wedge) of two tensor arguments.
    The algoritm is too find the asymmetric form of tensor product of
    two arguments. The resulted array is multiplied on coefficient which is
    coeff = factorial(p+s)/factorial(p)*factorial(s)

    Examples
    ========
    from tensor_analysis.arraypy import Arraypy, Tensor
    from tensor_analysis.tensor_methods import wedge
    >>> a = Arraypy((3,), 'A').to_tensor((-1))
    >>> b = Arraypy((3,), 'B').to_tensor((1,))
    >>> c = wedge(a, b)
    >>> print(c)
    0  10*A[0]*B[1] - 10*A[1]*B[0]  10*A[0]*B[2] - 10*A[2]*B[0]  
    -10*A[0]*B[1] + 10*A[1]*B[0]  0  10*A[1]*B[2] - 10*A[2]*B[1]  
    -10*A[0]*B[2] + 10*A[2]*B[0]  -10*A[1]*B[2] + 10*A[2]*B[1]  0 


    """
    if not isinstance(first_tensor, TensorArray):
        raise TypeError('Input must be of TensorArray type')
    if not isinstance(second_tensor, TensorArray):
        raise TypeError('Input must be of TensorArray type')

    p = len(first_tensor)
    s = len(second_tensor)

    coeff = factorial(p + s) / (factorial(p) * factorial(s))
    return coeff * asymmetric(tensor_product(first_tensor, second_tensor))

 def calc2(n,m):
     if m >= 0:
         return assoc_legendre._calc2(n,m)
     else:
         factorial = C.Factorial
         m = -m
         return (-1)**m *factorial(n-m)/factorial(n+m) * assoc_legendre._calc2(n, m)

def asymmetric(in_arr):
    """
    Creates the asymmetric form of input tensor.
    Input: Arraypy or TensorArray with equal axes (array shapes).
    Output: asymmetric array. Output type - Arraypy or TensorArray, depends of input

    Examples
    ========

    >>> from sympy.tensor.arraypy import list2arraypy
    >>> from sympy.tensor.tensor_methods import asymmetric
    >>> a = list2arraypy(list(range(9)), (3,3))
    >>> b = asymmetric(a)
    >>> print (b)
    0  -1.00000000000000  -2.00000000000000
    1.00000000000000  0  -1.00000000000000
    2.00000000000000  1.00000000000000  0
    """
    if not isinstance(in_arr, Arraypy):
        raise TypeError('Input must be Arraypy or TensorArray type')

    flag = 0
    for j in range(in_arr.rank):
        if (in_arr.shape[0] != in_arr.shape[j]):
            raise ValueError('Different size of arrays axes')

    # forming list of tuples for Arraypy constructor of type a = Arraypy( [(a,
    # b), (c, d), ... , (y, z)] )
    arg = [(in_arr.start_index[i], in_arr.end_index[i])
           for i in range(in_arr.rank)]

    # if in_arr TensorArray, then res_arr will be TensorArray, else it will be
    # Arraypy
    if isinstance(in_arr, TensorArray):
        res_arr = TensorArray(Arraypy(arg), in_arr.ind_char)
    else:
        res_arr = Arraypy(arg)

    signs = [0 for i in range(factorial(in_arr.rank))]
    temp_i = 0
    for p in permutations(range(in_arr.rank)):
        signs[temp_i] = perm_parity(list(p))
        temp_i += 1

    index = [in_arr.start_index[i] for i in range(in_arr.rank)]

    for i in range(len(in_arr)):
        perm = list(permutations(index))
        perm_number = 0
        for temp_index in perm:
            res_arr[tuple(index)] += signs[perm_number] * \
                in_arr[tuple(temp_index)]
            perm_number += 1
        if isinstance(res_arr[tuple(index)], int):
            res_arr[tuple(index)] = float(res_arr[tuple(index)])
        res_arr[tuple(index)] /= factorial(in_arr.rank)

        index = in_arr.next_index(index)

    return res_arr

 def as_real_imag(self, deep=True, **hints):
     # TODO: Handle deep and hints
     n, m, theta, phi = self.args
     re = (sqrt((2*n + 1)/(4*pi) * factorial(n - m)/factorial(n + m)) *
           cos(m*phi) * assoc_legendre(n, m, cos(theta)))
     im = (sqrt((2*n + 1)/(4*pi) * factorial(n - m)/factorial(n + m)) *
           sin(m*phi) * assoc_legendre(n, m, cos(theta)))
     return (re, im)

 def _eval_rewrite_as_polynomial(self, n, x):
     from sympy import Sum
     # Make sure n \in N_0
     if n.is_negative or n.is_integer is False:
         raise ValueError("Error: n should be a non-negative integer.")
     k = Dummy("k")
     kern = RisingFactorial(
         -n, k) / (gamma(k + alpha + 1) * factorial(k)) * x**k
     return gamma(n + alpha + 1) / factorial(n) * Sum(kern, (k, 0, n))

def interpolate(data, x):
    """
    Construct an interpolating polynomial for the data points.

    Examples
    ========

    >>> from sympy.polys.polyfuncs import interpolate
    >>> from sympy.abc import x

    A list is interpreted as though it were paired with a range starting
    from 1:

    >>> interpolate([1, 4, 9, 16], x)
    x**2

    This can be made explicit by giving a list of coordinates:

    >>> interpolate([(1, 1), (2, 4), (3, 9)], x)
    x**2

    The (x, y) coordinates can also be given as keys and values of a
    dictionary (and the points need not be equispaced):

    >>> interpolate([(-1, 2), (1, 2), (2, 5)], x)
    x**2 + 1
    >>> interpolate({-1: 2, 1: 2, 2: 5}, x)
    x**2 + 1

    """
    n = len(data)
    poly = None

    if isinstance(data, dict):
        X, Y = list(zip(*data.items()))
        poly = interpolating_poly(n, x, X, Y)
    else:
        if isinstance(data[0], tuple):
            X, Y = list(zip(*data))
            poly = interpolating_poly(n, x, X, Y)
        else:
            Y = list(data)

            numert = Mul(*[(x - i) for i in range(1, n + 1)])
            denom = -factorial(n - 1) if n%2 == 0 else factorial(n - 1)
            coeffs = []
            for i in range(1, n + 1):
                coeffs.append(numert/(x - i)/denom)
                denom = denom/(i - n)*i

            poly = Add(*[coeff*y for coeff, y in zip(coeffs, Y)])

    return poly.expand()

 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))) * factorial((n - S.One)/S(3)) /
                     ((n + S.One) * Abs(cos(2*pi*(n + S.Half)/S(3))) * 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))) /
                     factorial(n) * (root(3, 3)*x)**n)

def _nP(n, k=None, replacement=False):
    from sympy.functions.combinatorial.factorials import factorial
    from sympy.core.mul import prod

    if k == 0:
        return 1
    if isinstance(n, SYMPY_INTS):  # n different items
        # assert n >= 0
        if k is None:
            return sum(_nP(n, i, replacement) for i in range(n + 1))
        elif replacement:
            return n**k
        elif k > n:
            return 0
        elif k == n:
            return factorial(k)
        elif k == 1:
            return n
        else:
            # assert k >= 0
            return _product(n - k + 1, n)
    elif isinstance(n, _MultisetHistogram):
        if k is None:
            return sum(_nP(n, i, replacement) for i in range(n[_N] + 1))
        elif replacement:
            return n[_ITEMS]**k
        elif k == n[_N]:
            return factorial(k)/prod([factorial(i) for i in n[_M] if i > 1])
        elif k > n[_N]:
            return 0
        elif k == 1:
            return n[_ITEMS]
        else:
            # assert k >= 0
            tot = 0
            n = list(n)
            for i in range(len(n[_M])):
                if not n[i]:
                    continue
                n[_N] -= 1
                if n[i] == 1:
                    n[i] = 0
                    n[_ITEMS] -= 1
                    tot += _nP(_MultisetHistogram(n), k - 1)
                    n[_ITEMS] += 1
                    n[i] = 1
                else:
                    n[i] -= 1
                    tot += _nP(_MultisetHistogram(n), k - 1)
                    n[i] += 1
                n[_N] += 1
            return tot

 def eval(cls, n, m, x):
     if m.could_extract_minus_sign():
         # P^{-m}_n  --->  F * P^m_n
         return S.NegativeOne**(-m) * (factorial(m + n)/factorial(n - m)) * assoc_legendre(n, -m, x)
     if m == 0:
         # P^0_n  --->  L_n
         return legendre(n, x)
     if x == 0:
         return 2**m*sqrt(S.Pi) / (gamma((1 - m - n)/2)*gamma(1 - (m - n)/2))
     if n.is_Number and m.is_Number and n.is_integer and m.is_integer:
         if n.is_negative:
             raise ValueError("%s : 1st index must be nonnegative integer (got %r)" % (cls, n))
         if abs(m) > n:
             raise ValueError("%s : abs('2nd index') must be <= '1st index' (got %r, %r)" % (cls, n, m))
         return cls._eval_at_order(int(n), abs(int(m))).subs(_x, x)

    def _eval_expand_func(self, **hints):
        n, z = self.args

        if n.is_Integer and n.is_nonnegative:
            if z.is_Add:
                coeff = z.args[0]
                if coeff.is_Integer:
                    e = -(n + 1)
                    if coeff > 0:
                        tail = Add(*[Pow(
                            z - i, e) for i in range(1, int(coeff) + 1)])
                    else:
                        tail = -Add(*[Pow(
                            z + i, e) for i in range(0, int(-coeff))])
                    return polygamma(n, z - coeff) + (-1)**n*factorial(n)*tail

            elif z.is_Mul:
                coeff, z = z.as_two_terms()
                if coeff.is_Integer and coeff.is_positive:
                    tail = [ polygamma(n, z + Rational(
                        i, coeff)) for i in range(0, int(coeff)) ]
                    if n == 0:
                        return Add(*tail)/coeff + log(coeff)
                    else:
                        return Add(*tail)/coeff**(n + 1)
                z *= coeff

        return polygamma(n, z)

def _compute_formula(f, x, P, Q, k, m, k_max):
    """Computes the formula for f."""
    from sympy.polys import roots

    sol = []
    for i in range(k_max + 1, k_max + m + 1):
        r = f.diff(x, i).limit(x, 0) / factorial(i)
        if r is S.Zero:
            continue

        kterm = m*k + i
        res = r

        p = P.subs(k, kterm)
        q = Q.subs(k, kterm)
        c1 = p.subs(k, 1/k).leadterm(k)[0]
        c2 = q.subs(k, 1/k).leadterm(k)[0]
        res *= (-c1 / c2)**k

        for r, mul in roots(p, k).items():
            res *= rf(-r, k)**mul
        for r, mul in roots(q, k).items():
            res /= rf(-r, k)**mul

        sol.append((res, kterm))

    return sol

    def eval(cls, arg):
        if arg.is_Number:
            if arg is S.NaN:
                return S.NaN
            elif arg is S.Infinity:
                return S.Infinity
            elif arg.is_Integer:
                if arg.is_positive:
                    return factorial(arg - 1)
                else:
                    return S.ComplexInfinity
            elif arg.is_Rational:
                if arg.q == 2:
                    n = abs(arg.p) // arg.q

                    if arg.is_positive:
                        k, coeff = n, S.One
                    else:
                        n = k = n + 1

                        if n & 1 == 0:
                            coeff = S.One
                        else:
                            coeff = S.NegativeOne

                    for i in range(3, 2*k, 2):
                        coeff *= i

                    if arg.is_positive:
                        return coeff*sqrt(S.Pi) / 2**n
                    else:
                        return 2**n*sqrt(S.Pi) / coeff

        if arg.is_integer and arg.is_nonpositive:
            return S.ComplexInfinity

 def taylor_term(n, x, *previous_terms):
     from sympy.functions.combinatorial.numbers import euler
     if n < 0 or n % 2 == 1:
         return S.Zero
     else:
         x = sympify(x)
         return euler(n) / factorial(n) * x**(n)

def jacobi_normalized(n, a, b, x):
    r"""
    Jacobi polynomial :math:`P_n^{\left(\alpha, \beta\right)}(x)`

    jacobi_normalized(n, alpha, beta, x) gives the nth Jacobi polynomial
    in x, :math:`P_n^{\left(\alpha, \beta\right)}(x)`.

    The Jacobi polynomials are orthogonal on :math:`[-1, 1]` with respect
    to the weight :math:`\left(1-x\right)^\alpha \left(1+x\right)^\beta`.

    This functions returns the polynomials normilzed:

    .. math::

        \int_{-1}^{1}
          P_m^{\left(\alpha, \beta\right)}(x)
          P_n^{\left(\alpha, \beta\right)}(x)
          (1-x)^{\alpha} (1+x)^{\beta} \mathrm{d}x
        = \delta_{m,n}

    Examples
    ========

    >>> from sympy import jacobi_normalized
    >>> from sympy.abc import n,a,b,x

    >>> jacobi_normalized(n, a, b, x)
    jacobi(n, a, b, x)/sqrt(2**(a + b + 1)*gamma(a + n + 1)*gamma(b + n + 1)/((a + b + 2*n + 1)*factorial(n)*gamma(a + b + n + 1)))

    See Also
    ========

    gegenbauer,
    chebyshevt_root, chebyshevu, chebyshevu_root,
    legendre, assoc_legendre,
    hermite,
    laguerre, assoc_laguerre,
    sympy.polys.orthopolys.jacobi_poly,
    sympy.polys.orthopolys.gegenbauer_poly
    sympy.polys.orthopolys.chebyshevt_poly
    sympy.polys.orthopolys.chebyshevu_poly
    sympy.polys.orthopolys.hermite_poly
    sympy.polys.orthopolys.legendre_poly
    sympy.polys.orthopolys.laguerre_poly

    References
    ==========

    .. [1] http://en.wikipedia.org/wiki/Jacobi_polynomials
    .. [2] http://mathworld.wolfram.com/JacobiPolynomial.html
    .. [3] http://functions.wolfram.com/Polynomials/JacobiP/
    """
    nfactor = (
        S(2) ** (a + b + 1)
        * (gamma(n + a + 1) * gamma(n + b + 1))
        / (2 * n + a + b + 1)
        / (factorial(n) * gamma(n + a + b + 1))
    )

    return jacobi(n, a, b, x) / sqrt(nfactor)

    def eval(cls, a, x):
        # For lack of a better place, we use this one to extract branching
        # information. The following can be
        # found in the literature (c/f references given above), albeit scattered:
        # 1) For fixed x != 0, lowergamma(s, x) is an entire function of s
        # 2) For fixed positive integers s, lowergamma(s, x) is an entire
        #    function of x.
        # 3) For fixed non-positive integers s,
        #    lowergamma(s, exp(I*2*pi*n)*x) =
        #              2*pi*I*n*(-1)**(-s)/factorial(-s) + lowergamma(s, x)
        #    (this follows from lowergamma(s, x).diff(x) = x**(s-1)*exp(-x)).
        # 4) For fixed non-integral s,
        #    lowergamma(s, x) = x**s*gamma(s)*lowergamma_unbranched(s, x),
        #    where lowergamma_unbranched(s, x) is an entire function (in fact
        #    of both s and x), i.e.
        #    lowergamma(s, exp(2*I*pi*n)*x) = exp(2*pi*I*n*a)*lowergamma(a, x)
        from sympy import unpolarify, I
        if x == 0:
            return S.Zero
        nx, n = x.extract_branch_factor()
        if a.is_integer and a.is_positive:
            nx = unpolarify(x)
            if nx != x:
                return lowergamma(a, nx)
        elif a.is_integer and a.is_nonpositive:
            if n != 0:
                return 2*pi*I*n*(-1)**(-a)/factorial(-a) + lowergamma(a, nx)
        elif n != 0:
            return exp(2*pi*I*n*a)*lowergamma(a, nx)

        # Special values.
        if a.is_Number:
            if a is S.One:
                return S.One - exp(-x)
            elif a is S.Half:
                return sqrt(pi)*erf(sqrt(x))
            elif a.is_Integer or (2*a).is_Integer:
                b = a - 1
                if b.is_positive:
                    if a.is_integer:
                        return factorial(b) - exp(-x) * factorial(b) * Add(*[x ** k / factorial(k) for k in range(a)])
                    else:
                        return gamma(a) * (lowergamma(S.Half, x)/sqrt(pi) - exp(-x) * Add(*[x**(k-S.Half) / gamma(S.Half+k) for k in range(1, a+S.Half)]))

                if not a.is_Integer:
                    return (-1)**(S.Half - a) * pi*erf(sqrt(x)) / gamma(1-a) + exp(-x) * Add(*[x**(k+a-1)*gamma(a) / gamma(a+k) for k in range(1, S(3)/2-a)])

 def _eval_rewrite_as_polynomial(self, n, x, **kwargs):
     from sympy import Sum
     # Make sure n \in N_0
     if n.is_negative or n.is_integer is False:
         raise ValueError("Error: n should be a non-negative integer.")
     k = Dummy("k")
     kern = RisingFactorial(-n, k) / factorial(k)**2 * x**k
     return Sum(kern, (k, 0, n))

    def _eval_rewrite_as_Sum(self, n, k_sym=None, symbols=None):
        if (k_sym is not None) or (symbols is not None):
            return self

        # Dobinski's formula
        if not n.is_nonnegative:
            return self
        k = C.Dummy('k', integer=True, nonnegative=True)
        return 1 / E * C.Sum(k**n / factorial(k), (k, 0, S.Infinity))

    def f(rv):
        if not rv.is_Mul:
            return rv
        rvd = rv.as_powers_dict()
        nd_fact_args = [[], []] # numerator, denominator

        for k in rvd:
            if isinstance(k, factorial) and rvd[k].is_Integer:
                if rvd[k].is_positive:
                    nd_fact_args[0].extend([k.args[0]]*rvd[k])
                else:
                    nd_fact_args[1].extend([k.args[0]]*-rvd[k])
                rvd[k] = 0
        if not nd_fact_args[0] or not nd_fact_args[1]:
            return rv

        hit = False
        for m in range(2):
            i = 0
            while i < len(nd_fact_args[m]):
                ai = nd_fact_args[m][i]
                for j in range(i + 1, len(nd_fact_args[m])):
                    aj = nd_fact_args[m][j]

                    sum = ai + aj
                    if sum in nd_fact_args[1 - m]:
                        hit = True

                        nd_fact_args[1 - m].remove(sum)
                        del nd_fact_args[m][j]
                        del nd_fact_args[m][i]

                        rvd[binomial(sum, ai if count_ops(ai) <
                                count_ops(aj) else aj)] += (
                                -1 if m == 0 else 1)
                        break
                else:
                    i += 1

        if hit:
            return Mul(*([k**rvd[k] for k in rvd] + [factorial(k)
                    for k in nd_fact_args[0]]))/Mul(*[factorial(k)
                    for k in nd_fact_args[1]])
        return rv