def test_roots_binomial():
    assert roots_binomial(Poly(5*x, x)) == [0]
    assert roots_binomial(Poly(5*x**4, x)) == [0, 0, 0, 0]
    assert roots_binomial(Poly(5*x + 2, x)) == [-Rational(2, 5)]

    A = 10**Rational(3, 4)/10

    assert roots_binomial(Poly(5*x**4 + 2, x)) == \
        [-A - A*I, -A + A*I, A - A*I, A + A*I]

    a1 = Symbol('a1', nonnegative=True)
    b1 = Symbol('b1', nonnegative=True)

    r0 = roots_quadratic(Poly(a1*x**2 + b1, x))
    r1 = roots_binomial(Poly(a1*x**2 + b1, x))

    assert powsimp(r0[0]) == powsimp(r1[0])
    assert powsimp(r0[1]) == powsimp(r1[1])
    for a, b, s, n in cartes((1, 2), (1, 2), (-1, 1), (2, 3, 4, 5)):
        if a == b and a != 1:  # a == b == 1 is sufficient
            continue
        p = Poly(a*x**n + s*b)
        ans = roots_binomial(p)
        assert ans == _nsort(ans)

    # issue 8813
    assert roots(Poly(2*x**3 - 16*y**3, x)) == {
        2*y*(-S(1)/2 - sqrt(3)*I/2): 1,
        2*y: 1,
        2*y*(-S(1)/2 + sqrt(3)*I/2): 1}

def generate_Sx_Sy_Sz_operators(atom_list, Sa, Sb, Sn):
    """Generates Sx, Sy and Sz operators"""
    Sx_list = []; Sy_list = []; Sz_list = []
    N = len(atom_list)

    loc_vect = spm.Matrix([Sa,Sb,Sn])
    loc_vect = loc_vect.reshape(3,1)

    for i in range(N):
        rotmat = sp.Matrix(atom_list[i].spinRmatrix)
        glo_vect = rotmat * loc_vect

        Sx = sp.powsimp(glo_vect[0].expand())
        Sy = sp.powsimp(glo_vect[1].expand())
        Sz = sp.powsimp(glo_vect[2].expand())

        Sx_list.append(Sx)
        Sy_list.append(Sy)
        Sz_list.append(Sz)
            
    Sx_list.append(KAPXHAT_SYM)
    Sy_list.append(KAPYHAT_SYM)
    Sz_list.append(KAPZHAT_SYM)
    print "Operators Generated: Sx, Sy, Sz"
    return (Sx_list,Sy_list,Sz_list)

def generate_Sx_Sy_Sz_operators(atom_list, Sa, Sb, Sn):
    """Generates Sx, Sy and Sz operators"""
    Sx_list = []; Sy_list = []; Sz_list = []
    N = len(atom_list)
    S = sp.Symbol('S', commutative = True)
    loc_vect = spm.Matrix([Sa,Sb,Sn])
    loc_vect = loc_vect.reshape(3,1)

    for i in range(N):
        rotmat = sp.Matrix(atom_list[i].spinRmatrix)
        glo_vect = rotmat * loc_vect

        Sx = sp.powsimp(glo_vect[0].expand())
        Sy = sp.powsimp(glo_vect[1].expand())
        Sz = sp.powsimp(glo_vect[2].expand())

        Sx_list.append(Sx)
        Sy_list.append(Sy)
        Sz_list.append(Sz)
        
    #Unit vector markers
    kapxhat = sp.Symbol('kapxhat',real = True)#spm.Matrix([1,0,0])#
    kapyhat = sp.Symbol('kapyhat',real = True)#spm.Matrix([0,1,0])#
    kapzhat = sp.Symbol('kapzhat',real = True)#spm.Matrix([0,0,1])#
    
    Sx_list.append(kapxhat)
    Sy_list.append(kapyhat)
    Sz_list.append(kapzhat)
    print "Operators Generated: Sx, Sy, Sz"
    return (Sx_list,Sy_list,Sz_list)

def test_issue_3268():
    z = -5*sqrt(2)/(2*sqrt(2*sqrt(29) + 29)) + sqrt(-sqrt(29)/29 + S(1)/2)
    assert Mul(*[powsimp(a) for a in Mul.make_args(z.normal())]) == 0
    assert powsimp(z.normal()) == 0
    assert simplify(z) == 0
    assert powsimp(sqrt(2 + sqrt(3))*sqrt(2 - sqrt(3)) + 1) == 2
    assert powsimp(z) != 0

def test_issue_from_PR1599():
    n1, n2, n3, n4 = symbols('n1 n2 n3 n4', negative=True)
    assert (powsimp(sqrt(n1)*sqrt(n2)*sqrt(n3)) ==
        -I*sqrt(-n1)*sqrt(-n2)*sqrt(-n3))
    assert (powsimp(root(n1, 3)*root(n2, 3)*root(n3, 3)*root(n4, 3)) ==
        -(-1)**(S(1)/3)*
        (-n1)**(S(1)/3)*(-n2)**(S(1)/3)*(-n3)**(S(1)/3)*(-n4)**(S(1)/3))

def test_powsimp_negated_base():
    assert powsimp((-x + y)/sqrt(x - y)) == -sqrt(x - y)
    assert powsimp((-x + y)*(-z + y)/sqrt(x - y)/sqrt(z - y)) == sqrt(x - y)*sqrt(z - y)
    p = symbols('p', positive=True)
    assert powsimp((-p)**a/p**a) == (-1)**a
    n = symbols('n', negative=True)
    assert powsimp((-n)**a/n**a) == (-1)**a
    # if x is 0 then the lhs is 0**a*oo**a which is not (-1)**a
    assert powsimp((-x)**a/x**a) != (-1)**a

def test_issue_10195():
    a = Symbol('a', integer=True)
    l = Symbol('l', even=True, nonzero=True)
    n = Symbol('n', odd=True)
    e_x = (-1)**(n/2 - Rational(1, 2)) - (-1)**(3*n/2 - Rational(1, 2))
    assert powsimp((-1)**(l/2)) == I**l
    assert powsimp((-1)**(n/2)) == I**n
    assert powsimp((-1)**(3*n/2)) == -I**n
    assert powsimp(e_x) == (-1)**(n/2 - Rational(1, 2)) + (-1)**(3*n/2 +
            Rational(1,2))
    assert powsimp((-1)**(3*a/2)) == (-I)**a

def test_powsimp():
    x,y,n = symbols('xyn')
    f = Function('f')
    assert powsimp( 4**x * 2**(-x) * 2**(-x) ) == 1
    assert powsimp( (-4)**x * (-2)**(-x) * 2**(-x) ) == 1

    assert powsimp( f(4**x * 2**(-x) * 2**(-x)) )   == f(4**x * 2**(-x) * 2**(-x))
    assert powsimp( f(4**x * 2**(-x) * 2**(-x)), deep = True )  == f(1)
    x,y = symbols('xy', nonnegative=True)
    n = Symbol('n', real=True)
    assert powsimp( y**n * (y/x)**(-n) ) == x**n

def list_mult(lista, listb):
    "Defines a way to multiply two lists of the same length"
    if len(lista) != len(listb):
        print "lists not same length"
        return []
    else:
        temp = []
        for i in range(len(lista)):
            if isinstance(lista[i], int):
                temp.append(sp.powsimp((lista[i] * listb[i]).expand(deep=False)))
            elif isinstance(listb[i], int):
                temp.append(sp.powsimp((lista[i] * listb[i]).expand(deep=False)))
            else:
                temp.append(sp.powsimp((lista[i] * listb[i]).expand(deep=False)))
        return temp

 def _contains(self, expr):
     from sympy import powsimp, limit
     if expr is S.Zero:
         return True
     if expr is S.NaN:
         return False
     if expr.is_Order:
         if self.variables and expr.variables:
             common_symbols = tuple([s for s in self.variables if s in expr.variables])
         elif self.variables:
             common_symbols = self.variables
         else:
             common_symbols = expr.variables
         if not common_symbols:
             if not (self.variables or expr.variables): # O(1),O(1)
                 return True
             return None
         r = None
         for s in common_symbols:
             l = limit(powsimp(self.expr/expr.expr, deep=True,\
             combine='exp'), s, 0) != 0
             if r is None:
                 r = l
             else:
                 if r != l:
                     return
         return r
     obj = Order(expr, *self.variables)
     return self.contains(obj)

 def contains(self, expr):
     """
     Return True if expr belongs to Order(self.expr, *self.variables).
     Return False if self belongs to expr.
     Return None if the inclusion relation cannot be determined
     (e.g. when self and expr have different symbols).
     """
     from sympy import powsimp, limit
     if expr is S.Zero:
         return True
     if expr is S.NaN:
         return False
     if expr.is_Order:
         if self.variables and expr.variables:
             common_symbols = tuple([s for s in self.variables if s in expr.variables])
         elif self.variables:
             common_symbols = self.variables
         else:
             common_symbols = expr.variables
         if not common_symbols:
             if not (self.variables or expr.variables): # O(1),O(1)
                 return True
             return None
         r = None
         for s in common_symbols:
             l = limit(powsimp(self.expr/expr.expr, deep=True,\
             combine='exp'), s, 0) != 0
             if r is None:
                 r = l
             else:
                 if r != l:
                     return
         return r
     obj = Order(expr, *self.variables)
     return self.contains(obj)

 def contains(self, expr):
     """
     Return True if expr belongs to Order(self.expr, \*self.variables).
     Return False if self belongs to expr.
     Return None if the inclusion relation cannot be determined
     (e.g. when self and expr have different symbols).
     """
     from sympy import powsimp, PoleError
     if expr is S.Zero:
         return True
     if expr is S.NaN:
         return False
     if expr.is_Order:
         if not all(p == expr.point[0] for p in expr.point) and \
            not all(p == self.point[0] for p in self.point):
             raise NotImplementedError('Order at points other than 0 '
                 'or oo not supported, got %s as a point.' % point)
         else:
             # self and/or expr is O(1):
             if any(not p for p in [expr.point, self.point]):
                 point = self.point + expr.point
                 if point:
                     point = point[0]
                 else:
                     point = S.Zero
             else:
                 point = self.point[0]
         if expr.expr == self.expr:
             # O(1) + O(1), O(1) + O(1, x), etc.
             return all([x in self.args[1:] for x in expr.args[1:]])
         if expr.expr.is_Add:
             return all([self.contains(x) for x in expr.expr.args])
         if self.expr.is_Add:
             return any([self.func(x, *self.args[1:]).contains(expr)
                         for x in self.expr.args])
         if self.variables and expr.variables:
             common_symbols = tuple(
                 [s for s in self.variables if s in expr.variables])
         elif self.variables:
             common_symbols = self.variables
         else:
             common_symbols = expr.variables
         if not common_symbols:
             return None
         r = None
         ratio = self.expr/expr.expr
         ratio = powsimp(ratio, deep=True, combine='exp')
         for s in common_symbols:
             try:
                 l = ratio.limit(s, point) != 0
             except PoleError:
                 l = None
             if r is None:
                 r = l
             else:
                 if r != l:
                     return
         return r
     obj = self.func(expr, *self.args[1:])
     return self.contains(obj)

def _simplify(expr):
    u"""Simplifie une expression.

    Alias de simplify (sympy 0.6.4).
    Mais simplify n'est pas garanti d'être stable dans le temps.
    (Cf. simplify.__doc__)."""
    return together(expand(Poly.cancel(powsimp(expr))))

    def as_leading_term(self, *symbols):
        """
        Returns the leading term.

        Example:

        >>> from sympy.abc import x
        >>> (1+x+x**2).as_leading_term(x)
        1
        >>> (1/x**2+x+x**2).as_leading_term(x)
        x**(-2)

        Note:

        self is assumed to be the result returned by Basic.series().
        """
        from sympy import powsimp
        if len(symbols)>1:
            c = self
            for x in symbols:
                c = c.as_leading_term(x)
            return c
        elif not symbols:
            return self
        x = sympify(symbols[0])
        assert x.is_Symbol, `x`
        if not self.has(x):
            return self
        obj = self._eval_as_leading_term(x)
        if obj is not None:
            return powsimp(obj, deep=True, combine='exp')
        raise NotImplementedError('as_leading_term(%s, %s)' % (self, x))

def test_TR2i():
    # just a reminder that ratios of powers only simplify if both
    # numerator and denominator satisfy the condition that each
    # has a positive base or an integer exponent; e.g. the following,
    # at y=-1, x=1/2 gives sqrt(2)*I != -sqrt(2)*I
    assert powsimp(2**x/y**x) != (2/y)**x

    assert TR2i(sin(x)/cos(x)) == tan(x)
    assert TR2i(sin(x)*sin(y)/cos(x)) == tan(x)*sin(y)
    assert TR2i(1/(sin(x)/cos(x))) == 1/tan(x)
    assert TR2i(1/(sin(x)*sin(y)/cos(x))) == 1/tan(x)/sin(y)
    assert TR2i(sin(x)/2/(cos(x) + 1)) == sin(x)/(cos(x) + 1)/2

    assert TR2i(sin(x)/2/(cos(x) + 1), half=True) == tan(x/2)/2
    assert TR2i(sin(1)/(cos(1) + 1), half=True) == tan(S.Half)
    assert TR2i(sin(2)/(cos(2) + 1), half=True) == tan(1)
    assert TR2i(sin(4)/(cos(4) + 1), half=True) == tan(2)
    assert TR2i(sin(5)/(cos(5) + 1), half=True) == tan(5*S.Half)
    assert TR2i((cos(1) + 1)/sin(1), half=True) == 1/tan(S.Half)
    assert TR2i((cos(2) + 1)/sin(2), half=True) == 1/tan(1)
    assert TR2i((cos(4) + 1)/sin(4), half=True) == 1/tan(2)
    assert TR2i((cos(5) + 1)/sin(5), half=True) == 1/tan(5*S.Half)
    assert TR2i((cos(1) + 1)**(-a)*sin(1)**a, half=True) == tan(S.Half)**a
    assert TR2i((cos(2) + 1)**(-a)*sin(2)**a, half=True) == tan(1)**a
    assert TR2i((cos(4) + 1)**(-a)*sin(4)**a, half=True) == (cos(4) + 1)**(-a)*sin(4)**a
    assert TR2i((cos(5) + 1)**(-a)*sin(5)**a, half=True) == (cos(5) + 1)**(-a)*sin(5)**a
    assert TR2i((cos(1) + 1)**a*sin(1)**(-a), half=True) == tan(S.Half)**(-a)
    assert TR2i((cos(2) + 1)**a*sin(2)**(-a), half=True) == tan(1)**(-a)
    assert TR2i((cos(4) + 1)**a*sin(4)**(-a), half=True) == (cos(4) + 1)**a*sin(4)**(-a)
    assert TR2i((cos(5) + 1)**a*sin(5)**(-a), half=True) == (cos(5) + 1)**a*sin(5)**(-a)

    i = symbols('i', integer=True)
    assert TR2i(((cos(5) + 1)**i*sin(5)**(-i)), half=True) == tan(5*S.Half)**(-i)
    assert TR2i(1/((cos(5) + 1)**i*sin(5)**(-i)), half=True) == tan(5*S.Half)**i

def eqdiff(leq,lvars):
    '''Calculate the Jacobian of the system of equations with respect to a set of variables.'''
    from sympy import powsimp
    resp = []
    for eq in leq:
        el = [ powsimp(eq.diff(v)) for v in lvars]
        resp += [el]
    return resp

 def taylor_term(n, x, *previous_terms): # of log(1+x)
     """
     Returns the next term in the Taylor series expansion of log(1+x).
     """
     from sympy import powsimp
     if n<0: return S.Zero
     x = sympify(x)
     if n==0: return x
     if previous_terms:
         p = previous_terms[-1]
         if p is not None:
             return powsimp((-n) * p * x / (n+1), deep=True, combine='exp')
     return (1-2*(n%2)) * x**(n+1)/(n+1)

 def _eval_nseries(self, x, n):
     from sympy import limit, oo, powsimp
     arg = self.args[0]
     arg_series = arg._eval_nseries(x, n=n)
     if arg_series.is_Order:
         return 1 + arg_series
     arg0 = limit(arg_series.removeO(), x, 0)
     if arg0 in [-oo, oo]:
         return self
     t = Dummy("t")
     exp_series = exp(t)._taylor(t, n)
     r = exp(arg0)*exp_series.subs(t, arg_series - arg0)
     r = r.expand()
     return powsimp(r, deep=True, combine='exp')

def test_roots_binomial():
    assert roots_binomial(Poly(5 * x, x)) == [0]
    assert roots_binomial(Poly(5 * x ** 4, x)) == [0, 0, 0, 0]
    assert roots_binomial(Poly(5 * x + 2, x)) == [-Rational(2, 5)]

    A = 10 ** Rational(3, 4) / 10

    assert roots_binomial(Poly(5 * x ** 4 + 2, x)) == [-A - A * I, -A + A * I, A - A * I, A + A * I]

    a1 = Symbol("a1", nonnegative=True)
    b1 = Symbol("b1", nonnegative=True)

    r0 = roots_quadratic(Poly(a1 * x ** 2 + b1, x))
    r1 = roots_binomial(Poly(a1 * x ** 2 + b1, x))

    assert powsimp(r0[0]) == powsimp(r1[0])
    assert powsimp(r0[1]) == powsimp(r1[1])
    for a, b, s, n in cartes((1, 2), (1, 2), (-1, 1), (2, 3, 4, 5)):
        if a == b and a != 1:  # a == b == 1 is sufficient
            continue
        p = Poly(a * x ** n + s * b)
        roots = roots_binomial(p)
        assert roots == _nsort(roots)

def test_action_verbs():
    assert nsimplify((1/(exp(3*pi*x/5)+1))) == (1/(exp(3*pi*x/5)+1)).nsimplify()
    assert ratsimp(1/x + 1/y) == (1/x + 1/y).ratsimp()
    assert trigsimp(log(x), deep=True) == (log(x)).trigsimp(deep = True)
    assert radsimp(1/(2+sqrt(2))) == (1/(2+sqrt(2))).radsimp()
    assert powsimp(x**y*x**z*y**z, combine='all') == (x**y*x**z*y**z).powsimp(combine='all')
    assert simplify(x**y*x**z*y**z) == (x**y*x**z*y**z).simplify()
    assert together(1/x + 1/y) == (1/x + 1/y).together()
    assert separate((x*(y*z)**3)**2) == ((x*(y*z)**3)**2).separate()
    assert collect(a*x**2 + b*x**2 + a*x - b*x + c, x) == (a*x**2 + b*x**2 + a*x - b*x + c).collect(x)
    assert apart(y/(y+2)/(y+1), y) == (y/(y+2)/(y+1)).apart(y)
    assert combsimp(y/(x+2)/(x+1)) == (y/(x+2)/(x+1)).combsimp()
    assert factor(x**2+5*x+6) == (x**2+5*x+6).factor()
    assert refine(sqrt(x**2)) == sqrt(x**2).refine()
    assert cancel((x**2+5*x+6)/(x+2)) == ((x**2+5*x+6)/(x+2)).cancel()