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 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_rewrite_as_polynomial(self, n, x):
     from sympy import Sum
     # TODO: Should make sure n is in N_0
     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 _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 pdf(self, *args):
     from sympy.functions.special.gamma_functions import gamma
     mu, sigma = self.mu, self.shape_mat
     v = S(self.dof)
     k = S(len(mu))
     sigma_inv = sigma.inv()
     args = ImmutableMatrix(args)
     x = args - mu
     return gamma((k + v)/2)/(gamma(v/2)*(v*pi)**(k/2)*sqrt(det(sigma)))\
     *(1 + 1/v*(x.transpose()*sigma_inv*x)[0])**((-v - k)/2)

 def eval(cls, z):
     if z is S.Zero:
         return pi / 2
     elif z is S.Half:
         return 8 * pi ** (S(3) / 2) / gamma(-S(1) / 4) ** 2
     elif z is S.One:
         return S.ComplexInfinity
     elif z is S.NegativeOne:
         return gamma(S(1) / 4) ** 2 / (4 * sqrt(2 * pi))
     elif z in (S.Infinity, S.NegativeInfinity, I * S.Infinity, I * S.NegativeInfinity, S.ComplexInfinity):
         return S.Zero

 def marginal_distribution(self, indices, *sym):
     from sympy.functions.special.gamma_functions import gamma
     if len(indices) == 2:
         return self.pdf(*sym)
     if indices[0] == 0:
         #For marginal over `x`, return non-standardized Student-T's
         #distribution
         x = sym[0]
         v, mu, sigma = self.alpha - S(1)/2, self.mu, \
             S(self.beta)/(self.lamda * self.alpha)
         return Lambda(sym, gamma((v + 1)/2)/(gamma(v/2)*sqrt(pi*v)*sigma)*\
             (1 + 1/v*((x - mu)/sigma)**2)**((-v -1)/2))
     #For marginal over `tau`, return Gamma distribution as per construction
     from sympy.stats.crv_types import GammaDistribution
     return Lambda(sym, GammaDistribution(self.alpha, self.beta)(sym[0]))

 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(cls, arg):
     if arg.is_Number:
         if arg is S.NaN:
             return S.NaN
         elif arg is S.Infinity:
             return S.Zero
         elif arg is S.Zero:
             return -S.One / (3 ** Rational(1, 3) * gamma(Rational(1, 3)))

 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 eval(cls, n, x):
     if not n.is_Number:
         # Symbolic result L_n(x)
         # L_n(-x)  --->  (-1)**n * L_n(x)
         if x.could_extract_minus_sign():
             return S.NegativeOne**n * legendre(n, -x)
         # L_{-n}(x)  --->  L_{n-1}(x)
         if n.could_extract_minus_sign():
             return legendre(-n - S.One, x)
         # We can evaluate for some special values of x
         if x == S.Zero:
             return sqrt(S.Pi)/(gamma(S.Half - n/2)*gamma(S.One + n/2))
         elif x == S.One:
             return S.One
         elif x == S.Infinity:
             return S.Infinity
     else:
         # n is a given fixed integer, evaluate into polynomial;
         # L_{-n}(x)  --->  L_{n-1}(x)
         if n.is_negative:
             n = -n - S.One
         return cls._eval_at_order(n, x)

    def eval(cls, n, a, b, x):
        # Simplify to other polynomials
        # P^{a, a}_n(x)
        if a == b:
            if a == -S.Half:
                return RisingFactorial(S.Half, n) / factorial(n) * chebyshevt(n, x)
            elif a == S.Zero:
                return legendre(n, x)
            elif a == S.Half:
                return RisingFactorial(3*S.Half, n) / factorial(n + 1) * chebyshevu(n, x)
            else:
                return RisingFactorial(a + 1, n) / RisingFactorial(2*a + 1, n) * gegenbauer(n, a + S.Half, x)
        elif b == -a:
            # P^{a, -a}_n(x)
            return gamma(n + a + 1) / gamma(n + 1) * (1 + x)**(a/2) / (1 - x)**(a/2) * assoc_legendre(n, -a, x)


        if not n.is_Number:
            # Symbolic result P^{a,b}_n(x)
            # P^{a,b}_n(-x)  --->  (-1)**n * P^{b,a}_n(-x)
            if x.could_extract_minus_sign():
                return S.NegativeOne**n * jacobi(n, b, a, -x)
            # We can evaluate for some special values of x
            if x == S.Zero:
                return (2**(-n) * gamma(a + n + 1) / (gamma(a + 1) * factorial(n)) *
                        hyper([-b - n, -n], [a + 1], -1))
            if x == S.One:
                return RisingFactorial(a + 1, n) / factorial(n)
            elif x == S.Infinity:
                if n.is_positive:
                    # Make sure a+b+2*n \notin Z
                    if (a + b + 2*n).is_integer:
                        raise ValueError("Error. a + b + 2*n should not be an integer.")
                    return RisingFactorial(a + b + n + 1, n) * S.Infinity
        else:
            # n is a given fixed integer, evaluate into polynomial
            return jacobi_poly(n, a, b, x)

def gauss_gen_laguerre(n, alpha, n_digits):
    r"""
    Computes the generalized Gauss-Laguerre quadrature [1]_ points and weights.

    The generalized Gauss-Laguerre quadrature approximates the integral:

    .. math::
        \int_{0}^\infty x^{\alpha} e^{-x} f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)

    The nodes `x_i` of an order `n` quadrature rule are the roots of `L^{\alpha}_n`
    and the weights `w_i` are given by:

    .. math::
        w_i = \frac{\Gamma(\alpha+n)}{n \Gamma(n) L^{\alpha}_{n-1}(x_i) L^{\alpha+1}_{n-1}(x_i)}

    Parameters
    ==========

    n : the order of quadrature

    alpha : the exponent of the singularity, `\alpha > -1`

    n_digits : number of significant digits of the points and weights to return

    Returns
    =======

    (x, w) : the ``x`` and ``w`` are lists of points and weights as Floats.
             The points `x_i` and weights `w_i` are returned as ``(x, w)``
             tuple of lists.

    Examples
    ========

    >>> from sympy import S
    >>> from sympy.integrals.quadrature import gauss_gen_laguerre
    >>> x, w = gauss_gen_laguerre(3, -S.Half, 5)
    >>> x
    [0.19016, 1.7845, 5.5253]
    >>> w
    [1.4493, 0.31413, 0.00906]

    >>> x, w = gauss_gen_laguerre(4, 3*S.Half, 5)
    >>> x
    [0.97851, 2.9904, 6.3193, 11.712]
    >>> w
    [0.53087, 0.67721, 0.11895, 0.0023152]

    See Also
    ========

    gauss_legendre, gauss_laguerre, gauss_hermite, gauss_chebyshev_t, gauss_chebyshev_u, gauss_jacobi

    References
    ==========

    .. [1] http://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature
    .. [2] http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_laguerre_rule/gen_laguerre_rule.html
    """
    x = Dummy("x")
    p = laguerre_poly(n, x, alpha=alpha, polys=True)
    p1 = laguerre_poly(n-1, x, alpha=alpha, polys=True)
    p2 = laguerre_poly(n-1, x, alpha=alpha+1, polys=True)
    xi = []
    w  = []
    for r in p.real_roots():
        if isinstance(r, RootOf):
            r = r.eval_rational(S(1)/10**(n_digits+2))
        xi.append(r.n(n_digits))
        w.append((gamma(alpha+n)/(n*gamma(n)*p1.subs(x, r)*p2.subs(x, r))).n(n_digits))
    return xi, w

def gauss_jacobi(n, alpha, beta, n_digits):
    r"""
    Computes the Gauss-Jacobi quadrature [1]_ points and weights.

    The Gauss-Jacobi quadrature of the first kind approximates the integral:

    .. math::
        \int_{-1}^1 (1-x)^\alpha (1+x)^\beta f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)

    The nodes `x_i` of an order `n` quadrature rule are the roots of `P^{(\alpha,\beta)}_n`
    and the weights `w_i` are given by:

    .. math::
        w_i = -\frac{2n+\alpha+\beta+2}{n+\alpha+\beta+1}\frac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}
              {\Gamma(n+\alpha+\beta+1)(n+1)!} \frac{2^{\alpha+\beta}}{P'_n(x_i)
              P^{(\alpha,\beta)}_{n+1}(x_i)}

    Parameters
    ==========

    n : the order of quadrature

    alpha : the first parameter of the Jacobi Polynomial, `\alpha > -1`

    beta : the second parameter of the Jacobi Polynomial, `\beta > -1`

    n_digits : number of significant digits of the points and weights to return

    Returns
    =======

    (x, w) : the ``x`` and ``w`` are lists of points and weights as Floats.
             The points `x_i` and weights `w_i` are returned as ``(x, w)``
             tuple of lists.

    Examples
    ========

    >>> from sympy import S
    >>> from sympy.integrals.quadrature import gauss_jacobi
    >>> x, w = gauss_jacobi(3, S.Half, -S.Half, 5)
    >>> x
    [-0.90097, -0.22252, 0.62349]
    >>> w
    [1.7063, 1.0973, 0.33795]

    >>> x, w = gauss_jacobi(6, 1, 1, 5)
    >>> x
    [-0.87174, -0.5917, -0.2093, 0.2093, 0.5917, 0.87174]
    >>> w
    [0.050584, 0.22169, 0.39439, 0.39439, 0.22169, 0.050584]

    See Also
    ========

    gauss_legendre, gauss_laguerre, gauss_hermite, gauss_gen_laguerre, gauss_chebyshev_t, gauss_chebyshev_u

    References
    ==========

    .. [1] http://en.wikipedia.org/wiki/Gauss%E2%80%93Jacobi_quadrature
    .. [2] http://people.sc.fsu.edu/~jburkardt/cpp_src/jacobi_rule/jacobi_rule.html
    .. [3] http://people.sc.fsu.edu/~jburkardt/cpp_src/gegenbauer_rule/gegenbauer_rule.html
    """
    x = Dummy("x")
    p = jacobi_poly(n, alpha, beta, x, polys=True)
    pd = p.diff(x)
    pn = jacobi_poly(n+1, alpha, beta, x, polys=True)
    xi = []
    w  = []
    for r in p.real_roots():
        if isinstance(r, RootOf):
            r = r.eval_rational(S(1)/10**(n_digits+2))
        xi.append(r.n(n_digits))
        w.append((
            - (2*n+alpha+beta+2) / (n+alpha+beta+S.One)
            * (gamma(n+alpha+1)*gamma(n+beta+1)) / (gamma(n+alpha+beta+S.One)*gamma(n+2))
            * 2**(alpha+beta) / (pd.subs(x, r) * pn.subs(x, r))
        ).n(n_digits))
    return xi, w

 def _eval_rewrite_as_hyper(self, z):
     pf1 = z ** 2 / (2 * root(3, 6) * gamma(S(2) / 3))
     pf2 = root(3, 6) / gamma(S(1) / 3)
     return pf1 * hyper([], [S(5) / 3], z ** 3 / 9) + pf2 * hyper([], [S(1) / 3], z ** 3 / 9)

def test_sympy__functions__special__gamma_functions__gamma():
    from sympy.functions.special.gamma_functions import gamma
    assert _test_args(gamma(x))

 def _eval_expand_func(self, **hints):
     x, y = self.args
     return gamma(x)*gamma(y) / gamma(x + y)