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:
                coeff, terms = arg.as_coeff_terms()

                if coeff.is_negative:
                    return -cls(-arg)

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 _eval_expand_trig(self, **hints):
        from sympy import expand_mul
        arg = self.args[0]
        x = None
        if arg.is_Add:  # TODO, implement more if deep stuff here
            # TODO: Do this more efficiently for more than two terms
            x, y = arg.as_two_terms()
            sx = sin(x, evaluate=False)._eval_expand_trig()
            sy = sin(y, evaluate=False)._eval_expand_trig()
            cx = cos(x, evaluate=False)._eval_expand_trig()
            cy = cos(y, evaluate=False)._eval_expand_trig()
            return sx*cy + sy*cx
        else:
            n, x = arg.as_coeff_Mul(rational=True)
            if n.is_Integer:  # n will be positive because of .eval
                # canonicalization

                # See http://mathworld.wolfram.com/Multiple-AngleFormulas.html
                if n.is_odd:
                    return (-1)**((n - 1)/2)*C.chebyshevt(n, sin(x))
                else:
                    return expand_mul((-1)**(n/2 - 1)*cos(x)*C.chebyshevu(n -
                        1, sin(x)), deep=False)
            pi_coeff = _pi_coeff(arg)
            if pi_coeff is not None:
                if pi_coeff.is_Rational:
                    return self.rewrite(sqrt)
        return sin(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 _eval_expand_complex(self, *args):
     if self.args[0].is_real:
         return self
     re, im = self.args[0].as_real_imag()
     denom = sin(re)**2 + C.sinh(im)**2
     return (sin(re)*cos(re) - \
         S.ImaginaryUnit*C.sinh(im)*C.cosh(im))/denom

def solve_ODE_first_order(eq, f):
    """
    solves many kinds of first order odes, different methods are used
    depending on the form of the given equation. Now the linear
    and Bernoulli cases are implemented.
    """
    from sympy.integrals.integrals import integrate
    x = f.args[0]
    f = f.func

    #linear case: a(x)*f'(x)+b(x)*f(x)+c(x) = 0
    a = Wild('a', exclude=[f(x)])
    b = Wild('b', exclude=[f(x)])
    c = Wild('c', exclude=[f(x)])

    r = eq.match(a*diff(f(x),x) + b*f(x) + c)
    if r:
        t = C.exp(integrate(r[b]/r[a], x))
        tt = integrate(t*(-r[c]/r[a]), x)
        return (tt + Symbol("C1"))/t

    #Bernoulli case: a(x)*f'(x)+b(x)*f(x)+c(x)*f(x)^n = 0
    n = Wild('n', exclude=[f(x)])

    r = eq.match(a*diff(f(x),x) + b*f(x) + c*f(x)**n)
    if r:
        t = C.exp((1-r[n])*integrate(r[b]/r[a],x))
        tt = (r[n]-1)*integrate(t*r[c]/r[a],x)
        return ((tt + Symbol("C1"))/t)**(1/(1-r[n]))

    #other cases of first order odes will be implemented here

    raise NotImplementedError("solve_ODE_first_order: Cannot solve " + str(eq))

 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 _eval_expand_trig(self, **hints):
        arg = self.args[0]
        x = None
        if arg.is_Add:
            from sympy import symmetric_poly

            n = len(arg.args)
            CX = []
            for x in arg.args:
                cx = cot(x, evaluate=False)._eval_expand_trig()
                CX.append(cx)

            Yg = numbered_symbols("Y")
            Y = [Yg.next() for i in xrange(n)]

            p = [0, 0]
            for i in xrange(n, -1, -1):
                p[(n - i) % 2] += symmetric_poly(i, Y) * (-1) ** (((n - i) % 4) // 2)
            return (p[0] / p[1]).subs(zip(Y, CX))
        else:
            coeff, terms = arg.as_coeff_Mul(rational=True)
            if coeff.is_Integer and coeff > 1:
                I = S.ImaginaryUnit
                z = C.Symbol("dummy", real=True)
                P = ((z + I) ** coeff).expand()
                return (C.re(P) / C.im(P)).subs([(z, cot(terms))])
        return cot(arg)

 def _eval_rewrite_as_polynomial(self, n, m, x):
     k = C.Dummy("k")
     kern = (
         C.factorial(2 * n - 2 * k)
         / (2 ** n * C.factorial(n - k) * C.factorial(k) * C.factorial(n - 2 * k - m))
         * (-1) ** k
         * x ** (n - m - 2 * k)
     )
     return (1 - x ** 2) ** (m / 2) * C.Sum(kern, (k, 0, C.floor((n - m) * S.Half)))

 def _eval_expand_func(self, **hints):
     n, m, theta, phi = self.args
     rv = (
         sqrt((2 * n + 1) / (4 * pi) * C.factorial(n - m) / C.factorial(n + m))
         * C.exp(I * m * phi)
         * assoc_legendre(n, m, C.cos(theta))
     )
     # We can do this because of the range of theta
     return rv.subs(sqrt(-cos(theta) ** 2 + 1), sin(theta))

def solve_ODE_second_order(eq, f):
    """
    solves many kinds of second order odes, different methods are used
    depending on the form of the given equation. So far the constants
    coefficients case and a special case are implemented.
    """
    x = f.args[0]
    f = f.func

    #constant coefficients case: af''(x)+bf'(x)+cf(x)=0
    a = Wild('a', exclude=[x])
    b = Wild('b', exclude=[x])
    c = Wild('c', exclude=[x])

    r = eq.match(a*f(x).diff(x,x) + c*f(x))
    if r:
        return Symbol("C1")*C.sin(sqrt(r[c]/r[a])*x)+Symbol("C2")*C.cos(sqrt(r[c]/r[a])*x)

    r = eq.match(a*f(x).diff(x,x) + b*diff(f(x),x) + c*f(x))
    if r:
        r1 = solve(r[a]*x**2 + r[b]*x + r[c], x)
        if r1[0].is_real:
            if len(r1) == 1:
                return (Symbol("C1") + Symbol("C2")*x)*exp(r1[0]*x)
            else:
                return Symbol("C1")*exp(r1[0]*x) + Symbol("C2")*exp(r1[1]*x)
        else:
            r2 = abs((r1[0] - r1[1])/(2*S.ImaginaryUnit))
            return (Symbol("C2")*C.cos(r2*x) + Symbol("C1")*C.sin(r2*x))*exp((r1[0] + r1[1])*x/2)

    #other cases of the second order odes will be implemented here

    #special equations, that we know how to solve
    a = Wild('a')
    t = x*exp(f(x))
    tt = a*t.diff(x, x)/t
    r = eq.match(tt.expand())
    if r:
        return -solve_ODE_1(f(x), x)

    t = x*exp(-f(x))
    tt = a*t.diff(x, x)/t
    r = eq.match(tt.expand())
    if r:
        #check, that we've rewritten the equation correctly:
        #assert ( r[a]*t.diff(x,2)/t ) == eq.subs(f, t)
        return solve_ODE_1(f(x), x)

    neq = eq*exp(f(x))/exp(-f(x))
    r = neq.match(tt.expand())
    if r:
        #check, that we've rewritten the equation correctly:
        #assert ( t.diff(x,2)*r[a]/t ).expand() == eq
        return solve_ODE_1(f(x), x)

    raise NotImplementedError("solve_ODE_second_order: cannot solve " + str(eq))

    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))

 def _eval_expand_complex(self, deep=True, **hints):
     if self.args[0].is_real:
         if deep:
             hints['complex'] = False
             return self.expand(deep, **hints)
         else:
             return self
     if deep:
         re, im = self.args[0].expand(deep, **hints).as_real_imag()
     else:
         re, im = self.args[0].as_real_imag()
     return sin(re)*C.cosh(im) + S.ImaginaryUnit*cos(re)*C.sinh(im)

    def taylor_term(n, x, *previous_terms):
        if n == 0:
            return 1 / sympify(x)
        elif n < 0 or n % 2 == 0:
            return S.Zero
        else:
            x = sympify(x)

            B = C.bernoulli(n+1)
            F = C.factorial(n+1)

            return (-1)**((n+1)//2) * 2**(n+1) * B/F * x**n

 def as_real_imag(self, deep=True, **hints):
     if self.args[0].is_real:
         if deep:
             return (self.expand(deep, **hints), S.Zero)
         else:
             return (self, S.Zero)
     if deep:
         re, im = self.args[0].expand(deep, **hints).as_real_imag()
     else:
         re, im = self.args[0].as_real_imag()
     denom = sinh(re) ** 2 + C.sin(im) ** 2
     return (sinh(re) * cosh(re) / denom, -C.sin(im) * C.cos(im) / denom)

 def as_real_imag(self, deep=True, **hints):
     if self.args[0].is_real:
         if deep:
             hints['complex'] = False
             return (self.expand(deep, **hints), S.Zero)
         else:
             return (self, S.Zero)
     if deep:
         re, im = self.args[0].expand(deep, **hints).as_real_imag()
     else:
         re, im = self.args[0].as_real_imag()
     return (cos(re)*C.cosh(im), -sin(re)*C.sinh(im))

    def eval(cls, n, m, theta, phi):
        n, m, theta, phi = [sympify(x) for x in (n, m, theta, phi)]

        # Handle negative index m and arguments theta, phi
        if m.could_extract_minus_sign():
            m = -m
            return S.NegativeOne**m * C.exp(-2*I*m*phi) * Ynm(n, m, theta, phi)
        if theta.could_extract_minus_sign():
            theta = -theta
            return Ynm(n, m, theta, phi)
        if phi.could_extract_minus_sign():
            phi = -phi
            return C.exp(-2*I*m*phi) * Ynm(n, m, theta, phi)

 def _eval_expand_complex(self, deep=True, **hints):
     if self.args[0].is_real:
         if deep:
             return self.expand(deep, **hints)
         else:
             return self
     if deep:
         re, im = self.args[0].expand(deep, **hints).as_real_imag()
     else:
         re, im = self.args[0].as_real_imag()
     denom = sinh(re)**2 + C.sin(im)**2
     return (sinh(re)*cosh(re) - \
         S.ImaginaryUnit*C.sin(im)*C.cos(im))/denom

 def as_real_imag(self, deep=True, **hints):
     other = []
     coeff = S(1)
     for a in self.args:
         if a.is_real:
             coeff *= a
         else:
             other.append(a)
     m = Mul(*other)
     if hints.get('ignore') == m:
         return None
     else:
         return (coeff*C.re(m), coeff*C.im(m))

    def eval(cls, arg):
        if arg.is_integer:
            return arg
        if arg.is_imaginary or (S.ImaginaryUnit*arg).is_real:
            i = C.im(arg)
            if not i.has(S.ImaginaryUnit):
                return cls(i)*S.ImaginaryUnit
            return cls(arg, evaluate=False)

        v = cls._eval_number(arg)
        if v is not None:
            return v

        # Integral, numerical, symbolic part
        ipart = npart = spart = S.Zero

        # Extract integral (or complex integral) terms
        terms = Add.make_args(arg)

        for t in terms:
            if t.is_integer or (t.is_imaginary and C.im(t).is_integer):
                ipart += t
            elif t.has(C.Symbol):
                spart += t
            else:
                npart += t

        if not (npart or spart):
            return ipart

        # Evaluate npart numerically if independent of spart
        if npart and (
            not spart or
            npart.is_real and (spart.is_imaginary or (S.ImaginaryUnit*spart).is_real) or
                npart.is_imaginary and spart.is_real):
            try:
                re, im = get_integer_part(
                    npart, cls._dir, {}, return_ints=True)
                ipart += C.Integer(re) + C.Integer(im)*S.ImaginaryUnit
                npart = S.Zero
            except (PrecisionExhausted, NotImplementedError):
                pass

        spart += npart
        if not spart:
            return ipart
        elif spart.is_imaginary or (S.ImaginaryUnit*spart).is_real:
            return ipart + cls(C.im(spart), evaluate=False)*S.ImaginaryUnit
        else:
            return ipart + cls(spart, evaluate=False)