def test_as_real_imag():
    n = pi**1000
    # the special code for working out the real
    # and complex parts of a power with Integer exponent
    # should not run if there is no imaginary part, hence
    # this should not hang
    assert n.as_real_imag() == (n, 0)

    # issue 6261
    x = Symbol('x')
    assert sqrt(x).as_real_imag() == \
        ((re(x)**2 + im(x)**2)**(S(1)/4)*cos(atan2(im(x), re(x))/2),
     (re(x)**2 + im(x)**2)**(S(1)/4)*sin(atan2(im(x), re(x))/2))

    # issue 3853
    a, b = symbols('a,b', real=True)
    assert ((1 + sqrt(a + b*I))/2).as_real_imag() == \
           (
               (a**2 + b**2)**Rational(
                   1, 4)*cos(atan2(b, a)/2)/2 + Rational(1, 2),
               (a**2 + b**2)**Rational(1, 4)*sin(atan2(b, a)/2)/2)

    assert sqrt(a**2).as_real_imag() == (sqrt(a**2), 0)
    i = symbols('i', imaginary=True)
    assert sqrt(i**2).as_real_imag() == (0, abs(i))

def test_f_expand_complex():
    x = Symbol("x", real=True)

    assert f(x).expand(complex=True) == I * im(f(x)) + re(f(x))
    assert exp(x).expand(complex=True) == exp(x)
    assert exp(I * x).expand(complex=True) == cos(x) + I * sin(x)
    assert exp(z).expand(complex=True) == cos(im(z)) * exp(re(z)) + I * sin(im(z)) * exp(re(z))

    def _generic_mul(q1, q2):

        q1 = sympify(q1)
        q2 = sympify(q2)

        # None is a Quaternion:
        if not isinstance(q1, Quaternion) and not isinstance(q2, Quaternion):
            return q1 * q2

        # If q1 is a number or a sympy expression instead of a quaternion
        if not isinstance(q1, Quaternion):
            if q2.real_field:
                if q1.is_complex:
                    return q2 * Quaternion(re(q1), im(q1), 0, 0)
                else:
                    return Mul(q1, q2)
            else:
                return Quaternion(q1 * q2.a, q1 * q2.b, q1 * q2.c, q1 * q2.d)


        # If q2 is a number or a sympy expression instead of a quaternion
        if not isinstance(q2, Quaternion):
            if q1.real_field:
                if q2.is_complex:
                    return q1 * Quaternion(re(q2), im(q2), 0, 0)
                else:
                    return Mul(q1, q2)
            else:
                return Quaternion(q2 * q1.a, q2 * q1.b, q2 * q1.c, q2 * q1.d)

        return Quaternion(-q1.b*q2.b - q1.c*q2.c - q1.d*q2.d + q1.a*q2.a,
                          q1.b*q2.a + q1.c*q2.d - q1.d*q2.c + q1.a*q2.b,
                          -q1.b*q2.d + q1.c*q2.a + q1.d*q2.b + q1.a*q2.c,
                          q1.b*q2.c - q1.c*q2.b + q1.d*q2.a + q1.a * q2.d)

def velocity_field(psi): #takes a symbolic function and returns two lambda functions
#to evaluate the derivatives in both x and y.
   global w
   if velocity_components:
      u = lambdify((x,y), eval(x_velocity), modules='numpy')
      v = lambdify((x,y), eval(y_velocity), modules='numpy')
   else:
      if is_complex_potential:
         print "Complex potential, w(z) given"
         #define u, v symbolically as the imaginary part of the derivatives
         u = lambdify((x, y), sympy.im(psi.diff(y)), modules='numpy')
         v = lambdify((x, y), -sympy.im(psi.diff(x)), modules='numpy')
      else:
         #define u,v as the derivatives 
         print "Stream function, psi given"
         u = sympy.lambdify((x, y), psi.diff(y), 'numpy')
         v = sympy.lambdify((x, y), -psi.diff(x), 'numpy')
   if (branch_cuts): # If it's indicated that there are branch cuts in the mapping,
                      # then we need to return vectorized numpy functions to evaluate
                      # everything numerically, instead of symbolically 
                      # This of course results in a SIGNIFICANT time increase
                      #   (I don't know how to handle more than the primitive root
                      #   (symbolically in Sympy
      return np.vectorize(u),np.vectorize(v)
   else:
       # If there are no branch cuts, then return the symbolic lambda functions (MUCH faster)
      return u,v

def test_re():
    x, y = symbols('x,y')
    a, b = symbols('a,b', real=True)

    r = Symbol('r', real=True)
    i = Symbol('i', imaginary=True)

    assert re(nan) == nan

    assert re(oo) == oo
    assert re(-oo) == -oo

    assert re(0) == 0

    assert re(1) == 1
    assert re(-1) == -1

    assert re(E) == E
    assert re(-E) == -E

    assert re(x) == re(x)
    assert re(x*I) == -im(x)
    assert re(r*I) == 0
    assert re(r) == r
    assert re(i*I) == I * i
    assert re(i) == 0

    assert re(x + y) == re(x + y)
    assert re(x + r) == re(x) + r

    assert re(re(x)) == re(x)

    assert re(2 + I) == 2
    assert re(x + I) == re(x)

    assert re(x + y*I) == re(x) - im(y)
    assert re(x + r*I) == re(x)

    assert re(log(2*I)) == log(2)

    assert re((2 + I)**2).expand(complex=True) == 3

    assert re(conjugate(x)) == re(x)
    assert conjugate(re(x)) == re(x)

    assert re(x).as_real_imag() == (re(x), 0)

    assert re(i*r*x).diff(r) == re(i*x)
    assert re(i*r*x).diff(i) == I*r*im(x)

    assert re(
        sqrt(a + b*I)) == (a**2 + b**2)**Rational(1, 4)*cos(atan2(b, a)/2)
    assert re(a * (2 + b*I)) == 2*a

    assert re((1 + sqrt(a + b*I))/2) == \
        (a**2 + b**2)**Rational(1, 4)*cos(atan2(b, a)/2)/2 + Rational(1, 2)

    assert re(x).rewrite(im) == x - im(x)
    assert (x + re(y)).rewrite(re, im) == x + y - im(y)

def test_reduce_inequalities_multivariate():
    x = Symbol('x')
    y = Symbol('y')

    assert reduce_inequalities([Ge(x**2, 1), Ge(y**2, 1)]) == \
        And(Eq(im(x), 0), Eq(im(y), 0), Or(And(Le(1, re(x)), Lt(re(x), oo)),
                                           And(Le(re(x), -1), Lt(-oo, re(x)))),
            Or(And(Le(1, re(y)), Lt(re(y), oo)), And(Le(re(y), -1), Lt(-oo, re(y)))))

def test_derivatives_issue1658():
    x = Symbol('x')
    f = Function('f')
    assert re(f(x)).diff(x) == re(f(x).diff(x))
    assert im(f(x)).diff(x) == im(f(x).diff(x))

    x = Symbol('x', real=True)
    assert Abs(f(x)).diff(x).subs(f(x), 1+I*x).doit() == x/sqrt(1 + x**2)
    assert arg(f(x)).diff(x).subs(f(x), 1+I*x**2).doit() == 2*x/(1+x**4)

def test_f_expand_complex():
    f = Function('f')
    x = Symbol('x', real=True)
    z = Symbol('z')

    assert f(x).expand(complex=True)        == I*im(f(x)) + re(f(x))
    assert exp(x).expand(complex=True)      == exp(x)
    assert exp(I*x).expand(complex=True)    == cos(x) + I*sin(x)
    assert exp(z).expand(complex=True)      == cos(im(z))*exp(re(z)) + \
                                             I*sin(im(z))*exp(re(z))

def test_latex_functions():
    assert latex(exp(x)) == "e^{x}"
    assert latex(exp(1)+exp(2)) == "e + e^{2}"

    f = Function('f')
    assert latex(f(x)) == '\\operatorname{f}{\\left (x \\right )}'

    beta = Function('beta')

    assert latex(beta(x)) == r"\beta{\left (x \right )}"
    assert latex(sin(x)) == r"\sin{\left (x \right )}"
    assert latex(sin(x), fold_func_brackets=True) == r"\sin {x}"
    assert latex(sin(2*x**2), fold_func_brackets=True) == \
    r"\sin {2 x^{2}}"
    assert latex(sin(x**2), fold_func_brackets=True) == \
    r"\sin {x^{2}}"

    assert latex(asin(x)**2) == r"\operatorname{asin}^{2}{\left (x \right )}"
    assert latex(asin(x)**2,inv_trig_style="full") == \
        r"\arcsin^{2}{\left (x \right )}"
    assert latex(asin(x)**2,inv_trig_style="power") == \
        r"\sin^{-1}{\left (x \right )}^{2}"
    assert latex(asin(x**2),inv_trig_style="power",fold_func_brackets=True) == \
        r"\sin^{-1} {x^{2}}"

    assert latex(factorial(k)) == r"k!"
    assert latex(factorial(-k)) == r"\left(- k\right)!"

    assert latex(factorial2(k)) == r"k!!"
    assert latex(factorial2(-k)) == r"\left(- k\right)!!"

    assert latex(binomial(2,k)) == r"{\binom{2}{k}}"

    assert latex(FallingFactorial(3,k)) == r"{\left(3\right)}_{\left(k\right)}"
    assert latex(RisingFactorial(3,k)) == r"{\left(3\right)}^{\left(k\right)}"

    assert latex(floor(x)) == r"\lfloor{x}\rfloor"
    assert latex(ceiling(x)) == r"\lceil{x}\rceil"
    assert latex(Abs(x)) == r"\lvert{x}\rvert"
    assert latex(re(x)) == r"\Re{x}"
    assert latex(re(x+y)) == r"\Re {\left (x + y \right )}"
    assert latex(im(x)) == r"\Im{x}"
    assert latex(conjugate(x)) == r"\overline{x}"
    assert latex(gamma(x)) == r"\Gamma\left(x\right)"
    assert latex(Order(x)) == r"\mathcal{O}\left(x\right)"
    assert latex(lowergamma(x, y)) == r'\gamma\left(x, y\right)'
    assert latex(uppergamma(x, y)) == r'\Gamma\left(x, y\right)'

    assert latex(cot(x)) == r'\cot{\left (x \right )}'
    assert latex(coth(x)) == r'\coth{\left (x \right )}'
    assert latex(re(x)) == r'\Re{x}'
    assert latex(im(x)) == r'\Im{x}'
    assert latex(root(x,y)) == r'x^{\frac{1}{y}}'
    assert latex(arg(x)) == r'\arg{\left (x \right )}'
    assert latex(zeta(x)) == r'\zeta{\left (x \right )}'

 def _eval_evalf(self, prec):
     """ Careful! any evalf of polar numbers is flaky """
     from sympy import im, pi, re
     i = im(self.args[0])
     if i <= -pi or i > pi:
         return self  # cannot evalf for this argument
     res = exp(self.args[0])._eval_evalf(prec)
     if i > 0 and im(res) < 0:
         # i ~ pi, but exp(I*i) evaluated to argument slightly bigger than pi
         return re(res)
     return res

def test_map_coeffs():
    p = Poly(x**2 + 2*x*y, x, y)
    q = p.map_coeffs(lambda c: 2*c)

    assert q.as_basic() == 2*x**2 + 4*x*y

    p = Poly(u*x**2 + v*x*y, x, y)
    q = p.map_coeffs(expand, complex=True)

    assert q.as_basic() == x**2*(I*im(u) + re(u)) + x*y*(I*im(v) + re(v))

    raises(PolynomialError, "p.map_coeffs(lambda c: x*c)")

def test_derivatives_issue_4757():
    x = Symbol('x', real=True)
    y = Symbol('y', imaginary=True)
    f = Function('f')
    assert re(f(x)).diff(x) == re(f(x).diff(x))
    assert im(f(x)).diff(x) == im(f(x).diff(x))
    assert re(f(y)).diff(y) == -I*im(f(y).diff(y))
    assert im(f(y)).diff(y) == -I*re(f(y).diff(y))
    assert Abs(f(x)).diff(x).subs(f(x), 1 + I*x).doit() == x/sqrt(1 + x**2)
    assert arg(f(x)).diff(x).subs(f(x), 1 + I*x**2).doit() == 2*x/(1 + x**4)
    assert Abs(f(y)).diff(y).subs(f(y), 1 + y).doit() == -y/sqrt(1 - y**2)
    assert arg(f(y)).diff(y).subs(f(y), I + y**2).doit() == 2*y/(1 + y**4)

def test_as_real_imag():
    n = pi**1000
    # the special code for working out the real
    # and complex parts of a power with Integer exponent
    # should not run if there is no imaginary part, hence
    # this should not hang
    assert n.as_real_imag() == (n, 0)

    # issue 3162
    x = Symbol('x')
    assert sqrt(x).as_real_imag() == \
    ((re(x)**2 + im(x)**2)**(S(1)/4)*cos(atan2(im(x), re(x))/2), \
     (re(x)**2 + im(x)**2)**(S(1)/4)*sin(atan2(im(x), re(x))/2))

    def eval(cls, arg):
        from sympy import im
        if arg.is_integer or arg.is_finite is False:
            return arg
        if arg.is_imaginary or (S.ImaginaryUnit*arg).is_real:
            i = im(arg)
            if not i.has(S.ImaginaryUnit):
                return cls(i)*S.ImaginaryUnit
            return cls(arg, evaluate=False)

        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 im(t).is_integer):
                ipart += t
            elif t.has(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 (S.ImaginaryUnit*spart).is_real) or
                npart.is_imaginary and spart.is_real):
            try:
                r, i = get_integer_part(
                    npart, cls._dir, {}, return_ints=True)
                ipart += Integer(r) + Integer(i)*S.ImaginaryUnit
                npart = S.Zero
            except (PrecisionExhausted, NotImplementedError):
                pass

        spart += npart
        if not spart:
            return ipart
        elif spart.is_imaginary or (S.ImaginaryUnit*spart).is_real:
            return ipart + cls(im(spart), evaluate=False)*S.ImaginaryUnit
        else:
            return ipart + cls(spart, evaluate=False)

def test_zero_assumptions():
    nr = Symbol('nonreal', real=False)
    ni = Symbol('nonimaginary', imaginary=False)
    # imaginary implies not zero
    nzni = Symbol('nonzerononimaginary', zero=False, imaginary=False)

    assert re(nr).is_zero is None
    assert im(nr).is_zero is False

    assert re(ni).is_zero is None
    assert im(ni).is_zero is None

    assert re(nzni).is_zero is False
    assert im(nzni).is_zero is None

def test_diff_no_eval_derivative():
    class My(Expr):
        def __new__(cls, x):
            return Expr.__new__(cls, x)

    x, y = symbols('x y')
    # My doesn't have its own _eval_derivative method
    assert My(x).diff(x).func is Derivative
    assert My(x).diff(x, 3).func is Derivative
    assert re(x).diff(x, 2) == Derivative(re(x), (x, 2))  # issue 15518
    assert diff(NDimArray([re(x), im(x)]), (x, 2)) == NDimArray(
        [Derivative(re(x), (x, 2)), Derivative(im(x), (x, 2))])
    # it doesn't have y so it shouldn't need a method for this case
    assert My(x).diff(y) == 0

def test_real_imag():
    x, y, z = symbols('x, y, z')
    X, Y, Z = symbols('X, Y, Z', commutative=False)
    a = Symbol('a', real=True)
    assert (2*a*x).as_real_imag() == (2*a*re(x), 2*a*im(x))

    # issue 2296:
    assert (x*x.conjugate()).as_real_imag() == (Abs(x)**2, 0)
    assert im(x*x.conjugate()) == 0
    assert im(x*y.conjugate()*z*y) == im(x*z)*Abs(y)**2
    assert im(x*y.conjugate()*x*y) == im(x**2)*Abs(y)**2
    assert im(Z*y.conjugate()*X*y) == im(Z*X)*Abs(y)**2
    assert im(X*X.conjugate()) == im(X*X.conjugate(), evaluate=False)
    assert (sin(x)*sin(x).conjugate()).as_real_imag() == \
        (Abs(sin(x))**2, 0)

def test_re():
    x, y = symbols('x,y')

    r = Symbol('r', real=True)
    i = Symbol('i', imaginary=True)

    assert re(nan) == nan

    assert re(oo) == oo
    assert re(-oo) == -oo

    assert re(0) == 0

    assert re(1) == 1
    assert re(-1) == -1

    assert re(E) == E
    assert re(-E) == -E

    assert re(x) == re(x)
    assert re(x*I) == -im(x)
    assert re(r*I) == 0
    assert re(r) == r
    assert re(i*I) == I * i
    assert re(i) == 0

    assert re(x + y) == re(x + y)
    assert re(x + r) == re(x) + r

    assert re(re(x)) == re(x)

    assert re(2 + I) == 2
    assert re(x + I) == re(x)

    assert re(x + y*I) == re(x) - im(y)
    assert re(x + r*I) == re(x)

    assert re(log(2*I)) == log(2)

    assert re((2+I)**2).expand(complex=True) == 3

    assert re(conjugate(x)) == re(x)
    assert conjugate(re(x)) == re(x)

    assert re(x).as_real_imag() == (re(x), 0)

    assert re(i*r*x).diff(r) == re(i*x)
    assert re(i*r*x).diff(i) == -I * im(r*x)

def test_reduce_poly_inequalities_complex_relational():
    cond = Eq(im(x), 0)

    assert reduce_poly_inequalities(
        [[Eq(x**2, 0)]], x, relational=True) == And(Eq(re(x), 0), cond)
    assert reduce_poly_inequalities(
        [[Le(x**2, 0)]], x, relational=True) == And(Eq(re(x), 0), cond)
    assert reduce_poly_inequalities(
        [[Lt(x**2, 0)]], x, relational=True) is False
    assert reduce_poly_inequalities(
        [[Ge(x**2, 0)]], x, relational=True) == cond
    assert reduce_poly_inequalities([[Gt(x**2, 0)]], x, relational=True) == And(Or(Lt(re(x), 0), Lt(0, re(x))), cond)
    assert reduce_poly_inequalities([[Ne(x**2, 0)]], x, relational=True) == And(Or(Lt(re(x), 0), Lt(0, re(x))), cond)

    assert reduce_poly_inequalities([[Eq(x**2, 1)]], x, relational=True) == And(Or(Eq(re(x), -1), Eq(re(x), 1)), cond)
    assert reduce_poly_inequalities([[Le(x**2, 1)]], x, relational=True) == And(And(Le(-1, re(x)), Le(re(x), 1)), cond)
    assert reduce_poly_inequalities([[Lt(x**2, 1)]], x, relational=True) == And(And(Lt(-1, re(x)), Lt(re(x), 1)), cond)
    assert reduce_poly_inequalities([[Ge(x**2, 1)]], x, relational=True) == And(Or(Le(re(x), -1), Le(1, re(x))), cond)
    assert reduce_poly_inequalities([[Gt(x**2, 1)]], x, relational=True) == And(Or(Lt(re(x), -1), Lt(1, re(x))), cond)
    assert reduce_poly_inequalities([[Ne(x**2, 1)]], x, relational=True) == And(Or(Lt(re(x), -1), And(Lt(-1, re(x)), Lt(re(x), 1)), Lt(1, re(x))), cond)

    assert reduce_poly_inequalities([[Eq(x**2, 1.0)]], x, relational=True).evalf() == And(Or(Eq(re(x), -1.0), Eq(re(x), 1.0)), cond)
    assert reduce_poly_inequalities([[Le(x**2, 1.0)]], x, relational=True) == And(And(Le(-1.0, re(x)), Le(re(x), 1.0)), cond)
    assert reduce_poly_inequalities([[Lt(x**2, 1.0)]], x, relational=True) == And(And(Lt(-1.0, re(x)), Lt(re(x), 1.0)), cond)
    assert reduce_poly_inequalities([[Ge(x**2, 1.0)]], x, relational=True) == And(Or(Le(re(x), -1.0), Le(1.0, re(x))), cond)
    assert reduce_poly_inequalities([[Gt(x**2, 1.0)]], x, relational=True) == And(Or(Lt(re(x), -1.0), Lt(1.0, re(x))), cond)
    assert reduce_poly_inequalities([[Ne(x**2, 1.0)]], x, relational=True) == And(Or(Lt(re(x), -1.0), And(Lt(-1.0, re(x)), Lt(re(x), 1.0)), Lt(1.0, re(x))), cond)

def get_integer_part(expr, no, options, return_ints=False):
    """
    With no = 1, computes ceiling(expr)
    With no = -1, computes floor(expr)

    Note: this function either gives the exact result or signals failure.
    """
    import sympy
    # The expression is likely less than 2^30 or so
    assumed_size = 30
    ire, iim, ire_acc, iim_acc = evalf(expr, assumed_size, options)

    # We now know the size, so we can calculate how much extra precision
    # (if any) is needed to get within the nearest integer
    if ire and iim:
        gap = max(fastlog(ire) - ire_acc, fastlog(iim) - iim_acc)
    elif ire:
        gap = fastlog(ire) - ire_acc
    elif iim:
        gap = fastlog(iim) - iim_acc
    else:
        # ... or maybe the expression was exactly zero
        return None, None, None, None

    margin = 10

    if gap >= -margin:
        ire, iim, ire_acc, iim_acc = \
            evalf(expr, margin + assumed_size + gap, options)

    # We can now easily find the nearest integer, but to find floor/ceil, we
    # must also calculate whether the difference to the nearest integer is
    # positive or negative (which may fail if very close).
    def calc_part(expr, nexpr):
        from sympy import Add
        nint = int(to_int(nexpr, rnd))
        n, c, p, b = nexpr
        if (c != 1 and p != 0) or p < 0:
            expr = Add(expr, -nint, evaluate=False)
            x, _, x_acc, _ = evalf(expr, 10, options)
            try:
                check_target(expr, (x, None, x_acc, None), 3)
            except PrecisionExhausted:
                if not expr.equals(0):
                    raise PrecisionExhausted
                x = fzero
            nint += int(no*(mpf_cmp(x or fzero, fzero) == no))
        nint = from_int(nint)
        return nint, fastlog(nint) + 10

    re, im, re_acc, im_acc = None, None, None, None

    if ire:
        re, re_acc = calc_part(sympy.re(expr, evaluate=False), ire)
    if iim:
        im, im_acc = calc_part(sympy.im(expr, evaluate=False), iim)

    if return_ints:
        return int(to_int(re or fzero)), int(to_int(im or fzero))
    return re, im, re_acc, im_acc