def _eval_expand_func(self, **hints):
        from sympy import Sum
        n = self.args[0]
        m = self.args[1] if len(self.args) == 2 else 1

        if m == S.One:
            if n.is_Add:
                off = n.args[0]
                nnew = n - off
                if off.is_Integer and off.is_positive:
                    result = [S.One/(nnew + i) for i in range(off, 0, -1)] + [harmonic(nnew)]
                    return Add(*result)
                elif off.is_Integer and off.is_negative:
                    result = [-S.One/(nnew + i) for i in range(0, off, -1)] + [harmonic(nnew)]
                    return Add(*result)

            if n.is_Rational:
                # Expansions for harmonic numbers at general rational arguments (u + p/q)
                # Split n as u + p/q with p < q
                p, q = n.as_numer_denom()
                u = p // q
                p = p - u * q
                if u.is_nonnegative and p.is_positive and q.is_positive and p < q:
                    k = Dummy("k")
                    t1 = q * Sum(1 / (q * k + p), (k, 0, u))
                    t2 = 2 * Sum(cos((2 * pi * p * k) / S(q)) *
                                   log(sin((pi * k) / S(q))),
                                   (k, 1, floor((q - 1) / S(2))))
                    t3 = (pi / 2) * cot((pi * p) / q) + log(2 * q)
                    return t1 + t2 - t3

        return self

 def eval(cls, arg):
     if arg.is_Rational:
         return log(arg + S.One)
     elif not arg.is_Float:  # not safe to add 1 to Float
         return log.eval(arg + S.One)
     elif arg.is_number:
         return log(Rational(arg) + S.One)

def test_power_rewrite_exp():
    assert (I**I).rewrite(exp) == exp(-pi/2)

    expr = (2 + 3*I)**(4 + 5*I)
    assert expr.rewrite(exp) == exp((4 + 5*I)*(log(sqrt(13)) + I*atan(S(3)/2)))
    assert expr.rewrite(exp).expand() == \
        169*exp(5*I*log(13)/2)*exp(4*I*atan(S(3)/2))*exp(-5*atan(S(3)/2))

    assert ((6 + 7*I)**5).rewrite(exp) == 7225*sqrt(85)*exp(5*I*atan(S(7)/6))

    expr = 5**(6 + 7*I)
    assert expr.rewrite(exp) == exp((6 + 7*I)*log(5))
    assert expr.rewrite(exp).expand() == 15625*exp(7*I*log(5))

    assert Pow(123, 789, evaluate=False).rewrite(exp) == 123**789
    assert (1**I).rewrite(exp) == 1**I
    assert (0**I).rewrite(exp) == 0**I

    expr = (-2)**(2 + 5*I)
    assert expr.rewrite(exp) == exp((2 + 5*I)*(log(2) + I*pi))
    assert expr.rewrite(exp).expand() == 4*exp(-5*pi)*exp(5*I*log(2))

    assert ((-2)**S(-5)).rewrite(exp) == (-2)**S(-5)

    x, y = symbols('x y')
    assert (x**y).rewrite(exp) == exp(y*log(x))
    assert (7**x).rewrite(exp) == exp(x*log(7), evaluate=False)
    assert ((2 + 3*I)**x).rewrite(exp) == exp(x*(log(sqrt(13)) + I*atan(S(3)/2)))
    assert (y**(5 + 6*I)).rewrite(exp) == exp(log(y)*(5 + 6*I))

    assert all((1/func(x)).rewrite(exp) == 1/(func(x).rewrite(exp)) for func in
                    (sin, cos, tan, sec, csc, sinh, cosh, tanh))

    def eval(cls, arg):
        from sympy import asin
        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 log(sqrt(2) + 1)
            elif arg is S.NegativeOne:
                return log(sqrt(2) - 1)
            elif arg.is_negative:
                return -cls(-arg)
        else:
            if arg is S.ComplexInfinity:
                return S.ComplexInfinity

            i_coeff = arg.as_coefficient(S.ImaginaryUnit)

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

def get_math_macros():
    """ Returns a dictionary with math-related macros from math.h/cmath

    Note that these macros are not strictly required by the C/C++-standard.
    For MSVC they are enabled by defining "_USE_MATH_DEFINES" (preferably
    via a compilation flag).

    Returns
    =======

    Dictionary mapping sympy expressions to strings (macro names)

    """
    from sympy.codegen.cfunctions import log2, Sqrt
    from sympy.functions.elementary.exponential import log
    from sympy.functions.elementary.miscellaneous import sqrt

    return {
        S.Exp1: 'M_E',
        log2(S.Exp1): 'M_LOG2E',
        1/log(2): 'M_LOG2E',
        log(2): 'M_LN2',
        log(10): 'M_LN10',
        S.Pi: 'M_PI',
        S.Pi/2: 'M_PI_2',
        S.Pi/4: 'M_PI_4',
        1/S.Pi: 'M_1_PI',
        2/S.Pi: 'M_2_PI',
        2/sqrt(S.Pi): 'M_2_SQRTPI',
        2/Sqrt(S.Pi): 'M_2_SQRTPI',
        sqrt(2): 'M_SQRT2',
        Sqrt(2): 'M_SQRT2',
        1/sqrt(2): 'M_SQRT1_2',
        1/Sqrt(2): 'M_SQRT1_2'
    }

    def eval(cls, arg):
        from sympy import acos
        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.Infinity
            elif arg is S.Zero:
                return S.Pi*S.ImaginaryUnit / 2
            elif arg is S.One:
                return S.Zero
            elif arg is S.NegativeOne:
                return S.Pi*S.ImaginaryUnit

        if arg.is_number:
            cst_table = {
                S.ImaginaryUnit: log(S.ImaginaryUnit*(1 + sqrt(2))),
                -S.ImaginaryUnit: log(-S.ImaginaryUnit*(1 + sqrt(2))),
                S.Half: S.Pi/3,
                -S.Half: 2*S.Pi/3,
                sqrt(2)/2: S.Pi/4,
                -sqrt(2)/2: 3*S.Pi/4,
                1/sqrt(2): S.Pi/4,
                -1/sqrt(2): 3*S.Pi/4,
                sqrt(3)/2: S.Pi/6,
                -sqrt(3)/2: 5*S.Pi/6,
                (sqrt(3) - 1)/sqrt(2**3): 5*S.Pi/12,
                -(sqrt(3) - 1)/sqrt(2**3): 7*S.Pi/12,
                sqrt(2 + sqrt(2))/2: S.Pi/8,
                -sqrt(2 + sqrt(2))/2: 7*S.Pi/8,
                sqrt(2 - sqrt(2))/2: 3*S.Pi/8,
                -sqrt(2 - sqrt(2))/2: 5*S.Pi/8,
                (1 + sqrt(3))/(2*sqrt(2)): S.Pi/12,
                -(1 + sqrt(3))/(2*sqrt(2)): 11*S.Pi/12,
                (sqrt(5) + 1)/4: S.Pi/5,
                -(sqrt(5) + 1)/4: 4*S.Pi/5
            }

            if arg in cst_table:
                if arg.is_real:
                    return cst_table[arg]*S.ImaginaryUnit
                return cst_table[arg]

        if arg is S.ComplexInfinity:
            return S.Infinity

        i_coeff = arg.as_coefficient(S.ImaginaryUnit)

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

    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 test_trigintegrate_mixed():
    assert trigintegrate(sin(x)*sec(x), x) == -log(sin(x)**2 - 1)/2
    assert trigintegrate(sin(x)*csc(x), x) == x
    assert trigintegrate(sin(x)*cot(x), x) == sin(x)

    assert trigintegrate(cos(x)*sec(x), x) == x
    assert trigintegrate(cos(x)*csc(x), x) == log(cos(x)**2 - 1)/2
    assert trigintegrate(cos(x)*tan(x), x) == -cos(x)
    assert trigintegrate(cos(x)*cot(x), x) == log(cos(x) - 1)/2 \
        - log(cos(x) + 1)/2 + cos(x)

def composite(nth):
    """ Return the nth composite number, with the composite numbers indexed as
        composite(1) = 4, composite(2) = 6, etc....

        Examples
        ========

        >>> from sympy import composite
        >>> composite(36)
        52
        >>> composite(1)
        4
        >>> composite(17737)
        20000

        See Also
        ========

        sympy.ntheory.primetest.isprime : Test if n is prime
        primerange : Generate all primes in a given range
        primepi : Return the number of primes less than or equal to n
        prime : Return the nth prime
        compositepi : Return the number of positive composite numbers less than or equal to n
    """
    n = as_int(nth)
    if n < 1:
        raise ValueError("nth must be a positive integer; composite(1) == 4")
    composite_arr = [4, 6, 8, 9, 10, 12, 14, 15, 16, 18]
    if n <= 10:
        return composite_arr[n - 1]

    from sympy.functions.special.error_functions import li
    from sympy.functions.elementary.exponential import log

    a = 4 # Lower bound for binary search
    b = int(n*(log(n) + log(log(n)))) # Upper bound for the search.

    while a < b:
        mid = (a + b) >> 1
        if mid - li(mid) - 1 > n:
            b = mid
        else:
            a = mid + 1

    n_composites = a - primepi(a) - 1
    while n_composites > n:
        if not isprime(a):
            n_composites -= 1
        a -= 1
    if isprime(a):
        a -= 1
    return a

 def _eval_aseries(self, n, args0, x, logx):
     if args0[0] != oo:
         return super(loggamma, self)._eval_aseries(n, args0, x, logx)
     z = self.args[0]
     m = min(n, C.ceiling((n + S(1)) / 2))
     r = log(z) * (z - S(1) / 2) - z + log(2 * pi) / 2
     l = [bernoulli(2 * k) / (2 * k * (2 * k - 1) * z ** (2 * k - 1)) for k in range(1, m)]
     o = None
     if m == 0:
         o = C.Order(1, x)
     else:
         o = C.Order(1 / z ** (2 * m - 1), x)
     # It is very inefficient to first add the order and then do the nseries
     return (r + Add(*l))._eval_nseries(x, n, logx) + o

    def _eval_expand_func(self, deep=True, **hints):
        if deep:
            hints['func'] = False
            n = self.args[0].expand(deep, **hints)
            z = self.args[1].expand(deep, **hints)
        else:
            n, z = self.args[0], self.args[1].expand(deep, func=True)

        if n.is_Integer and n.is_nonnegative:
            if z.is_Add:
                coeff, factors = z.as_coeff_factors()

                if coeff.is_Integer:
                    tail = Add(*[ z + i for i in xrange(0, int(coeff)) ])
                    return polygamma(n, z-coeff) + (-1)**n*C.Factorial(n)*tail
            elif z.is_Mul:
                coeff, terms = z.as_coeff_terms()

                if coeff.is_Integer and coeff.is_positive:
                    tail = [ polygamma(n, z + i//coeff) for i in xrange(0, int(coeff)) ]

                    if n is S.Zero:
                        return log(coeff) + Add(*tail)/coeff**(n+1)
                    else:
                        return Add(*tail)/coeff**(n+1)

        return polygamma(n, z)

    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(*[C.Pow(z - i, e) for i in xrange(1, int(coeff) + 1)])
                    else:
                        tail = -Add(*[C.Pow(z + i, e) for i in xrange(0, int(-coeff))])
                    return polygamma(n, z - coeff) + (-1) ** n * C.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 + C.Rational(i, coeff)) for i in xrange(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 _eval_as_leading_term(self, x):
     n, z = [a.as_leading_term(x) for a in self.args]
     o = C.Order(z, x)
     if n == 0 and o.contains(1 / x):
         return o.getn() * log(x)
     else:
         return self.func(n, z)

def test_issue_7638():
    f = pi/log(sqrt(2))
    assert ((1 + I)**(I*f/2))**0.3 == (1 + I)**(0.15*I*f)
    # if 1/3 -> 1.0/3 this should fail since it cannot be shown that the
    # sign will be +/-1; for the previous "small arg" case, it didn't matter
    # that this could not be proved
    assert (1 + I)**(4*I*f) == ((1 + I)**(12*I*f))**(S(1)/3)

    assert (((1 + I)**(I*(1 + 7*f)))**(S(1)/3)).exp == S(1)/3
    r = symbols('r', real=True)
    assert sqrt(r**2) == abs(r)
    assert cbrt(r**3) != r
    assert sqrt(Pow(2*I, 5*S.Half)) != (2*I)**(5/S(4))
    p = symbols('p', positive=True)
    assert cbrt(p**2) == p**(2/S(3))
    assert NS(((0.2 + 0.7*I)**(0.7 + 1.0*I))**(0.5 - 0.1*I), 1) == '0.4 + 0.2*I'
    assert sqrt(1/(1 + I)) == sqrt(1 - I)/sqrt(2)  # or 1/sqrt(1 + I)
    e = 1/(1 - sqrt(2))
    assert sqrt(e) == I/sqrt(-1 + sqrt(2))
    assert e**-S.Half == -I*sqrt(-1 + sqrt(2))
    assert sqrt((cos(1)**2 + sin(1)**2 - 1)**(3 + I)).exp == S.Half
    assert sqrt(r**(4/S(3))) != r**(2/S(3))
    assert sqrt((p + I)**(4/S(3))) == (p + I)**(2/S(3))
    assert sqrt((p - p**2*I)**2) == p - p**2*I
    assert sqrt((p + r*I)**2) != p + r*I
    e = (1 + I/5)
    assert sqrt(e**5) == e**(5*S.Half)
    assert sqrt(e**6) == e**3
    assert sqrt((1 + I*r)**6) != (1 + I*r)**3

 def fdiff(self, argindex=1):
     """
     Returns the first derivative of this function.
     """
     if argindex == 1:
         return S.One/(log(_Ten)*self.args[0])
     else:
         raise ArgumentIndexError(self, argindex)

    def is_convergent(self):
        r"""
        See docs of Sum.is_convergent() for explanation of convergence
        in SymPy.

        The infinite product:

        .. math::

            \prod_{1 \leq i < \infty} f(i)

        is defined by the sequence of partial products:

        .. math::

            \prod_{i=1}^{n} f(i) = f(1) f(2) \cdots f(n)

        as n increases without bound. The product converges to a non-zero
        value if and only if the sum:

        .. math::

            \sum_{1 \leq i < \infty} \log{f(n)}

        converges.

        References
        ==========

        .. [1] https://en.wikipedia.org/wiki/Infinite_product

        Examples
        ========

        >>> from sympy import Interval, S, Product, Symbol, cos, pi, exp, oo
        >>> n = Symbol('n', integer=True)
        >>> Product(n/(n + 1), (n, 1, oo)).is_convergent()
        False
        >>> Product(1/n**2, (n, 1, oo)).is_convergent()
        False
        >>> Product(cos(pi/n), (n, 1, oo)).is_convergent()
        True
        >>> Product(exp(-n**2), (n, 1, oo)).is_convergent()
        False
        """
        from sympy.concrete.summations import Sum

        sequence_term = self.function
        log_sum = log(sequence_term)
        lim = self.limits
        try:
            is_conv = Sum(log_sum, *lim).is_convergent()
        except NotImplementedError:
            if Sum(sequence_term - 1, *lim).is_absolutely_convergent() is S.true:
                return S.true
            raise NotImplementedError("The algorithm to find the product convergence of %s "
                                        "is not yet implemented" % (sequence_term))
        return is_conv

    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

        if n == 0 and z.is_Rational:
            p, q = z.as_numer_denom()

            # Reference:
            #   Values of the polygamma functions at rational arguments, J. Choi, 2007
            part_1 = -S.EulerGamma - pi * cot(p * pi / q) / 2 - log(q) + Add(
                *[cos(2 * k * pi * p / q) * log(2 * sin(k * pi / q)) for k in range(1, q)])

            if z > 0:
                n = floor(z)
                z0 = z - n
                return part_1 + Add(*[1 / (z0 + k) for k in range(n)])
            elif z < 0:
                n = floor(1 - z)
                z0 = z + n
                return part_1 - Add(*[1 / (z0 - 1 - k) for k in range(n)])

        return polygamma(n, z)

 def fdiff(self, argindex=2):
     from sympy import meijerg
     if argindex == 2:
         a, z = self.args
         return -C.exp(-z)*z**(a-1)
     elif argindex == 1:
         a, z = self.args
         return uppergamma(a, z)*log(z) + meijerg([], [1, 1], [0, 0, a], [], z)
     else:
         raise ArgumentIndexError(self, argindex)

def test_issue_11463():
    numpy = import_module('numpy')
    if not numpy:
        skip("numpy not installed.")
    x = Symbol('x')
    f = lambdify(x, real_root((log(x/(x-2))), 3), 'numpy')
    # numpy.select evaluates all options before considering conditions,
    # so it raises a warning about root of negative number which does
    # not affect the outcome. This warning is suppressed here
    with ignore_warnings(RuntimeWarning):
        assert f(numpy.array(-1)) < -1

    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 == 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(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