def _rational_case(cls, poly, func):
        """Handle the rational function case. """
        roots = symbols('r:%d' % poly.degree())
        var, expr = func.variables[0], func.expr

        f = sum(expr.subs(var, r) for r in roots)
        p, q = together(f).as_numer_denom()

        domain = QQ[roots]

        p = p.expand()
        q = q.expand()

        try:
            p = Poly(p, domain=domain, expand=False)
        except GeneratorsNeeded:
            p, p_coeff = None, (p,)
        else:
            p_monom, p_coeff = zip(*p.terms())

        try:
            q = Poly(q, domain=domain, expand=False)
        except GeneratorsNeeded:
            q, q_coeff = None, (q,)
        else:
            q_monom, q_coeff = zip(*q.terms())

        coeffs, mapping = symmetrize(p_coeff + q_coeff, formal=True)
        formulas, values = viete(poly, roots), []

        for (sym, _), (_, val) in zip(mapping, formulas):
            values.append((sym, val))

        for i, (coeff, _) in enumerate(coeffs):
            coeffs[i] = coeff.subs(values)

        n = len(p_coeff)

        p_coeff = coeffs[:n]
        q_coeff = coeffs[n:]

        if p is not None:
            p = Poly(dict(zip(p_monom, p_coeff)), *p.gens).as_expr()
        else:
            (p,) = p_coeff

        if q is not None:
            q = Poly(dict(zip(q_monom, q_coeff)), *q.gens).as_expr()
        else:
            (q,) = q_coeff

        return factor(p/q)

def _quintic_simplify(expr):
    expr = powsimp(expr)
    expr = cancel(expr)
    return together(expr)

def test_together():
    assert together(0) == 0
    assert together(1) == 1

    assert together(x*y*z) == x*y*z
    assert together(x + y) == x + y

    assert together(1/x) == 1/x

    assert together(1/x + 1) == (x + 1)/x
    assert together(1/x + 3) == (3*x + 1)/x
    assert together(1/x + x) == (x**2 + 1)/x

    assert together(1/x + Rational(1, 2)) == (x + 2)/(2*x)
    assert together(Rational(1, 2) + x/2) == Mul(S.Half, x + 1, evaluate=False)

    assert together(1/x + 2/y) == (2*x + y)/(y*x)
    assert together(1/(1 + 1/x)) == x/(1 + x)
    assert together(x/(1 + 1/x)) == x**2/(1 + x)

    assert together(1/x + 1/y + 1/z) == (x*y + x*z + y*z)/(x*y*z)
    assert together(1/(1 + x + 1/y + 1/z)) == y*z/(y + z + y*z + x*y*z)

    assert together(1/(x*y) + 1/(x*y)**2) == y**(-2)*x**(-2)*(1 + x*y)
    assert together(1/(x*y) + 1/(x*y)**4) == y**(-4)*x**(-4)*(1 + x**3*y**3)
    assert together(1/(x**7*y) + 1/(x*y)**4) == y**(-4)*x**(-7)*(x**3 + y**3)

    assert together(5/(2 + 6/(3 + 7/(4 + 8/(5 + 9/x))))) == \
         (S(5)/2)*((171 + 119*x)/(279 + 203*x))

    assert together(1 + 1/(x + 1)**2) == (1 + (x + 1)**2)/(x + 1)**2
    assert together(1 + 1/(x*(1 + x))) == (1 + x*(1 + x))/(x*(1 + x))
    assert together(1/(x*(x + 1)) + 1/(x*(x + 2))) == (3 + 2*x)/(x*(1 + x)*(2 + x))
    assert together(1 + 1/(2*x + 2)**2) == (4*(x + 1)**2 + 1)/(4*(x + 1)**2)

    assert together(sin(1/x + 1/y)) == sin(1/x + 1/y)
    assert together(sin(1/x + 1/y), deep=True) == sin((x + y)/(x*y))

    assert together(1/exp(x) + 1/(x*exp(x))) == (1+x)/(x*exp(x))
    assert together(1/exp(2*x) + 1/(x*exp(3*x))) == (1+exp(x)*x)/(x*exp(3*x))

    assert together(Integral(1/x + 1/y, x)) == Integral((x + y)/(x*y), x)
    assert together(Eq(1/x + 1/y, 1 + 1/z)) == Eq((x + y)/(x*y), (z + 1)/z)

    assert together((A*B)**-1 + (B*A)**-1) == (A*B)**-1 + (B*A)**-1

def continuous_domain(f, symbol, domain):
    """
    Returns the intervals in the given domain for which the function
    is continuous.
    This method is limited by the ability to determine the various
    singularities and discontinuities of the given function.

    Parameters
    ==========

    f : Expr
        The concerned function.
    symbol : Symbol
        The variable for which the intervals are to be determined.
    domain : Interval
        The domain over which the continuity of the symbol has to be checked.

    Examples
    ========

    >>> from sympy import Symbol, S, tan, log, pi, sqrt
    >>> from sympy.sets import Interval
    >>> from sympy.calculus.util import continuous_domain
    >>> x = Symbol('x')
    >>> continuous_domain(1/x, x, S.Reals)
    Union(Interval.open(-oo, 0), Interval.open(0, oo))
    >>> continuous_domain(tan(x), x, Interval(0, pi))
    Union(Interval.Ropen(0, pi/2), Interval.Lopen(pi/2, pi))
    >>> continuous_domain(sqrt(x - 2), x, Interval(-5, 5))
    Interval(2, 5)
    >>> continuous_domain(log(2*x - 1), x, S.Reals)
    Interval.open(1/2, oo)

    Returns
    =======

    Interval
        Union of all intervals where the function is continuous.

    Raises
    ======
    NotImplementedError
        If the method to determine continuity of such a function
        has not yet been developed.

    """
    from sympy.solvers.inequalities import solve_univariate_inequality
    from sympy.solvers.solveset import solveset, _has_rational_power

    if domain.is_subset(S.Reals):
        constrained_interval = domain
        for atom in f.atoms(Pow):
            predicate, denomin = _has_rational_power(atom, symbol)
            constraint = S.EmptySet
            if predicate and denomin == 2:
                constraint = solve_univariate_inequality(atom.base >= 0,
                                                         symbol).as_set()
                constrained_interval = Intersection(constraint,
                                                    constrained_interval)

        for atom in f.atoms(log):
            constraint = solve_univariate_inequality(atom.args[0] > 0,
                                                     symbol).as_set()
            constrained_interval = Intersection(constraint,
                                                constrained_interval)

        domain = constrained_interval

    try:
        sings = S.EmptySet
        if f.has(Abs):
            sings = solveset(1/f, symbol, domain) + \
                solveset(denom(together(f)), symbol, domain)
        else:
            for atom in f.atoms(Pow):
                predicate, denomin = _has_rational_power(atom, symbol)
                if predicate and denomin == 2:
                    sings = solveset(1/f, symbol, domain) +\
                        solveset(denom(together(f)), symbol, domain)
                    break
            else:
                sings = Intersection(solveset(1/f, symbol), domain) + \
                    solveset(denom(together(f)), symbol, domain)

    except NotImplementedError:
        import sys
        raise (NotImplementedError("Methods for determining the continuous domains"
                                   " of this function have not been developed."),
               None,
               sys.exc_info()[2])

    return domain - sings