def test_functions_basic():
    assert oclcode(abs(i)) == "abs(i)"
    assert oclcode(abs(x)) == "fabs(x)"
    assert oclcode(log(x)) == "log(x)"
    assert oclcode(log(x, 2)) == "log2(x)"
    assert oclcode(log(x, 10)) == "log10(x)"
    assert oclcode(log(x, 3)) == "log(x)/log(3)"

def test_2():
    test = [
        [(a + b*x)/x, x, S(2), b*x + a*log(x)],
        [x**S(5)/(a + b*x), x, S(2), a**S(4)*x/b**S(5) - S(1)/S(2)*a**S(3)*x**S(2)/b**S(4) + S(1)/S(3)*a**S(2)*x**S(3)/b**S(3) - S(1)/S(4)*a*x**S(4)/b**S(2) + S(1)/S(5)*x**S(5)/b - a**S(5)*log(a + b*x)/b**S(6)],
        [S(1)/(a + b*x)**S(2), x, S(1), ( - S(1))/(b*(a + b*x))],
        [S(1)/(x*(a + b*x)**S(3)), x, S(2), S(1)/S(2)/(a*(a + b*x)**S(2)) + S(1)/(a**S(2)*(a + b*x)) + log(x)/a**S(3) - log(a + b*x)/a**S(3)],
        [S(1)/(S(2) + S(2)*x), x, S(1), S(1)/S(2)*log(S(1) + x)],
        [S(1)/(x*(S(1) + b*x)), x, S(3), log(x) - log(S(1) + b*x)],
        [x**S(3)*sqrt(a + b*x), x, S(2), - S(2)/S(3)*a**S(3)*(a + b*x)**(S(3)/S(2))/b**S(4) + S(6)/S(5)*a**S(2)*(a + b*x)**(S(5)/S(2))/b**S(4) - S(6)/S(7)*a*(a + b*x)**(S(7)/S(2))/b**S(4) + S(2)/S(9)*(a + b*x)**(S(9)/S(2))/b**S(4)],
        [(a + b*x)**(S(3)/S(2)), x, S(1), S(2)/S(5)*(a + b*x)**(S(5)/S(2))/b],
        [x**S(4)/sqrt(a + b*x), x, S(2), - S(8)/S(3)*a**S(3)*(a + b*x)**(S(3)/S(2))/b**S(5) + S(12)/S(5)*a**S(2)*(a + b*x)**(S(5)/S(2))/b**S(5) - S(8)/S(7)*a*(a + b*x)**(S(7)/S(2))/b**S(5) + S(2)/S(9)*(a + b*x)**(S(9)/S(2))/b**S(5) + S(2)*a**S(4)*sqrt(a + b*x)/b**S(5)],
        [S(1)/sqrt(a + b*x), x, S(1), S(2)*sqrt(a + b*x)/b],
        [S(1)/(x*(a + b*x)**(S(3)/S(2))), x, S(3), - S(2)*arctanh(sqrt(a + b*x)/sqrt(a))/a**(S(3)/S(2)) + S(2)/(a*sqrt(a + b*x))],
        [S(1)/(x**S(2)*( - a + b*x)**(S(3)/S(2))), x, S(4), - S(3)*b*arctan(sqrt( - a + b*x)/sqrt(a))/a**(S(5)/S(2)) + ( - S(2))/(a*x*sqrt( - a + b*x)) - S(3)*sqrt( - a + b*x)/(a**S(2)*x)],
        [x**S(3)*(a + b*x)**(S(1)/S(3)), x, S(2), - S(3)/S(4)*a**S(3)*(a + b*x)**(S(4)/S(3))/b**S(4) + S(9)/S(7)*a**S(2)*(a + b*x)**(S(7)/S(3))/b**S(4) - S(9)/S(10)*a*(a + b*x)**(S(10)/S(3))/b**S(4) + S(3)/S(13)*(a + b*x)**(S(13)/S(3))/b**S(4)],
        [x**S(2)*(a + b*x)**(S(2)/S(3)), x, S(2), S(3)/S(5)*a**S(2)*(a + b*x)**(S(5)/S(3))/b**S(3) - S(3)/S(4)*a*(a + b*x)**(S(8)/S(3))/b**S(3) + S(3)/S(11)*(a + b*x)**(S(11)/S(3))/b**S(3)],
        [x**S(2)/(a + b*x)**(S(1)/S(3)), x, S(2), S(3)/S(2)*a**S(2)*(a + b*x)**(S(2)/S(3))/b**S(3) - S(6)/S(5)*a*(a + b*x)**(S(5)/S(3))/b**S(3) + S(3)/S(8)*(a + b*x)**(S(8)/S(3))/b**S(3)],
        [x**S(3)/( - a + b*x)**(S(1)/S(3)), x, S(2), S(3)/S(2)*a**S(3)*( - a + b*x)**(S(2)/S(3))/b**S(4) + S(9)/S(5)*a**S(2)*( - a + b*x)**(S(5)/S(3))/b**S(4) + S(9)/S(8)*a*( - a + b*x)**(S(8)/S(3))/b**S(4) + S(3)/S(11)*( - a + b*x)**(S(11)/S(3))/b**S(4)],
    ]

    for i in test:
        r = rubi_integrate(i[0], i[1])
        if len(i) == 5:
            assert rubi_test(r, i[1], i[3], expand=True, _diff=True) or rubi_test(r, i[1], i[4], expand=True, _diff=True)
        else:
            assert rubi_test(r, i[1], i[3], expand=True, _diff=True)

def _generate_patterns():
    """
    Generates patterns for transcendental equations.

    This is lazily calculated (called) in the tsolve() function and stored in
    the patterns global variable.
    """

    tmp1 = _f ** (_h-(_c*_g/_b))
    tmp2 = (-_e*tmp1/_a)**(1/_d)
    global _patterns
    _patterns = [
        (_a*(_b*_x+_c)**_d + _e   ,
            ((-(_e/_a))**(1/_d)-_c)/_b),
        (_b+_c*exp(_d*_x+_e) ,
            (log(-_b/_c)-_e)/_d),
        (_a*_x+_b+_c*exp(_d*_x+_e) ,
            -_b/_a-LambertW(_c*_d*exp(_e-_b*_d/_a)/_a)/_d),
        (_b+_c*_f**(_d*_x+_e) ,
            (log(-_b/_c)-_e*log(_f))/_d/log(_f)),
        (_a*_x+_b+_c*_f**(_d*_x+_e) ,
            -_b/_a-LambertW(_c*_d*_f**(_e-_b*_d/_a)*log(_f)/_a)/_d/log(_f)),
        (_b+_c*log(_d*_x+_e) ,
            (exp(-_b/_c)-_e)/_d),
        (_a*_x+_b+_c*log(_d*_x+_e) ,
            -_e/_d+_c/_a*LambertW(_a/_c/_d*exp(-_b/_c+_a*_e/_c/_d))),
        (_a*(_b*_x+_c)**_d + _e*_f**(_g*_x+_h) ,
            -_c/_b-_d*LambertW(-tmp2*_g*log(_f)/_b/_d)/_g/log(_f))
    ]

def test_log():
    R, x = ring('x', QQ)
    p = 1 + x
    p1 = rs_log(p, x, 4)/x**2
    assert p1 == S(1)/3*x - S(1)/2 + x**(-1)
    p = 1 + x +2*x**2/3
    p1 = rs_log(p, x, 9)
    assert p1 == -17*x**8/648 + 13*x**7/189 - 11*x**6/162 - x**5/45 + \
      7*x**4/36 - x**3/3 + x**2/6 + x
    p2 = rs_series_inversion(p, x, 9)
    p3 = rs_log(p2, x, 9)
    assert p3 == -p1

    R, x, y = ring('x, y', QQ)
    p = 1 + x + 2*y*x**2
    p1 = rs_log(p, x, 6)
    assert p1 == (4*x**5*y**2 - 2*x**5*y - 2*x**4*y**2 + x**5/5 + 2*x**4*y -
                  x**4/4 - 2*x**3*y + x**3/3 + 2*x**2*y - x**2/2 + x)

    # Constant term in series
    a = symbols('a')
    R, x, y = ring('x, y', EX)
    assert rs_log(x + a, x, 5) == -EX(1/(4*a**4))*x**4 + EX(1/(3*a**3))*x**3 \
        - EX(1/(2*a**2))*x**2 + EX(1/a)*x + EX(log(a))
    assert rs_log(x + x**2*y + a, x, 4) == -EX(a**(-2))*x**3*y + \
        EX(1/(3*a**3))*x**3 + EX(1/a)*x**2*y - EX(1/(2*a**2))*x**2 + \
        EX(1/a)*x + EX(log(a))

    p = x + x**2 + 3
    assert rs_log(p, x, 10).compose(x, 5) == EX(log(3) + S(19281291595)/9920232)

def mrv(e, x):
    """Returns a SubsSet of most rapidly varying (mrv) subexpressions of 'e',
       and e rewritten in terms of these"""
    e = powsimp(e, deep=True, combine='exp')
    assert isinstance(e, Basic)
    if not e.has(x):
        return SubsSet(), e
    elif e == x:
        s = SubsSet()
        return s, s[x]
    elif e.is_Mul or e.is_Add:
        i, d = e.as_independent(x)  # throw away x-independent terms
        if d.func != e.func:
            s, expr = mrv(d, x)
            return s, e.func(i, expr)
        a, b = d.as_two_terms()
        s1, e1 = mrv(a, x)
        s2, e2 = mrv(b, x)
        return mrv_max1(s1, s2, e.func(i, e1, e2), x)
    elif e.is_Pow:
        b, e = e.as_base_exp()
        if e.has(x):
            return mrv(exp(e * log(b)), x)
        else:
            s, expr = mrv(b, x)
            return s, expr**e
    elif e.func is log:
        s, expr = mrv(e.args[0], x)
        return s, log(expr)
    elif e.func is exp:
        # We know from the theory of this algorithm that exp(log(...)) may always
        # be simplified here, and doing so is vital for termination.
        if e.args[0].func is log:
            return mrv(e.args[0].args[0], x)
        if limitinf(e.args[0], x).is_unbounded:
            s1 = SubsSet()
            e1 = s1[e]
            s2, e2 = mrv(e.args[0], x)
            su = s1.union(s2)[0]
            su.rewrites[e1] = exp(e2)
            return mrv_max3(s1, e1, s2, exp(e2), su, e1, x)
        else:
            s, expr = mrv(e.args[0], x)
            return s, exp(expr)
    elif e.is_Function:
        l = [mrv(a, x) for a in e.args]
        l2 = [s for (s, _) in l if s != SubsSet()]
        if len(l2) != 1:
            # e.g. something like BesselJ(x, x)
            raise NotImplementedError("MRV set computation for functions in"
                                      " several variables not implemented.")
        s, ss = l2[0], SubsSet()
        args = [ss.do_subs(x[1]) for x in l]
        return s, e.func(*args)
    elif e.is_Derivative:
        raise NotImplementedError("MRV set computation for derviatives"
                                  " not implemented yet.")
        return mrv(e.args[0], x)
    raise NotImplementedError(
        "Don't know how to calculate the mrv of '%s'" % e)

    def dig(self):
        """ Number of decimal digits that are guaranteed to be preserved in text.

        When converting text -> float -> text, you are guaranteed that at least ``dig``
        number of digits are preserved with respect to rounding or overflow.
        """
        from sympy.functions import floor, log
        return floor(self.nmant * log(2)/log(10))

def compare(a,b,x):
    """Returns "<" if a<b, "=" for a==b, ">" for a>b"""
    c = limitinf(log(a)/log(b), x)
    if c == 0:
        return "<"
    elif c in [oo,-oo]:
        return ">"
    else:
        return "="

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 decimal_dig(self):
        """ Number of digits needed to store & load without loss.

        Number of decimal digits needed to guarantee that two consecutive conversions
        (float -> text -> float) to be idempotent. This is useful when one do not want
        to loose precision due to rounding errors when storing a floating point value
        as text.
        """
        from sympy.functions import ceiling, log
        return ceiling((self.nmant + 1) * log(2)/log(10) + 1)

def rs_log(p, x, prec):
    """
    The Logarithm of ``p`` modulo ``O(x**prec)``

    Notes
    =====

    truncation of ``integral dx p**-1*d p/dx`` is used.

    Examples
    ========

    >>> from sympy.polys.domains import QQ
    >>> from sympy.polys.rings import ring
    >>> from sympy.polys.ring_series import rs_log
    >>> R, x = ring('x', QQ)
    >>> rs_log(1 + x, x, 8)
    1/7*x**7 - 1/6*x**6 + 1/5*x**5 - 1/4*x**4 + 1/3*x**3 - 1/2*x**2 + x
    >>> rs_log(x**QQ(3, 2) + 1, x, 5)
    1/3*x**(9/2) - 1/2*x**3 + x**(3/2)
    """
    if rs_is_puiseux(p, x):
        return rs_puiseux(rs_log, p, x, prec)
    R = p.ring
    if p == 1:
        return R.zero
    if _has_constant_term(p, x):
        const = 0
        zm = R.zero_monom
        c = p[zm]
        if c == 1:
            pass
        else:
            c_expr = c.as_expr()
            if R.domain is EX:
                const = log(c_expr)
            elif isinstance(c, PolyElement):
                try:
                    const = R(log(c_expr))
                except ValueError:
                    raise DomainError("The given series can't be expanded in "
                                      "this domain.")
            else:
                try:
                    const = R(log(c))
                except ValueError:
                    raise DomainError("The given series can't be expanded in "
                                      "this domain.")

        dlog = p.diff(x)
        dlog = rs_mul(dlog, _series_inversion1(p, x, prec), x, prec - 1)
        return rs_integrate(dlog, x) + const
    else:
        raise NotImplementedError

def test_variable_moment():
    E = Symbol('E')
    I = Symbol('I')

    b = Beam(4, E, 2*(4 - x))
    b.apply_load(20, 4, -1)
    R, M = symbols('R, M')
    b.apply_load(R, 0, -1)
    b.apply_load(M, 0, -2)
    b.bc_deflection = [(0, 0)]
    b.bc_slope = [(0, 0)]
    b.solve_for_reaction_loads(R, M)
    assert b.slope().expand() == ((10*x*SingularityFunction(x, 0, 0)
        - 10*(x - 4)*SingularityFunction(x, 4, 0))/E).expand()
    assert b.deflection().expand() == ((5*x**2*SingularityFunction(x, 0, 0)
        - 10*Piecewise((0, Abs(x)/4 < 1), (16*meijerg(((3, 1), ()), ((), (2, 0)), x/4), True))
        + 40*SingularityFunction(x, 4, 1))/E).expand()

    b = Beam(4, E - x, I)
    b.apply_load(20, 4, -1)
    R, M = symbols('R, M')
    b.apply_load(R, 0, -1)
    b.apply_load(M, 0, -2)
    b.bc_deflection = [(0, 0)]
    b.bc_slope = [(0, 0)]
    b.solve_for_reaction_loads(R, M)
    assert b.slope().expand() == ((-80*(-log(-E) + log(-E + x))*SingularityFunction(x, 0, 0)
        + 80*(-log(-E + 4) + log(-E + x))*SingularityFunction(x, 4, 0) + 20*(-E*log(-E)
        + E*log(-E + x) + x)*SingularityFunction(x, 0, 0) - 20*(-E*log(-E + 4) + E*log(-E + x)
        + x - 4)*SingularityFunction(x, 4, 0))/I).expand()

def test_rubi_integrate():
    assert rubi_integrate(x, x) == x**2/2
    assert rubi_integrate(x**a, x) == x**(a + S(1))/(a + S(1))
    assert rubi_integrate(S(1)/x, x) == log(x)
    assert rubi_integrate(a*x, x) == a*(S(1)/S(2))*x**S(2)
    assert rubi_integrate(1/(x**2*(a + b*x)**2), x) == -b/(a**2*(a + b*x)) - 1/(a**2*x) - 2*b*log(x)/a**3 + 2*b*log(a + b*x)/a**3
    assert rubi_integrate(x**6/(a + b*x)**2, x) == (-a**6/(b**7*(a + b*x)) - S(6)*a**5*log(a + b*x)/b**7 + 5*a**4*x/b**6 - S(2)*a**3*x**2/b**5 + a**2*x**3/b**4 - a*x**4/(S(2)*b**3) + x**5/(S(5)*b**2))
    assert rubi_integrate(1/(x**2*(a + b*x)**2), x) == -b/(a**2*(a + b*x)) - 1/(a**2*x) - 2*b*log(x)/a**3 + 2*b*log(a + b*x)/a**3
    assert rubi_integrate(a + S(1)/x, x) == a*x + log(x)
    assert rubi_integrate((a + b*x)**2/x**3, x) == -a**2/(2*x**2) - 2*a*b/x + b**2*log(x)
    assert rubi_integrate(a**3*x, x) == S(1)/S(2)*a**3*x**2
    assert rubi_integrate((a + b*x)**3/x**3, x) == -a**3/(2*x**2) - 3*a**2*b/x + 3*a*b**2*log(x) + b**3*x
    assert rubi_integrate(x**3*(a + b*x), x) == a*x**4/4 + b*x**5/5
    assert rubi_integrate((b*x)**m*(d*x + 2)**n, x) == 2**n*(b*x)**(m + 1)*hyper((-n, m + 1), (m + 2,), -d*x/2)/(b*(m + 1))
    assert rubi_test(rubi_integrate(1/(1 + x**5), x), x, log(x + S(1))/S(5) + S(2)*Sum(-log((S(2)*x - S(2)*cos(pi*(S(2)*k/S(5) + S(-1)/5)))**S(2) - S(4)*sin(S(2)*pi*k/S(5) + S(3)*pi/S(10))**S(2) + S(4))*cos(pi*(S(2)*k/S(5) + S(-1)/5))/S(2) - (-S(2)*cos(pi*(S(2)*k/S(5) + S(-1)/5))**S(2) + S(2))*atan((-x/cos(pi*(S(2)*k/S(5) + S(-1)/5)) + S(1))/sqrt(-(cos(S(2)*pi*k/S(5) - pi/S(5)) + S(-1))*(cos(S(2)*pi*k/S(5) - pi/S(5)) + S(1))/cos(S(2)*pi*k/S(5) - pi/S(5))**S(2)))/(S(2)*sqrt(-(cos(S(2)*pi*k/S(5) - pi/S(5)) + S(-1))*(cos(S(2)*pi*k/S(5) - pi/S(5)) + S(1))/cos(S(2)*pi*k/S(5) - pi/S(5))**S(2))*cos(pi*(S(2)*k/S(5) + S(-1)/5))), (k, S(1), S(2)))/S(5), _numerical=True)

def test_bounded():
    x, y = symbols('x,y')
    assert ask(Q.bounded(x)) == False
    assert ask(Q.bounded(x), Q.bounded(x)) == True
    assert ask(Q.bounded(x), Q.bounded(y)) == False
    assert ask(Q.bounded(x), Q.complex(x)) == False

    assert ask(Q.bounded(x+1)) == False
    assert ask(Q.bounded(x+1), Q.bounded(x)) == True
    assert ask(Q.bounded(x+y)) == None
    assert ask(Q.bounded(x+y), Q.bounded(x)) == False
    assert ask(Q.bounded(x+1), Q.bounded(x) & Q.bounded(y)) == True

    assert ask(Q.bounded(2*x)) == False
    assert ask(Q.bounded(2*x), Q.bounded(x)) == True
    assert ask(Q.bounded(x*y)) == None
    assert ask(Q.bounded(x*y), Q.bounded(x)) == False
    assert ask(Q.bounded(x*y), Q.bounded(x) & Q.bounded(y)) == True

    assert ask(Q.bounded(x**2)) == False
    assert ask(Q.bounded(2**x)) == False
    assert ask(Q.bounded(2**x), Q.bounded(x)) == True
    assert ask(Q.bounded(x**x)) == False
    assert ask(Q.bounded(Rational(1,2) ** x)) == None
    assert ask(Q.bounded(Rational(1,2) ** x), Q.positive(x)) == True
    assert ask(Q.bounded(Rational(1,2) ** x), Q.negative(x)) == False
    assert ask(Q.bounded(sqrt(x))) == False

    # sign function
    assert ask(Q.bounded(sign(x))) == True
    assert ask(Q.bounded(sign(x)), ~Q.bounded(x)) == True

    # exponential functions
    assert ask(Q.bounded(log(x))) == False
    assert ask(Q.bounded(log(x)), Q.bounded(x)) == True
    assert ask(Q.bounded(exp(x))) == False
    assert ask(Q.bounded(exp(x)), Q.bounded(x)) == True
    assert ask(Q.bounded(exp(2))) == True

    # trigonometric functions
    assert ask(Q.bounded(sin(x))) == True
    assert ask(Q.bounded(sin(x)), ~Q.bounded(x)) == True
    assert ask(Q.bounded(cos(x))) == True
    assert ask(Q.bounded(cos(x)), ~Q.bounded(x)) == True
    assert ask(Q.bounded(2*sin(x))) == True
    assert ask(Q.bounded(sin(x)**2)) == True
    assert ask(Q.bounded(cos(x)**2)) == True
    assert ask(Q.bounded(cos(x) + sin(x))) == True

def test_bounded():
    x, y = symbols('xy')
    assert ask(x, Q.bounded) == False
    assert ask(x, Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(x, Q.bounded, Assume(y, Q.bounded)) == False
    assert ask(x, Q.bounded, Assume(x, Q.complex)) == False

    assert ask(x+1, Q.bounded) == False
    assert ask(x+1, Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(x+y, Q.bounded) == None
    assert ask(x+y, Q.bounded, Assume(x, Q.bounded)) == False
    assert ask(x+1, Q.bounded, Assume(x, Q.bounded) & \
                Assume(y, Q.bounded)) == True

    assert ask(2*x, Q.bounded) == False
    assert ask(2*x, Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(x*y, Q.bounded) == None
    assert ask(x*y, Q.bounded, Assume(x, Q.bounded)) == False
    assert ask(x*y, Q.bounded, Assume(x, Q.bounded) & \
                Assume(y, Q.bounded)) == True

    assert ask(x**2, Q.bounded) == False
    assert ask(2**x, Q.bounded) == False
    assert ask(2**x, Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(x**x, Q.bounded) == False
    assert ask(Rational(1,2) ** x, Q.bounded) == True
    assert ask(x ** Rational(1,2), Q.bounded) == False

    # sign function
    assert ask(sign(x), Q.bounded) == True
    assert ask(sign(x), Q.bounded, Assume(x, Q.bounded, False)) == True

    # exponential functions
    assert ask(log(x), Q.bounded) == False
    assert ask(log(x), Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(exp(x), Q.bounded) == False
    assert ask(exp(x), Q.bounded, Assume(x, Q.bounded)) == True
    assert ask(exp(2), Q.bounded) == True

    # trigonometric functions
    assert ask(sin(x), Q.bounded) == True
    assert ask(sin(x), Q.bounded, Assume(x, Q.bounded, False)) == True
    assert ask(cos(x), Q.bounded) == True
    assert ask(cos(x), Q.bounded, Assume(x, Q.bounded, False)) == True
    assert ask(2*sin(x), Q.bounded) == True
    assert ask(sin(x)**2, Q.bounded) == True
    assert ask(cos(x)**2, Q.bounded) == True
    assert ask(cos(x) + sin(x), Q.bounded) == True

def test_catalan():
    n = Symbol('n', integer=True)
    m = Symbol('n', integer=True, positive=True)

    catalans = [1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786]
    for i, c in enumerate(catalans):
        assert catalan(i) == c
        assert catalan(n).rewrite(factorial).subs(n, i) == c
        assert catalan(n).rewrite(Product).subs(n, i).doit() == c

    assert catalan(x) == catalan(x)
    assert catalan(2*x).rewrite(binomial) == binomial(4*x, 2*x)/(2*x + 1)
    assert catalan(Rational(1, 2)).rewrite(gamma) == 8/(3*pi)
    assert catalan(Rational(1, 2)).rewrite(factorial).rewrite(gamma) ==\
        8 / (3 * pi)
    assert catalan(3*x).rewrite(gamma) == 4**(
        3*x)*gamma(3*x + Rational(1, 2))/(sqrt(pi)*gamma(3*x + 2))
    assert catalan(x).rewrite(hyper) == hyper((-x + 1, -x), (2,), 1)

    assert catalan(n).rewrite(factorial) == factorial(2*n) / (factorial(n + 1)
                                                              * factorial(n))
    assert isinstance(catalan(n).rewrite(Product), catalan)
    assert isinstance(catalan(m).rewrite(Product), Product)

    assert diff(catalan(x), x) == (polygamma(
        0, x + Rational(1, 2)) - polygamma(0, x + 2) + log(4))*catalan(x)

    assert catalan(x).evalf() == catalan(x)
    c = catalan(S.Half).evalf()
    assert str(c) == '0.848826363156775'
    c = catalan(I).evalf(3)
    assert str((re(c), im(c))) == '(0.398, -0.0209)'

def mrv(e, x):
    "Returns a python set of  most rapidly varying (mrv) subexpressions of 'e'"
    assert isinstance(e, Basic)
    if not e.has(x):
        return set([])
    elif e == x:
        return set([x])
    elif e.is_Mul:
        a, b = e.as_two_terms()
        return mrv_max(mrv(a,x), mrv(b,x), x)
    elif e.is_Add:
        a, b = e.as_two_terms()
        return mrv_max(mrv(a,x), mrv(b,x), x)
    elif e.is_Pow:
        if e.exp.has(x):
            return mrv(exp(e.exp * log(e.base)), x)
        else:
            return mrv(e.base, x)
    elif e.func is log:
        return mrv(e.args[0], x)
    elif e.func is exp:
        if limitinf(e.args[0], x) in [oo,-oo]:
            return mrv_max(set([e]), mrv(e.args[0], x), x)
        else:
            return mrv(e.args[0], x)
    elif e.is_Function:
        if len(e.args) == 1:
            return mrv(e.args[0], x)
        #only functions of 1 argument currently implemented
        raise NotImplementedError("Functions with more arguments: '%s'" % e)
    raise NotImplementedError("Don't know how to calculate the mrv of '%s'" % e)

def mrv_leadterm(e, x):
    """Returns (c0, e0) for e."""
    Omega = SubsSet()
    if not e.has(x):
        return (e, S.Zero)
    if Omega == SubsSet():
        Omega, exps = mrv(e, x)
    if not Omega:
        # e really does not depend on x after simplification
        series = calculate_series(e, x)
        c0, e0 = series.leadterm(x)
        assert e0 == 0
        return c0, e0
    if x in Omega:
        #move the whole omega up (exponentiate each term):
        Omega_up = moveup2(Omega, x)
        e_up = moveup([e], x)[0]
        exps_up = moveup([exps], x)[0]
        # NOTE: there is no need to move this down!
        e = e_up
        Omega = Omega_up
        exps = exps_up
    #
    # The positive dummy, w, is used here so log(w*2) etc. will expand;
    # a unique dummy is needed in this algorithm
    #
    # For limits of complex functions, the algorithm would have to be
    # improved, or just find limits of Re and Im components separately.
    #
    w = Dummy("w", real=True, positive=True, bounded=True)
    f, logw = rewrite(exps, Omega, x, w)
    series = calculate_series(f, w, logx=logw)
    series = series.subs(log(w), logw) # this should not be necessary
    return series.leadterm(w)

def mrv_leadterm(e, x, Omega=[]):
    """Returns (c0, e0) for e."""
    if not e.has(x):
        return (e, S.Zero)
    Omega = [t for t in Omega if subexp(e, t)]
    if Omega == []:
        Omega = mrv(e, x)
    if not Omega:
        # e really does not depend on x after simplification
        series = calculate_series(e, x)
        c0, e0 = series.leadterm(x)
        assert e0 == 0
        return c0, e0
    if x in set(Omega):
        # move the whole omega up (exponentiate each term):
        Omega_up = set(moveup(Omega, x))
        e_up = moveup([e], x)[0]
        # calculate the lead term
        mrv_leadterm_up = mrv_leadterm(e_up, x, Omega_up)
        # move the result (c0, e0) down
        return tuple(movedown(mrv_leadterm_up, x))
    #
    # The positive dummy, w, is used here so log(w*2) etc. will expand;
    # a unique dummy is needed in this algorithm
    #
    # For limits of complex functions, the algorithm would have to be
    # improved, or just find limits of Re and Im components separately.
    #
    w = Dummy("w", real=True, positive=True)
    f, logw = rewrite(e, set(Omega), x, w)
    series = calculate_series(f, w)
    series = series.subs(log(w), logw)
    return series.leadterm(w)

def compare(a, b, x):
    """Returns "<" if a<b, "=" for a == b, ">" for a>b"""
    # log(exp(...)) must always be simplified here for termination
    la, lb = log(a), log(b)
    if isinstance(a, Basic) and a.func is exp:
        la = a.args[0]
    if isinstance(b, Basic) and b.func is exp:
        lb = b.args[0]

    c = limitinf(la/lb, x)
    if c == 0:
        return "<"
    elif c.is_unbounded:
        return ">"
    else:
        return "="

def test_trigintegrate_odd():
    assert trigintegrate(Rational(1), x) == x
    assert trigintegrate(x, x) is None
    assert trigintegrate(x**2, x) is None

    assert trigintegrate(sin(x), x) == -cos(x)
    assert trigintegrate(cos(x), x) == sin(x)

    assert trigintegrate(sin(3*x), x) == -cos(3*x)/3
    assert trigintegrate(cos(3*x), x) == sin(3*x)/3

    y = Symbol('y')
    assert trigintegrate(sin(y*x), x) == \
        Piecewise((0, Eq(y, 0)), (-cos(y*x)/y, True))
    assert trigintegrate(cos(y*x), x) == \
        Piecewise((x, Eq(y, 0)), (sin(y*x)/y, True))
    assert trigintegrate(sin(y*x)**2, x) == \
        Piecewise((0, Eq(y, 0)), ((x*y/2 - sin(x*y)*cos(x*y)/2)/y, True))
    assert trigintegrate(sin(y*x)*cos(y*x), x) == \
        Piecewise((0, Eq(y, 0)), (sin(x*y)**2/(2*y), True))
    assert trigintegrate(cos(y*x)**2, x) == \
        Piecewise((x, Eq(y, 0)), ((x*y/2 + sin(x*y)*cos(x*y)/2)/y, True))

    y = Symbol('y', positive=True)
    # TODO: remove conds='none' below. For this to work we would have to rule
    #       out (e.g. by trying solve) the condition y = 0, incompatible with
    #       y.is_positive being True.
    assert trigintegrate(sin(y*x), x, conds='none') == -cos(y*x)/y
    assert trigintegrate(cos(y*x), x, conds='none') == sin(y*x)/y

    assert trigintegrate(sin(x)*cos(x), x) == sin(x)**2/2
    assert trigintegrate(sin(x)*cos(x)**2, x) == -cos(x)**3/3
    assert trigintegrate(sin(x)**2*cos(x), x) == sin(x)**3/3

    # check if it selects right function to substitute,
    # so the result is kept simple
    assert trigintegrate(sin(x)**7 * cos(x), x) == sin(x)**8/8
    assert trigintegrate(sin(x) * cos(x)**7, x) == -cos(x)**8/8

    assert trigintegrate(sin(x)**7 * cos(x)**3, x) == \
        -sin(x)**10/10 + sin(x)**8/8
    assert trigintegrate(sin(x)**3 * cos(x)**7, x) == \
        cos(x)**10/10 - cos(x)**8/8

    # both n, m are odd and -ve, and not necessarily equal
    assert trigintegrate(sin(x)**-1*cos(x)**-1, x) == \
        -log(sin(x)**2 - 1)/2 + log(sin(x))