def monomial_count(V, N):
    r"""
    Computes the number of monomials.

    The number of monomials is given by the following formula:

    .. math::

        \frac{(\#V + N)!}{\#V! N!}

    where `N` is a total degree and `V` is a set of variables.

    **Examples**

    >>> 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 C.factorial(V + N) / C.factorial(V) / C.factorial(N)

    def arbitrary_point(self, parameter='t'):
        """A parameterized point on the ellipse.

        Parameters
        ----------
        parameter : str, optional
            Default value is 't'.

        Returns
        -------
        arbitrary_point : Point

        Raises
        ------
        ValueError
            When `parameter` already appears in the functions.

        See Also
        --------
        Point

        Examples
        --------
        >>> from sympy import Point, Ellipse
        >>> e1 = Ellipse(Point(0, 0), 3, 2)
        >>> e1.arbitrary_point()
        Point(3*cos(t), 2*sin(t))

        """
        t = _symbol(parameter)
        if t.name in (f.name for f in self.free_symbols):
            raise ValueError('Symbol %s already appears in object and cannot be used as a parameter.' % t.name)
        return Point(self.center[0] + self.hradius*C.cos(t),
                self.center[1] + self.vradius*C.sin(t))

def monomial_count(V, N):
    r"""
    Computes the number of monomials.

    The number of monomials is given by the following formula:

    .. math::

        \frac{(\#V + N)!}{\#V! N!}

    where `N` is a total degree and `V` is a set of variables.

    Examples
    ========

    >>> from sympy.polys.monomials import itermonomials, monomial_count
    >>> from sympy.polys.orderings import monomial_key
    >>> from sympy.abc import x, y

    >>> monomial_count(2, 2)
    6

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

    >>> sorted(M, key=monomial_key('grlex', [y, x]))
    [1, x, y, x**2, x*y, y**2]
    >>> len(M)
    6

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

    def vertices(self):
        """The vertices of the regular polygon.

        Returns
        -------
        vertices : list
            Each vertex is a Point.

        See Also
        --------
        Point

        Examples
        --------
        >>> from sympy.geometry import RegularPolygon, Point
        >>> rp = RegularPolygon(Point(0, 0), 5, 4)
        >>> rp.vertices
        [Point(5, 0), Point(0, 5), Point(-5, 0), Point(0, -5)]

        """
        points = []
        c, r, n = self
        v = 2*S.Pi/n
        for k in xrange(0, n):
            points.append(Point(c[0] + r*C.cos(k*v), c[1] + r*C.sin(k*v)))
        return points

    def eval(cls, z, a=S.One):
        z, a = map(sympify, (z, a))

        if a.is_Number:
            if a is S.NaN:
                return S.NaN
            elif a is S.Zero:
                return cls(z)

        if z.is_Number:
            if z is S.NaN:
                return S.NaN
            elif z is S.Infinity:
                return S.One
            elif z is S.Zero:
                if a.is_negative:
                    return S.Half - a - 1
                else:
                    return S.Half - a
            elif z is S.One:
                return S.ComplexInfinity
            elif z.is_Integer:
                if a.is_Integer:
                    if z.is_negative:
                        zeta = (-1)**z * C.bernoulli(-z+1)/(-z+1)
                    elif z.is_even:
                        B, F = C.bernoulli(z), C.factorial(z)
                        zeta = 2**(z-1) * abs(B) * pi**z / F
                    else:
                        return

                    if a.is_negative:
                        return zeta + C.harmonic(abs(a), z)
                    else:
                        return zeta - C.harmonic(a-1, z)

    def eval(cls, arg):
        arg = sympify(arg)

        if arg.is_Number:
            if arg is S.NaN:
                return S.NaN
            elif arg is S.Infinity:
                return S.Infinity
            elif arg is S.NegativeInfinity:
                return S.NegativeInfinity
            elif arg is S.Zero:
                return S.Zero
            elif arg is S.One:
                return C.log(2**S.Half + 1)
            elif arg is S.NegativeOne:
                return C.log(2**S.Half - 1)
            elif arg.is_negative:
                return -cls(-arg)
        else:
            i_coeff = arg.as_coefficient(S.ImaginaryUnit)

            if i_coeff is not None:
                return S.ImaginaryUnit * C.asin(i_coeff)
            else:
                if arg.as_coeff_mul()[0].is_negative:
                    return -cls(-arg)

 def vertices(self):
     points = []
     c, r, n = self
     v = 2*S.Pi/n
     for k in xrange(0, n):
         points.append( Point(c[0] + r*C.cos(k*v), c[1] + r*C.sin(k*v)) )
     return points

    def arbitrary_point(self, parameter_name='t'):
        """A parametric point on the ellipse.

        Parameters
        ----------
        parameter_name : str, optional
            Default value is 't'.

        Returns
        -------
        arbitrary_point : Point

        See Also
        --------
        Point

        Examples
        --------
        >>> from sympy import Point, Ellipse
        >>> e1 = Ellipse(Point(0, 0), 3, 2)
        >>> e1.arbitrary_point()
        Point(3*cos(t), 2*sin(t))

        """
        t = C.Symbol(parameter_name, real=True)
        return Point(self.center[0] + self.hradius*C.cos(t),
                self.center[1] + self.vradius*C.sin(t))

 def as_real_imag(self, deep=True, **hints):
     re, im = self.args[0].as_real_imag()
     if deep:
         re = re.expand(deep, **hints)
         im = im.expand(deep, **hints)
     cos, sin = C.cos(im), C.sin(im)
     return (exp(re) * cos, exp(re) * sin)

 def _eval_expand_complex(self, deep=True, **hints):
     re, im = self.args[0].as_real_imag()
     if deep:
         re = re.expand(deep, **hints)
         im = im.expand(deep, **hints)
     cos, sin = C.cos(im), C.sin(im)
     return exp(re) * cos + S.ImaginaryUnit * exp(re) * sin

    def as_real_imag(self, deep=True, **hints):
        """
        Returns this function as a 2-tuple representing a complex number.

        Examples
        ========

        >>> from sympy import I
        >>> from sympy.abc import x
        >>> from sympy.functions import exp
        >>> exp(x).as_real_imag()
        (exp(re(x))*cos(im(x)), exp(re(x))*sin(im(x)))
        >>> exp(1).as_real_imag()
        (E, 0)
        >>> exp(I).as_real_imag()
        (cos(1), sin(1))
        >>> exp(1+I).as_real_imag()
        (E*cos(1), E*sin(1))

        See Also
        ========

        sympy.functions.elementary.complexes.re
        sympy.functions.elementary.complexes.im
        """
        re, im = self.args[0].as_real_imag()
        if deep:
            re = re.expand(deep, **hints)
            im = im.expand(deep, **hints)
        cos, sin = C.cos(im), C.sin(im)
        return (exp(re)*cos, exp(re)*sin)

    def eval(cls, arg):
        arg = sympify(arg)

        if arg.is_Number:
            if arg is S.NaN:
                return S.NaN
            elif arg is S.Zero:
                return S.Zero
            elif arg is S.One:
                return S.Infinity
            elif arg is S.NegativeOne:
                return S.NegativeInfinity
            elif arg is S.Infinity:
                return -S.ImaginaryUnit * C.atan(arg)
            elif arg is S.NegativeInfinity:
                return S.ImaginaryUnit * C.atan(-arg)
            elif arg.is_negative:
                return -cls(-arg)
        else:
            if arg is S.ComplexInfinity:
                return S.NaN

            i_coeff = arg.as_coefficient(S.ImaginaryUnit)

            if i_coeff is not None:
                return S.ImaginaryUnit * C.atan(i_coeff)
            else:
                if _coeff_isneg(arg):
                    return -cls(-arg)

def real_root(arg, n=None):
    """Return the real nth-root of arg if possible. If n is omitted then
    all instances of (-n)**(1/odd) will be changed to -n**(1/odd); this
    will only create a real root of a principle root -- the presence of
    other factors may cause the result to not be real.

    Examples
    ========

    >>> from sympy import root, real_root, Rational
    >>> from sympy.abc import x, n

    >>> real_root(-8, 3)
    -2
    >>> root(-8, 3)
    2*(-1)**(1/3)
    >>> real_root(_)
    -2

    If one creates a non-principle root and applies real_root, the
    result will not be real (so use with caution):

    >>> root(-8, 3, 2)
    -2*(-1)**(2/3)
    >>> real_root(_)
    -2*(-1)**(2/3)


    See Also
    ========

    sympy.polys.rootoftools.RootOf
    sympy.core.power.integer_nthroot
    root, sqrt
    """
    if n is not None:
        try:
            n = as_int(n)
            arg = sympify(arg)
            if arg.is_positive or arg.is_negative:
                rv = root(arg, n)
            else:
                raise ValueError
        except ValueError:
            return root(arg, n)*C.Piecewise(
                (S.One, ~C.Equality(C.im(arg), 0)),
                (C.Pow(S.NegativeOne, S.One/n)**(2*C.floor(n/2)), C.And(
                    C.Equality(n % 2, 1),
                    arg < 0)),
                (S.One, True))
    else:
        rv = sympify(arg)
    n1pow = Transform(lambda x: -(-x.base)**x.exp,
                      lambda x:
                      x.is_Pow and
                      x.base.is_negative and
                      x.exp.is_Rational and
                      x.exp.p == 1 and x.exp.q % 2)
    return rv.xreplace(n1pow)

 def eval(cls, arg):
     from sympy.simplify.simplify import signsimp
     if hasattr(arg, '_eval_Abs'):
         obj = arg._eval_Abs()
         if obj is not None:
             return obj
     # handle what we can
     arg = signsimp(arg, evaluate=False)
     if arg.is_Mul:
         known = []
         unk = []
         for t in arg.args:
             tnew = cls(t)
             if tnew.func is cls:
                 unk.append(tnew.args[0])
             else:
                 known.append(tnew)
         known = Mul(*known)
         unk = cls(Mul(*unk), evaluate=False) if unk else S.One
         return known*unk
     if arg is S.NaN:
         return S.NaN
     if arg.is_Pow:
         base, exponent = arg.as_base_exp()
         if base.is_real:
             if exponent.is_integer:
                 if exponent.is_even:
                     return arg
                 if base is S.NegativeOne:
                     return S.One
                 if base.func is cls and exponent is S.NegativeOne:
                     return arg
                 return Abs(base)**exponent
             if base.is_positive == True:
                 return base**re(exponent)
             return (-base)**re(exponent)*C.exp(-S.Pi*im(exponent))
     if isinstance(arg, C.exp):
         return C.exp(re(arg.args[0]))
     if arg.is_zero:  # it may be an Expr that is zero
         return S.Zero
     if arg.is_nonnegative:
         return arg
     if arg.is_nonpositive:
         return -arg
     if arg.is_imaginary:
         arg2 = -S.ImaginaryUnit * arg
         if arg2.is_nonnegative:
             return arg2
     if arg.is_Add:
         if arg.has(S.Infinity, S.NegativeInfinity):
             if any(a.is_infinite for a in arg.as_real_imag()):
                 return S.Infinity
         if arg.is_real is None and arg.is_imaginary is None:
             if all(a.is_real or a.is_imaginary or (S.ImaginaryUnit*a).is_real for a in arg.args):
                 from sympy import expand_mul
                 return sqrt(expand_mul(arg*arg.conjugate()))
     if arg.is_real is False and arg.is_imaginary is False:
         from sympy import expand_mul
         return sqrt(expand_mul(arg*arg.conjugate()))

    def eval(cls, n, z):
        n, z = list(map(sympify, (n, z)))
        from sympy import unpolarify

        if n.is_integer:
            if n.is_nonnegative:
                nz = unpolarify(z)
                if z != nz:
                    return polygamma(n, nz)

            if n == -1:
                return loggamma(z)
            else:
                if z.is_Number:
                    if z is S.NaN:
                        return S.NaN
                    elif z is S.Infinity:
                        if n.is_Number:
                            if n is S.Zero:
                                return S.Infinity
                            else:
                                return S.Zero
                    elif z.is_Integer:
                        if z.is_nonpositive:
                            return S.ComplexInfinity
                        else:
                            if n is S.Zero:
                                return -S.EulerGamma + C.harmonic(z - 1, 1)
                            elif n.is_odd:
                                return (-1) ** (n + 1) * C.factorial(n) * zeta(n + 1, z)

        if n == 0:
            if z is S.NaN:
                return S.NaN
            elif z.is_Rational:
                # TODO actually *any* n/m can be done, but that is messy
                lookup = {
                    S(1) / 2: -2 * log(2) - S.EulerGamma,
                    S(1) / 3: -S.Pi / 2 / sqrt(3) - 3 * log(3) / 2 - S.EulerGamma,
                    S(1) / 4: -S.Pi / 2 - 3 * log(2) - S.EulerGamma,
                    S(3) / 4: -3 * log(2) - S.EulerGamma + S.Pi / 2,
                    S(2) / 3: -3 * log(3) / 2 + S.Pi / 2 / sqrt(3) - S.EulerGamma,
                }
                if z > 0:
                    n = floor(z)
                    z0 = z - n
                    if z0 in lookup:
                        return lookup[z0] + Add(*[1 / (z0 + k) for k in range(n)])
                elif z < 0:
                    n = floor(1 - z)
                    z0 = z + n
                    if z0 in lookup:
                        return lookup[z0] - Add(*[1 / (z0 - 1 - k) for k in range(n)])
            elif z in (S.Infinity, S.NegativeInfinity):
                return S.Infinity
            else:
                t = z.extract_multiplicatively(S.ImaginaryUnit)
                if t in (S.Infinity, S.NegativeInfinity):
                    return S.Infinity

    def _eval_aseries(self, n, args0, x, logx):
        if args0[0] != S.Infinity:
            return super(_erfs, self)._eval_aseries(n, args0, x, logx)

        z = self.args[0]
        l = [ 1/sqrt(S.Pi) * C.factorial(2*k)*(-S(4))**(-k)/C.factorial(k) * (1/z)**(2*k+1) for k in xrange(0,n) ]
        o = C.Order(1/z**(2*n+1), x)
        # It is very inefficient to first add the order and then do the nseries
        return (Add(*l))._eval_nseries(x, n, logx) + o

    def eval(cls, a, x):
        if a.is_Number:
            if a is S.One:
                return S.One - C.exp(-x)
            elif a.is_Integer:
                b = a - 1

                if b.is_positive:
                    return b*cls(b, x) - x**b * C.exp(-x)

    def eval(cls, n):
        if (n.is_Integer and n.is_nonnegative) or \
           (n.is_noninteger and n.is_negative):
            return 4**n*C.gamma(n + S.Half)/(C.gamma(S.Half)*C.gamma(n + 2))

        if (n.is_integer and n.is_negative):
            if (n + 1).is_negative:
                return S.Zero
            if (n + 1).is_zero:
                return -S.Half

 def taylor_term(n, x, *previous_terms):
     if n < 0 or n % 2 == 0:
         return S.Zero
     else:
         x = sympify(x)
         k = C.floor((n - 1)/S(2))
         if len(previous_terms) > 2:
             return -previous_terms[-2] * x**2 * (n - 2)/(n*k)
         else:
             return 2*(-1)**k * x**n/(n*C.factorial(k)*sqrt(S.Pi))

    def eval(cls, r, k):
        r, k = map(sympify, (r, k))

        if k.is_Number:
            if k is S.Zero:
                return S.One
            elif k.is_Integer:
                if k.is_negative:
                    return S.Zero
                else:
                    if r.is_Integer and r.is_nonnegative:
                        r, k = int(r), int(k)

                        if k > r:
                            return S.Zero
                        elif k > r // 2:
                            k = r - k

                        M, result = int(sqrt(r)), 1

                        for prime in sieve.primerange(2, r + 1):
                            if prime > r - k:
                                result *= prime
                            elif prime > r // 2:
                                continue
                            elif prime > M:
                                if r % prime < k % prime:
                                    result *= prime
                            else:
                                R, K = r, k
                                exp = a = 0

                                while R > 0:
                                    a = int((R % prime) < (K % prime + a))
                                    R, K = R // prime, K // prime
                                    exp = a + exp

                                if exp > 0:
                                    result *= prime ** exp

                        return C.Integer(result)
                    else:
                        result = r - k + 1

                        for i in xrange(2, k + 1):
                            result *= r - k + i
                            result /= i

                        return result

        if k.is_integer:
            if k.is_negative:
                return S.Zero
        else:
            return C.gamma(r + 1) / (C.gamma(r - k + 1) * C.gamma(k + 1))