def _separatevars_dict(expr, *symbols):
    if symbols:
        assert all((t.is_Atom for t in symbols)), "symbols must be Atoms."
    ret = dict(((i,sympify(1)) for i in symbols))
    ret['coeff'] = sympify(1)
    if expr.is_Mul:
        for i in expr.args:
            expsym = i.atoms(Symbol)
            if len(set(symbols).intersection(expsym)) > 1:
                return None
            if len(set(symbols).intersection(expsym)) == 0:
                # There are no symbols, so it is part of the coefficient
                ret['coeff'] *= i
            else:
                ret[expsym.pop()] *= i
    else:
        expsym = expr.atoms(Symbol)
        if len(set(symbols).intersection(expsym)) > 1:
            return None
        if len(set(symbols).intersection(expsym)) == 0:
            # There are no symbols, so it is part of the coefficient
            ret['coeff'] *= expr
        else:
            ret[expsym.pop()] *= expr

    return ret

def bisect(f, a, b, tol):
    """
    Implements bisection. This function is used in RootOf.eval_rational() and
    it needs to be robust.

    Examples
    ========

    >>> from sympy import S
    >>> from sympy.polys.rootoftools import bisect
    >>> bisect(lambda x: x**2-1, -10, 0, S(1)/10**2)
    -1025/1024
    >>> bisect(lambda x: x**2-1, -10, 0, S(1)/10**4)
    -131075/131072

    """
    a = sympify(a)
    b = sympify(b)
    fa = f(a)
    fb = f(b)
    if fa * fb >= 0:
        raise ValueError("bisect: f(a) and f(b) must have opposite signs")
    while (b - a > tol):
        c = (a + b)/2
        fc = f(c)
        if (fc == 0):
            return c  # We need to make sure f(c) is not zero below
        if (fa * fc < 0):
            b = c
            fb = fc
        else:
            a = c
            fa = fc
    return (a + b)/2

 def _eval_expand_log(self, deep=True, **hints):
     if deep:
         arg = self.args[0].expand(deep=deep, **hints)
     else:
         arg = self.args[0]
     if arg.is_Mul:
         expr = sympify(0)
         nonpos = sympify(1)
         for x in arg.args:
             if deep:
                 x = x.expand(deep=deep, **hints)
             if x.is_positive:
                 expr += self.func(x)._eval_expand_log(deep=deep, **hints)
             else:
                 nonpos *= x
         return expr + log(nonpos)
     elif arg.is_Pow:
         if arg.exp.is_real and arg.base.is_positive:
             if deep:
                 b = arg.base.expand(deep=deep, **hints)
                 e = arg.exp.expand(deep=deep, **hints)
             else:
                 b = arg.base
                 e = arg.exp
             return e * self.func(b)._eval_expand_log(deep=deep,\
             **hints)
     return self.func(arg)

    def __getitem__(self, key):
        if not isinstance(key, tuple) and isinstance(key, slice):
            from sympy.matrices.expressions.slice import MatrixSlice

            return MatrixSlice(self, key, (0, None, 1))
        if isinstance(key, tuple) and len(key) == 2:
            i, j = key
            if isinstance(i, slice) or isinstance(j, slice):
                from sympy.matrices.expressions.slice import MatrixSlice

                return MatrixSlice(self, i, j)
            i, j = sympify(i), sympify(j)
            if self.valid_index(i, j) != False:
                return self._entry(i, j)
            else:
                raise IndexError("Invalid indices (%s, %s)" % (i, j))
        elif isinstance(key, (int, Integer)):
            # row-wise decomposition of matrix
            rows, cols = self.shape
            if not (isinstance(rows, Integer) and isinstance(cols, Integer)):
                raise IndexError("Single index only supported for " "non-symbolic matrix shapes.")
            key = sympify(key)
            i = key // cols
            j = key % cols
            if self.valid_index(i, j) != False:
                return self._entry(i, j)
            else:
                raise IndexError("Invalid index %s" % key)
        elif isinstance(key, (Symbol, Expr)):
            raise IndexError("Single index only supported for " "non-symbolic indices.")
        raise IndexError("Invalid index, wanted %s[i,j]" % self)

def real_root(arg, n=None):
    """Return the real nth-root of arg if possible. If n is omitted then
    all instances of (-n)**(1/odd) will be changed to -n**(1/odd); this
    will only create a real root of a principle root -- the presence of
    other factors may cause the result to not be real.

    Examples
    ========

    >>> from sympy import root, real_root, Rational
    >>> from sympy.abc import x, n

    >>> real_root(-8, 3)
    -2
    >>> root(-8, 3)
    2*(-1)**(1/3)
    >>> real_root(_)
    -2

    If one creates a non-principle root and applies real_root, the
    result will not be real (so use with caution):

    >>> root(-8, 3, 2)
    -2*(-1)**(2/3)
    >>> real_root(_)
    -2*(-1)**(2/3)


    See Also
    ========

    sympy.polys.rootoftools.RootOf
    sympy.core.power.integer_nthroot
    root, sqrt
    """
    from sympy import im, Piecewise
    if n is not None:
        try:
            n = as_int(n)
            arg = sympify(arg)
            if arg.is_positive or arg.is_negative:
                rv = root(arg, n)
            else:
                raise ValueError
        except ValueError:
            return root(arg, n)*Piecewise(
                (S.One, ~Equality(im(arg), 0)),
                (Pow(S.NegativeOne, S.One/n)**(2*floor(n/2)), And(
                    Equality(n % 2, 1),
                    arg < 0)),
                (S.One, True))
    else:
        rv = sympify(arg)
    n1pow = Transform(lambda x: -(-x.base)**x.exp,
                      lambda x:
                      x.is_Pow and
                      x.base.is_negative and
                      x.exp.is_Rational and
                      x.exp.p == 1 and x.exp.q % 2)
    return rv.xreplace(n1pow)

    def eval(cls, A):
        from .customized_commands import MatrixAsVector
        from numpy import array
        from numpy.linalg import eig
        from sympy.core import sympify
        
        # convert A to a numpy array of type float64
        try:
            A_array = array(A.tolist(), dtype='float64')
        except (TypeError, AttributeError):
            raise ValueError("Argument to eigenvects_tuple must be a matrix of numerical entries")

        [evals,evects] = eig(A_array)

        eigtuplelist = []
        for i in range(evals.size):
            eigtuplelist.append([sympify(evals[i]),
                                 sympify(evects[:,i].tolist())])


        eigtuplelist.sort(key=lambda w: customized_sort_key(w[0]))

        eiglist=[]
        for t in eigtuplelist:
            eiglist.append(MatrixAsVector(t[1]))
        return TupleNoParen(*eiglist)

    def __new__(
        cls, center=None, hradius=None, vradius=None, eccentricity=None,
            **kwargs):
        hradius = sympify(hradius)
        vradius = sympify(vradius)

        eccentricity = sympify(eccentricity)

        if center is None:
            center = Point(0, 0)
        else:
            center = Point(center, dim=2)

        if len(center) != 2:
            raise ValueError('The center of "{0}" must be a two dimensional point'.format(cls))

        if len(list(filter(None, (hradius, vradius, eccentricity)))) != 2:
            raise ValueError('Exactly two arguments of "hradius", '
                '"vradius", and "eccentricity" must not be None."')

        if eccentricity is not None:
            if hradius is None:
                hradius = vradius / sqrt(1 - eccentricity**2)
            elif vradius is None:
                vradius = hradius * sqrt(1 - eccentricity**2)

        if hradius == vradius:
            return Circle(center, hradius, **kwargs)

        return GeometryEntity.__new__(cls, center, hradius, vradius, **kwargs)

 def _eval_Eq(self, other):
     # CRootOf represents a Root, so if other is that root, it should set
     # the expression to zero *and* it should be in the interval of the
     # CRootOf instance. It must also be a number that agrees with the
     # is_real value of the CRootOf instance.
     if type(self) == type(other):
         return sympify(self == other)
     if not (other.is_number and not other.has(AppliedUndef)):
         return S.false
     if not other.is_finite:
         return S.false
     z = self.expr.subs(self.expr.free_symbols.pop(), other).is_zero
     if z is False:  # all roots will make z True but we don't know
                     # whether this is the right root if z is True
         return S.false
     o = other.is_real, other.is_imaginary
     s = self.is_real, self.is_imaginary
     assert None not in s  # this is part of initial refinement
     if o != s and None not in o:
         return S.false
     re, im = other.as_real_imag()
     if self.is_real:
         if im:
             return S.false
         i = self._get_interval()
         a, b = [Rational(str(_)) for _ in (i.a, i.b)]
         return sympify(a <= other and other <= b)
     i = self._get_interval()
     r1, r2, i1, i2 = [Rational(str(j)) for j in (
         i.ax, i.bx, i.ay, i.by)]
     return sympify((
         r1 <= re and re <= r2) and (
         i1 <= im and im <= i2))

    def eval(cls, x, k):
        x = sympify(x)
        k = sympify(k)

        if x is S.NaN:
            return S.NaN
        elif k.is_Integer:
            if k is S.NaN:
                return S.NaN
            elif k is S.Zero:
                return S.One
            else:
                if k.is_positive:
                    if x is S.Infinity:
                        return S.Infinity
                    elif x is S.NegativeInfinity:
                        if k.is_odd:
                            return S.NegativeInfinity
                        else:
                            return S.Infinity
                    else:
                        return reduce(lambda r, i: r*(x - i), xrange(0, int(k)), 1)
                else:
                    if x is S.Infinity:
                        return S.Infinity
                    elif x is S.NegativeInfinity:
                        return S.Infinity
                    else:
                        return 1/reduce(lambda r, i: r*(x + i), xrange(1, abs(int(k)) + 1), 1)

 def __new__(cls, radius=1, center=[0,0,0], direction=[0,0,1], closed=False, **kwargs):
     """
     >>> from sympy import *
     >>> from symplus.strplus import init_mprinting
     >>> init_mprinting()
     >>> InfiniteCylinder()
     InfiniteCylinder(1, [0 0 0]', [0 0 1]', False)
     >>> InfiniteCylinder(2, [0,0,0], [0,1,1])
     InfiniteCylinder(2, [0 0 0]', [0 -sqrt(2)/2 -sqrt(2)/2]', False)
     >>> InfiniteCylinder().contains((1,1,1))
     False
     >>> InfiniteCylinder(2, [0,0,0], [0,1,1]).contains((1,1,1))
     True
     """
     normalization = kwargs.pop("normalization", True)
     radius = sympify(abs(radius))
     direction = Mat(direction)
     if normalization:
         if norm(direction) == 0:
             raise ValueError
         direction = simplify(normalize(direction))
     direction = max(direction, -direction, key=hash)
     center = Mat(center)
     if normalization:
         center = simplify(center - project(center, direction))
     closed = sympify(bool(closed))
     return Basic.__new__(cls, radius, center, direction, closed)

def jn(n, z):
    """
    Spherical Bessel function of the first kind.

    Examples:

        >>> from sympy import Symbol, jn, sin, cos
        >>> z = Symbol("z")
        >>> print jn(0, z)
        sin(z)/z
        >>> jn(1, z) == sin(z)/z**2 - cos(z)/z
        True
        >>> jn(3, z) ==(1/z - 15/z**3)*cos(z) + (15/z**4 - 6/z**2)*sin(z)
        True

    The spherical Bessel functions are calculated using the formula:

    jn(n, z) == fn(n, z) * sin(z) + (-1)**(n+1) * fn(-n-1, z) * cos(z)

    where fn(n, z) are the coefficients, see fn()'s sourcecode for more
    information.
    """

    n = sympify(n)
    z = sympify(z)
    return fn(n, z) * sin(z) + (-1)**(n+1) * fn(-n-1, z) * cos(z)

    def eval(cls, arg, base=None):
        from sympy import unpolarify
        if base is not None:
            base = sympify(base)

            if arg.is_positive and arg.is_Integer and \
               base.is_positive and base.is_Integer:
                base = int(base)
                arg = int(arg)
                n = multiplicity(base, arg)
                return S(n) + log(arg // base ** n) / log(base)
            if base is not S.Exp1:
                return cls(arg)/cls(base)
            else:
                return cls(arg)

        arg = sympify(arg)

        if arg.is_Number:
            if arg is S.Zero:
                return S.ComplexInfinity
            elif arg is S.One:
                return S.Zero
            elif arg is S.Infinity:
                return S.Infinity
            elif arg is S.NegativeInfinity:
                return S.Infinity
            elif arg is S.NaN:
                return S.NaN
            elif arg.is_negative:
                return S.Pi * S.ImaginaryUnit + cls(-arg)
            elif arg.is_Rational:
                if arg.q != 1:
                    return cls(arg.p) - cls(arg.q)
                # remove perfect powers automatically
                p = perfect_power(int(arg))
                if p is not False:
                    return p[1]*cls(p[0])
        elif arg is S.ComplexInfinity:
            return S.ComplexInfinity
        elif arg is S.Exp1:
            return S.One
        elif arg.func is exp and arg.args[0].is_real:
            return arg.args[0]
        elif arg.func is exp_polar:
            return unpolarify(arg.exp)
        #don't autoexpand Pow or Mul (see the issue 252):
        elif not arg.is_Add:
            coeff = arg.as_coefficient(S.ImaginaryUnit)

            if coeff is not None:
                if coeff is S.Infinity:
                    return S.Infinity
                elif coeff is S.NegativeInfinity:
                    return S.Infinity
                elif coeff.is_Rational:
                    if coeff.is_nonnegative:
                        return S.Pi * S.ImaginaryUnit * S.Half + cls(coeff)
                    else:
                        return -S.Pi * S.ImaginaryUnit * S.Half + cls(-coeff)

 def __new__(cls, function, limits):
     fun = sympify(function)
     if not ordered_iter(fun) or len(fun) != 2:
         raise ValueError("Function argument should be (x(t), y(t)) but got %s" % str(function))
     if not ordered_iter(limits) or len(limits) != 3:
         raise ValueError("Limit argument should be (t, tmin, tmax) but got %s" % str(limits))
     return GeometryEntity.__new__(cls, tuple(sympify(fun)), tuple(sympify(limits)))

    def __new__(cls, center=None, hradius=None, vradius=None, eccentricity=None,
                **kwargs):
        hradius = sympify(hradius)
        vradius = sympify(vradius)
        eccentricity = sympify(eccentricity)

        if len(filter(None, (hradius, vradius, eccentricity))) != 2:
            raise ValueError, 'Exactly two arguments between "hradius", '\
                '"vradius", and "eccentricity" must be not None."'

        if eccentricity is not None:
            if hradius is None:
                hradius = vradius / sqrt(1 - eccentricity**2)
            elif vradius is None:
                vradius = hradius * sqrt(1 - eccentricity**2)
        else:
            if hradius is None and vradius is None:
                raise ValueError("At least two arguments between hradius, "
                    "vradius and eccentricity must not be none.")

        if center is None:
            center = Point(0, 0)

        if not isinstance(center, Point):
            raise TypeError("center must be a Point")

        if hradius == vradius:
            return Circle(center, hradius, **kwargs)
        return GeometryEntity.__new__(cls, center, hradius, vradius, **kwargs)

 def __new__(cls, e, z, z0, dir="+"):
     e = sympify(e)
     z = sympify(z)
     z0 = sympify(z0)
     obj = Expr.__new__(cls)
     obj._args = (e, z, z0, dir)
     return obj

def swinnerton_dyer_poly(n, x=None, **args):
    """Generates n-th Swinnerton-Dyer polynomial in `x`.  """
    from numberfields import minimal_polynomial
    if n <= 0:
        raise ValueError(
            "can't generate Swinnerton-Dyer polynomial of order %s" % n)

    if x is not None:
        sympify(x)
    else:
        x = Dummy('x')

    if n > 3:
        p = 2
        a = [sqrt(2)]
        for i in xrange(2, n + 1):
            p = nextprime(p)
            a.append(sqrt(p))
        return minimal_polynomial(Add(*a), x, polys=args.get('polys', False))

    if n == 1:
        ex = x**2 - 2
    elif n == 2:
        ex = x**4 - 10*x**2 + 1
    elif n == 3:
        ex = x**8 - 40*x**6 + 352*x**4 - 960*x**2 + 576
    if not args.get('polys', False):
        return ex
    else:
        return PurePoly(ex, x)

    def __new__(cls, f, x, index=None, radicals=False, expand=True):
        """ Construct an indexed complex root of a polynomial.

        See ``rootof`` for the parameters.

        The default value of ``radicals`` is ``False`` to satisfy
        ``eval(srepr(expr) == expr``.
        """
        x = sympify(x)

        if index is None and x.is_Integer:
            x, index = None, x
        else:
            index = sympify(index)

        if index is not None and index.is_Integer:
            index = int(index)
        else:
            raise ValueError("expected an integer root index, got %s" % index)

        poly = PurePoly(f, x, greedy=False, expand=expand)

        if not poly.is_univariate:
            raise PolynomialError("only univariate polynomials are allowed")

        if not poly.gen.is_Symbol:
            # PurePoly(sin(x) + 1) == PurePoly(x + 1) but the roots of
            # x for each are not the same: issue 8617
            raise PolynomialError("generator must be a Symbol")

        degree = poly.degree()

        if degree <= 0:
            raise PolynomialError("can't construct CRootOf object for %s" % f)

        if index < -degree or index >= degree:
            raise IndexError("root index out of [%d, %d] range, got %d" %
                             (-degree, degree - 1, index))
        elif index < 0:
            index += degree

        dom = poly.get_domain()

        if not dom.is_Exact:
            poly = poly.to_exact()

        roots = cls._roots_trivial(poly, radicals)

        if roots is not None:
            return roots[index]

        coeff, poly = preprocess_roots(poly)
        dom = poly.get_domain()

        if not dom.is_ZZ:
            raise NotImplementedError("CRootOf is not supported over %s" % dom)

        root = cls._indexed_root(poly, index)
        return coeff * cls._postprocess_root(root, radicals)

    def eval(cls, arg, base=None):
        if base is not None:
            base = sympify(base)

            if arg.is_positive and arg.is_Integer and \
               base.is_positive and base.is_Integer:
                base = int(base)
                arg = int(arg)
                n = multiplicity(base, arg)
                return S(n) + log(arg // base ** n) / log(base)
            if base is not S.Exp1:
                return cls(arg)/cls(base)
            else:
                return cls(arg)

        arg = sympify(arg)

        if arg.is_Number:
            if arg is S.Zero:
                return S.ComplexInfinity
            elif arg is S.One:
                return S.Zero
            elif arg is S.Infinity:
                return S.Infinity
            elif arg is S.NegativeInfinity:
                return S.Infinity
            elif arg is S.NaN:
                return S.NaN
            elif arg.is_negative:
                return S.Pi * S.ImaginaryUnit + cls(-arg)
        elif arg is S.ComplexInfinity:
            return S.ComplexInfinity
        elif arg is S.Exp1:
            return S.One
        #this doesn't work due to caching: :(
        #elif arg.func is exp and arg.args[0].is_real:
        #using this one instead:
        elif arg.func is exp and arg.args[0].is_real:
            return arg.args[0]
        #this shouldn't happen automatically (see the issue 252):
        #elif arg.is_Pow:
        #    if arg.exp.is_Number or arg.exp.is_NumberSymbol or \
        #        arg.exp.is_number:
        #        return arg.exp * self(arg.base)
        #elif arg.is_Mul and arg.is_real:
        #    return Add(*[self(a) for a in arg])
        elif not arg.is_Add:
            coeff = arg.as_coefficient(S.ImaginaryUnit)

            if coeff is not None:
                if coeff is S.Infinity:
                    return S.Infinity
                elif coeff is S.NegativeInfinity:
                    return S.Infinity
                elif coeff.is_Rational:
                    if coeff.is_nonnegative:
                        return S.Pi * S.ImaginaryUnit * S.Half + cls(coeff)
                    else:
                        return -S.Pi * S.ImaginaryUnit * S.Half + cls(-coeff)

    def eval(cls, arg, k=0):
        """
        Returns a simplified form or a value of DiracDelta depending on the
        argument passed by the DiracDelta object.

        The ``eval()`` method is automatically called when the ``DiracDelta`` class
        is about to be instantiated and it returns either some simplified instance
        or the unevaluated instance depending on the argument passed. In other words,
        ``eval()`` method is not needed to be called explicitly, it is being called
        and evaluated once the object is called.

        Examples
        ========

        >>> from sympy import DiracDelta, S, Subs
        >>> from sympy.abc import x

        >>> DiracDelta(x)
        DiracDelta(x)

        >>> DiracDelta(x,1)
        DiracDelta(x, 1)

        >>> DiracDelta(1)
        0

        >>> DiracDelta(5,1)
        0

        >>> DiracDelta(0)
        DiracDelta(0)

        >>> DiracDelta(-1)
        0

        >>> DiracDelta(S.NaN)
        nan

        >>> DiracDelta(x).eval(1)
        0

        >>> DiracDelta(x - 100).subs(x, 5)
        0

        >>> DiracDelta(x - 100).subs(x, 100)
        DiracDelta(0)

        """
        k = sympify(k)
        if not k.is_Integer or k.is_negative:
            raise ValueError("Error: the second argument of DiracDelta must be \
            a non-negative integer, %s given instead." % (k,))
        arg = sympify(arg)
        if arg is S.NaN:
            return S.NaN
        if arg.is_nonzero:
            return S.Zero
        if fuzzy_not(im(arg).is_zero):
            raise ValueError("Function defined only for Real Values. Complex part: %s  found in %s ." % (repr(im(arg)), repr(arg)) )

 def __new__(cls, expr1, expr2):
     expr1 = sympify(expr1)
     expr2 = sympify(expr2)
     expr1, expr2 = sorted([expr1, expr2], key=default_sort_key)
     obj = Expr.__new__(cls, expr1, expr2)
     obj._expr1 = expr1
     obj._expr2 = expr2
     return obj