def _sqrtdenest34(expr):
    """Helper that denests the square root of three or four surds.

    Examples
    ========
    >>> from sympy import sqrt
    >>> from sympy.simplify.sqrtdenest import _sqrtdenest34
    >>> _sqrtdenest34(sqrt(-72*sqrt(2) + 158*sqrt(5) + 498))
    -sqrt(10) + sqrt(2) + 9 + 9*sqrt(5)
    >>> _sqrtdenest34(sqrt(12 + 2*sqrt(6) + 2*sqrt(14) + 2*sqrt(21)))
    sqrt(2) + sqrt(3) + sqrt(7)

    References
    ==========
    - D. J. Jeffrey and A. D. Rich, 'Symplifying Square Roots of Square Roots
      by Denesting'

    """
    from sympy.simplify.simplify import radsimp
    if expr.base < 0:
        return sqrt(-1)*_sqrtdenest34(sqrt(-expr.base))
    # a should be > b so we sort; this also makes the process canonical
    args = sorted(expr.base.args, reverse=True)
    a = Add._from_args(args[:2])
    b = Add._from_args(args[2:])
    c = _sqrtdenest1(sqrt(_mexpand(a**2 - b**2)))
    if sqrt_depth(c) > 1:
        return expr
    d = _sqrtdenest1(sqrt(a + c))
    if sqrt_depth(d) > 1:
        return expr
    return _mexpand(radsimp((d/sqrt(2) + b/(sqrt(2)*d))))

def _dict_from_expr(expr, opt):
    """Transform an expression into a multinomial form. """
    if expr.is_commutative is False:
        raise PolynomialError('non-commutative expressions are not supported')

    def _is_expandable_pow(expr):
        return (expr.is_Pow and expr.exp.is_positive and expr.exp.is_Integer
                and expr.base.is_Add)

    if opt.expand is not False:
        expr = expr.expand()
        # TODO: Integrate this into expand() itself
        while any(_is_expandable_pow(i) or i.is_Mul and
            any(_is_expandable_pow(j) for j in i.args) for i in
                Add.make_args(expr)):

            expr = expand_multinomial(expr)
        while any(i.is_Mul and any(j.is_Add for j in i.args) for i in Add.make_args(expr)):
            expr = expand_mul(expr)

    if opt.gens:
        rep, gens = _dict_from_expr_if_gens(expr, opt)
    else:
        rep, gens = _dict_from_expr_no_gens(expr, opt)

    return rep, opt.clone({'gens': gens})

def interpolate(data, x):
    """
    Construct an interpolating polynomial for the data points.

    Examples
    ========

    >>> from sympy.polys.polyfuncs import interpolate
    >>> from sympy.abc import x

    A list is interpreted as though it were paired with a range starting
    from 1:

    >>> interpolate([1, 4, 9, 16], x)
    x**2

    This can be made explicit by giving a list of coordinates:

    >>> interpolate([(1, 1), (2, 4), (3, 9)], x)
    x**2

    The (x, y) coordinates can also be given as keys and values of a
    dictionary (and the points need not be equispaced):

    >>> interpolate([(-1, 2), (1, 2), (2, 5)], x)
    x**2 + 1
    >>> interpolate({-1: 2, 1: 2, 2: 5}, x)
    x**2 + 1

    """
    n = len(data)
    poly = None

    if isinstance(data, dict):
        X, Y = list(zip(*data.items()))
        poly = interpolating_poly(n, x, X, Y)
    else:
        if isinstance(data[0], tuple):
            X, Y = list(zip(*data))
            poly = interpolating_poly(n, x, X, Y)
        else:
            Y = list(data)

            numert = Mul(*[(x - i) for i in range(1, n + 1)])
            denom = -factorial(n - 1) if n%2 == 0 else factorial(n - 1)
            coeffs = []
            for i in range(1, n + 1):
                coeffs.append(numert/(x - i)/denom)
                denom = denom/(i - n)*i

            poly = Add(*[coeff*y for coeff, y in zip(coeffs, Y)])

    return poly.expand()

    def eval(cls, arg):
        from sympy import AccumBounds, im

        def _eval(arg):
            if arg is S.Infinity or arg is S.NegativeInfinity:
                return AccumBounds(0, 1)
            if arg.is_integer:
                return S.Zero
            if arg.is_number:
                if arg is S.NaN:
                    return S.NaN
                elif arg is S.ComplexInfinity:
                    return None
                else:
                    return arg - floor(arg)
            return cls(arg, evaluate=False)

        terms = Add.make_args(arg)
        real, imag = S.Zero, S.Zero
        for t in terms:
            # Two checks are needed for complex arguments
            # see issue-7649 for details
            if t.is_imaginary or (S.ImaginaryUnit*t).is_real:
                i = im(t)
                if not i.has(S.ImaginaryUnit):
                    imag += i
                else:
                    real += t
            else:
                real += t

        real = _eval(real)
        imag = _eval(imag)
        return real + S.ImaginaryUnit*imag

    def eval(cls, arg):
        if arg is S.NaN:
            return S.NaN
        elif arg.is_real:
            return S.Zero
        elif arg.is_imaginary:
            return -S.ImaginaryUnit * arg
        elif arg.is_Function and arg.func == conjugate:
            return -im(arg.args[0])
        else:
            included, reverted, excluded = [], [], []
            args = Add.make_args(arg)
            for term in args:
                coeff = term.as_coefficient(S.ImaginaryUnit)

                if coeff is not None:
                    if not coeff.is_real:
                        reverted.append(coeff)
                    else:
                        excluded.append(coeff)
                elif term.has(S.ImaginaryUnit) or not term.is_real:
                    # Try to do some advanced expansion.  If
                    # impossible, don't try to do im(arg) again
                    # (because this is what we are trying to do now).
                    real_imag = term.as_real_imag(ignore=arg)
                    if real_imag:
                        excluded.append(real_imag[1])
                    else:
                        included.append(term)

            if len(args) != len(included):
                a, b, c = map(lambda xs: Add(*xs),
                    [included, reverted, excluded])

                return cls(a) + re(b) + c

    def eval(cls, arg):
        if arg is S.NaN:
            return S.NaN
        elif arg.is_real:
            return arg
        elif arg.is_Function and arg.func == conjugate:
            return re(arg.args[0])
        else:

            included, reverted, excluded = [], [], []
            arg = Add.make_args(arg)
            for term in arg:
                coeff = term.as_coefficient(S.ImaginaryUnit)

                if coeff is not None:
                    if not coeff.is_real:
                        reverted.append(coeff)
                elif not term.has(S.ImaginaryUnit) and term.is_real:
                    excluded.append(term)
                else:
                    included.append(term)

            if len(arg) != len(included):
                a, b, c = map(lambda xs: Add(*xs),
                    [included, reverted, excluded])

                return cls(a) - im(b) + c

def _set_function(f, self):
    expr = f.expr
    if not isinstance(expr, Expr):
        return

    if len(f.variables) > 1:
        return

    n = f.variables[0]

    # f(x) + c and f(-x) + c cover the same integers
    # so choose the form that has the fewest negatives
    c = f(0)
    fx = f(n) - c
    f_x = f(-n) - c
    neg_count = lambda e: sum(_coeff_isneg(_) for _ in Add.make_args(e))
    if neg_count(f_x) < neg_count(fx):
        expr = f_x + c

    a = Wild('a', exclude=[n])
    b = Wild('b', exclude=[n])
    match = expr.match(a*n + b)
    if match and match[a]:
        # canonical shift
        expr = match[a]*n + match[b] % match[a]

    if expr != f.expr:
        return ImageSet(Lambda(n, expr), S.Integers)

    def doit(self, **hints):
        """Evaluates limit"""
        e, z, z0, dir = self.args

        if hints.get('deep', True):
            e = e.doit(**hints)
            z = z.doit(**hints)
            z0 = z0.doit(**hints)

        if e == z:
            return z0

        if not e.has(z):
            return e

        # gruntz fails on factorials but works with the gamma function
        # If no factorial term is present, e should remain unchanged.
        # factorial is defined to be zero for negative inputs (which
        # differs from gamma) so only rewrite for positive z0.
        if z0.is_positive:
            e = e.rewrite(factorial, gamma)

        if e.is_Mul:
            if abs(z0) is S.Infinity:
                # XXX todo: this should probably be stated in the
                # negative -- i.e. to exclude expressions that should
                # not be handled this way but I'm not sure what that
                # condition is; when ok is True it means that the leading
                # term approach is going to succeed (hopefully)
                ok = lambda w: (z in w.free_symbols and
                                any(a.is_polynomial(z) or
                                    any(z in m.free_symbols and m.is_polynomial(z)
                                        for m in Mul.make_args(a))
                                    for a in Add.make_args(w)))
                if all(ok(w) for w in e.as_numer_denom()):
                    u = Dummy(positive=(z0 is S.Infinity))
                    inve = e.subs(z, 1/u)
                    r = limit(inve.as_leading_term(u), u,
                              S.Zero, "+" if z0 is S.Infinity else "-")
                    if isinstance(r, Limit):
                        return self
                    else:
                        return r

        if e.is_Order:
            return Order(limit(e.expr, z, z0), *e.args[1:])

        try:
            r = gruntz(e, z, z0, dir)
            if r is S.NaN:
                raise PoleError()
        except (PoleError, ValueError):
            r = heuristics(e, z, z0, dir)
            if r is None:
                return self

        return r

def apart(f, x=None, full=False, **options):
    """
    Compute partial fraction decomposition of a rational function.

    Given a rational function ``f`` compute partial fraction decomposition
    of ``f``. Two algorithms are available: one is based on undetermined
    coefficients method and the other is Bronstein's full partial fraction
    decomposition algorithm.

    Examples
    ========

    >>> from sympy.polys.partfrac import apart
    >>> from sympy.abc import x, y

    >>> apart(y/(x + 2)/(x + 1), x)
    -y/(x + 2) + y/(x + 1)

    """
    allowed_flags(options, [])

    f = sympify(f)

    if f.is_Atom:
        return f
    else:
        P, Q = f.as_numer_denom()

    options = set_defaults(options, extension=True)
    (P, Q), opt = parallel_poly_from_expr((P, Q), x, **options)

    if P.is_multivariate:
        raise NotImplementedError(
            "multivariate partial fraction decomposition")

    common, P, Q = P.cancel(Q)

    poly, P = P.div(Q, auto=True)
    P, Q = P.rat_clear_denoms(Q)

    if Q.degree() <= 1:
        partial = P/Q
    else:
        if not full:
            partial = apart_undetermined_coeffs(P, Q)
        else:
            partial = apart_full_decomposition(P, Q)

    terms = S.Zero

    for term in Add.make_args(partial):
        if term.has(RootSum):
            terms += term
        else:
            terms += factor(term)

    return common*(poly.as_expr() + terms)

    def __new__(cls, expr, *symbols, **assumptions):

        expr = sympify(expr).expand()
        if expr is S.NaN:
            return S.NaN

        if symbols:
            symbols = map(sympify, symbols)
            if not all(isinstance(s, Symbol) for s in symbols):
                raise NotImplementedError('Order at points other than 0 not supported.')
        else:
            symbols = list(expr.free_symbols)

        if expr.is_Order:

            new_symbols = list(expr.variables)
            for s in symbols:
                if s not in new_symbols:
                    new_symbols.append(s)
            if len(new_symbols) == len(expr.variables):
                return expr
            symbols = new_symbols

        elif symbols:

            if expr.is_Add:
                lst = expr.extract_leading_order(*symbols)
                expr = Add(*[f.expr for (e,f) in lst])
            elif expr:
                expr = expr.as_leading_term(*symbols)
                coeff, terms = expr.as_coeff_mul()
                expr = Mul(*[t for t in terms if t.has(*symbols)])

        if expr is S.Zero:
            return expr
        elif not expr.has(*symbols):
            expr = S.One

        # create Order instance:
        symbols.sort(Basic.compare)
        obj = Expr.__new__(cls, expr, *symbols, **assumptions)

        return obj

    def __new__(cls, *args, **kwargs):
        args = list(map(sympify, args))
        check = kwargs.get('check', False)

        obj = Basic.__new__(cls, *args)
        if check:
            if all(not isinstance(i, MatrixExpr) for i in args):
                return Add.fromiter(args)
            validate(*args)
        return obj

def sum_simplify(s):
    """Main function for Sum simplification"""
    from sympy.concrete.summations import Sum
    from sympy.core.function import expand

    terms = Add.make_args(expand(s))
    s_t = [] # Sum Terms
    o_t = [] # Other Terms

    for term in terms:
        if isinstance(term, Mul):
            other = 1
            sum_terms = []

            if not term.has(Sum):
                o_t.append(term)
                continue

            mul_terms = Mul.make_args(term)
            for mul_term in mul_terms:
                if isinstance(mul_term, Sum):
                    r = mul_term._eval_simplify()
                    sum_terms.extend(Add.make_args(r))
                else:
                    other = other * mul_term
            if len(sum_terms):
                #some simplification may have happened
                #use if so
                s_t.append(Mul(*sum_terms) * other)
            else:
                o_t.append(other)
        elif isinstance(term, Sum):
            #as above, we need to turn this into an add list
            r = term._eval_simplify()
            s_t.extend(Add.make_args(r))
        else:
            o_t.append(term)


    result = Add(sum_combine(s_t), *o_t)

    return result

def sum_simplify(s):
    """Main function for Sum simplification"""
    from sympy.concrete.summations import Sum

    terms = Add.make_args(s)
    s_t = [] # Sum Terms
    o_t = [] # Other Terms

    for term in terms:
        if isinstance(term, Mul):
            constant = 1
            other = 1
            s = 0
            n_sum_terms = 0
            for j in range(len(term.args)):
                if isinstance(term.args[j], Sum):
                    s = term.args[j]
                    n_sum_terms = n_sum_terms + 1
                elif term.args[j].is_number == True:
                    constant = constant * term.args[j]
                else:
                    other = other * term.args[j]
            if other == 1 and n_sum_terms == 1:
                # Insert the constant inside the Sum
                s_t.append(Sum(constant * s.function, *s.limits))
            elif other != 1 and n_sum_terms == 1:
                o_t.append(other * Sum(constant * s.function, *s.limits))
            else:
                o_t.append(term)
        elif isinstance(term, Sum):
            s_t.append(term)
        else:
            o_t.append(term)

    used = [False] * len(s_t)

    for method in range(2):
        for i, s_term1 in enumerate(s_t):
            if not used[i]:
                for j, s_term2 in enumerate(s_t):
                    if not used[j] and i != j:
                        temp = sum_add(s_term1, s_term2, method)
                        if isinstance(temp, Sum):
                            s_t[i] = temp
                            s_term1 = s_t[i]
                            used[j] = True

    result = Add(*o_t)

    for i, s_term in enumerate(s_t):
        if not used[i]:
            result = Add(result, s_term)

    return result

    def as_terms(self):
        """Transform an expression to a list of terms. """
        from sympy.core import Add, Mul, S
        from sympy.core.exprtools import decompose_power

        gens, terms = set([]), []

        for term in Add.make_args(self):
            coeff, _term = term.as_coeff_Mul()

            coeff = complex(coeff)
            cpart, ncpart = {}, []

            if _term is not S.One:
                for factor in Mul.make_args(_term):
                    if factor.is_number:
                        try:
                            coeff *= complex(factor)
                        except ValueError:
                            pass
                        else:
                            continue

                    if factor.is_commutative:
                        base, exp = decompose_power(factor)

                        cpart[base] = exp
                        gens.add(base)
                    else:
                        ncpart.append(factor)

            coeff = coeff.real, coeff.imag
            ncpart = tuple(ncpart)

            terms.append((term, (coeff, cpart, ncpart)))

        gens = sorted(gens, key=Basic.sorted_key)

        k, indices = len(gens), {}

        for i, g in enumerate(gens):
            indices[g] = i

        result = []

        for term, (coeff, cpart, ncpart) in terms:
            monom = [0]*k

            for base, exp in cpart.iteritems():
                monom[indices[base]] = exp

            result.append((term, (coeff, tuple(monom), ncpart)))

        return result, gens

def _parallel_dict_from_expr_if_gens(exprs, opt):
    """Transform expressions into a multinomial form given generators. """
    k, indices = len(opt.gens), {}

    for i, g in enumerate(opt.gens):
        indices[g] = i

    polys = []

    for expr in exprs:
        poly = {}

        if expr.is_Equality:
            expr = expr.lhs - expr.rhs

        for term in Add.make_args(expr):
            coeff, monom = [], [0]*k

            for factor in Mul.make_args(term):
                if not _not_a_coeff(factor) and factor.is_Number:
                    coeff.append(factor)
                else:
                    try:
                        if opt.series is False:
                            base, exp = decompose_power(factor)

                            if exp < 0:
                                exp, base = -exp, Pow(base, -S.One)
                        else:
                            base, exp = decompose_power_rat(factor)

                        monom[indices[base]] = exp
                    except KeyError:
                        if not factor.free_symbols.intersection(opt.gens):
                            coeff.append(factor)
                        else:
                            raise PolynomialError("%s contains an element of "
                                                  "the set of generators." % factor)

            monom = tuple(monom)

            if monom in poly:
                poly[monom] += Mul(*coeff)
            else:
                poly[monom] = Mul(*coeff)

        polys.append(poly)

    return polys, opt.gens

def _sqrt_match(p):
    """Return [a, b, r] for p.match(a + b*sqrt(r)) where, in addition to
    matching, sqrt(r) also has then maximal sqrt_depth among addends of p.

    Examples
    ========
    >>> from sympy.functions.elementary.miscellaneous import sqrt
    >>> from sympy.simplify.sqrtdenest import _sqrt_match
    >>> _sqrt_match(1 + sqrt(2) + sqrt(2)*sqrt(3) +  2*sqrt(1+sqrt(5)))
    [1 + sqrt(2) + sqrt(6), 2, 1 + sqrt(5)]
    """

    p = _mexpand(p)
    if p.is_Number:
        res = (p, S.Zero, S.Zero)
    elif p.is_Add:
        pargs = list(p.args)
        # to make the process canonical, the argument is included in the tuple
        # so when the max is selected, it will be the largest arg having a
        # given depth
        v = [(sqrt_depth(x), x, i) for i, x in enumerate(pargs)]
        nmax = max(v)
        if nmax[0] == 0:
            res = []
        else:
            depth, _, i = nmax
            r = pargs.pop(i)
            a = Add._from_args(pargs)
            b = S.One
            if r.is_Mul:
                bv = []
                rv = []
                for x in r.args:
                    if sqrt_depth(x) < depth:
                        bv.append(x)
                    else:
                        rv.append(x)

                b = Mul._from_args(bv)
                r = Mul._from_args(rv)
            res = (a, b, r**2)
    else:
        b, r = p.as_coeff_Mul()
        if is_sqrt(r):
            res = (S.Zero, b, r**2)
        else:
            res = []
    return list(res)

    def __new__(cls, *args, **kwargs):
        if not args:
            return GenericZeroMatrix()

        # This must be removed aggressively in the constructor to avoid
        # TypeErrors from GenericZeroMatrix().shape
        args = filter(lambda i: GenericZeroMatrix() != i, args)
        args = list(map(sympify, args))
        check = kwargs.get('check', False)

        obj = Basic.__new__(cls, *args)
        if check:
            if all(not isinstance(i, MatrixExpr) for i in args):
                return Add.fromiter(args)
            validate(*args)
        return obj

    def eval(cls, arg):
        if arg.is_integer:
            return arg
        if arg.is_imaginary:
            return cls(C.im(arg))*S.ImaginaryUnit

        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
                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 + spart
        if not spart:
            return ipart
        elif spart.is_imaginary:
            return ipart + cls(C.im(spart), evaluate=False)*S.ImaginaryUnit
        else:
            return ipart + cls(spart, evaluate=False)

    def __new__(cls, expr, *symbols, **assumptions):

        expr = sympify(expr).expand()
        if expr is S.NaN:
            return S.NaN

        if symbols:
            symbols = map(sympify, symbols)
            if not all(isinstance(s, Symbol) for s in symbols):
                raise NotImplementedError(
                    'Order at points other than 0 not supported.')
        else:
            symbols = list(expr.free_symbols)

        if expr.is_Order:

            new_symbols = list(expr.variables)
            for s in symbols:
                if s not in new_symbols:
                    new_symbols.append(s)
            if len(new_symbols) == len(expr.variables):
                return expr
            symbols = new_symbols

        elif symbols:

            if expr.is_Add:
                lst = expr.extract_leading_order(*symbols)
                expr = Add(*[f.expr for (e, f) in lst])
            elif expr:
                if len(symbols) > 1 or expr.is_commutative is False:
                    # TODO
                    # We cannot use compute_leading_term because that only
                    # works in one symbol.
                    expr = expr.as_leading_term(*symbols)
                else:
                    expr = expr.compute_leading_term(symbols[0])
                terms = expr.as_coeff_mul(*symbols)[1]
                s = set(symbols)
                expr = Mul(*[t for t in terms if s & t.free_symbols])

        if expr is S.Zero:
            return expr
        elif not expr.has(*symbols):
            expr = S.One

        # create Order instance:
        symbols.sort(key=cmp_to_key(Basic.compare))
        obj = Expr.__new__(cls, expr, *symbols, **assumptions)

        return obj

def _parallel_dict_from_expr_if_gens(exprs, opt):
    """Transform expressions into a multinomial form given generators. """
    k, indices = len(opt.gens), {}

    for i, g in enumerate(opt.gens):
        indices[g] = i

    polys = []

    for expr in exprs:
        poly = {}

        for term in Add.make_args(expr):
            coeff, monom = [], [0]*k

            for factor in Mul.make_args(term):
                if factor.is_Number:
                    coeff.append(factor)
                else:
                    try:
                        base, exp = decompose_power(factor)

                        if exp < 0:
                            exp, base = -exp, Pow(base, -S.One)

                        monom[indices[base]] = exp
                    except KeyError:
                        if not factor.has(*opt.gens):
                            coeff.append(factor)
                        else:
                            raise PolynomialError("%s contains an element of the generators set" % factor)

            monom = tuple(monom)

            if monom in poly:
                poly[monom] += Mul(*coeff)
            else:
                poly[monom] = Mul(*coeff)

        polys.append(poly)

    return polys, opt.gens