def eval(cls, nu, z):
        if z.is_zero:
            if nu.is_zero:
                return S.One
            elif (nu.is_integer and nu.is_zero is False) or re(nu).is_positive:
                return S.Zero
            elif re(nu).is_negative and not (nu.is_integer is True):
                return S.ComplexInfinity
            elif nu.is_imaginary:
                return S.NaN
        if z is S.Infinity or (z is S.NegativeInfinity):
            return S.Zero

        if z.could_extract_minus_sign():
            return (z)**nu*(-z)**(-nu)*besselj(nu, -z)
        if nu.is_integer:
            if nu.could_extract_minus_sign():
                return S(-1)**(-nu)*besselj(-nu, z)
            newz = z.extract_multiplicatively(I)
            if newz:  # NOTE we don't want to change the function if z==0
                return I**(nu)*besseli(nu, newz)

        # branch handling:
        from sympy import unpolarify, exp
        if nu.is_integer:
            newz = unpolarify(z)
            if newz != z:
                return besselj(nu, newz)
        else:
            newz, n = z.extract_branch_factor()
            if n != 0:
                return exp(2*n*pi*nu*I)*besselj(nu, newz)
        nnu = unpolarify(nu)
        if nu != nnu:
            return besselj(nnu, z)

    def eval(cls, nu, z):
        if z.is_zero:
            if nu.is_zero:
                return S.NegativeInfinity
            elif re(nu).is_zero is False:
                return S.ComplexInfinity
            elif re(nu).is_zero:
                return S.NaN
        if z is S.Infinity or z is S.NegativeInfinity:
            return S.Zero

        if nu.is_integer:
            if nu.could_extract_minus_sign():
                return S(-1) ** (-nu) * bessely(-nu, z)

def test_solve_sqrt_3():
    R = Symbol("R")
    eq = sqrt(2) * R * sqrt(1 / (R + 1)) + (R + 1) * (sqrt(2) * sqrt(1 / (R + 1)) - 1)
    sol = solveset_complex(eq, R)

    assert sol == FiniteSet(
        *[
            S(5) / 3 + 4 * sqrt(10) * cos(atan(3 * sqrt(111) / 251) / 3) / 3,
            -sqrt(10) * cos(atan(3 * sqrt(111) / 251) / 3) / 3
            + 40 * re(1 / ((-S(1) / 2 - sqrt(3) * I / 2) * (S(251) / 27 + sqrt(111) * I / 9) ** (S(1) / 3))) / 9
            + sqrt(30) * sin(atan(3 * sqrt(111) / 251) / 3) / 3
            + S(5) / 3
            + I
            * (
                -sqrt(30) * cos(atan(3 * sqrt(111) / 251) / 3) / 3
                - sqrt(10) * sin(atan(3 * sqrt(111) / 251) / 3) / 3
                + 40 * im(1 / ((-S(1) / 2 - sqrt(3) * I / 2) * (S(251) / 27 + sqrt(111) * I / 9) ** (S(1) / 3))) / 9
            ),
        ]
    )

    # the number of real roots will depend on the value of m: for m=1 there are 4
    # and for m=-1 there are none.
    eq = -sqrt((m - q) ** 2 + (-m / (2 * q) + S(1) / 2) ** 2) + sqrt(
        (-m ** 2 / 2 - sqrt(4 * m ** 4 - 4 * m ** 2 + 8 * m + 1) / 4 - S(1) / 4) ** 2
        + (m ** 2 / 2 - m - sqrt(4 * m ** 4 - 4 * m ** 2 + 8 * m + 1) / 4 - S(1) / 4) ** 2
    )
    raises(NotImplementedError, lambda: solveset_real(eq, q))

    def as_real_imag(self):
        """Returns a tuple of the real part of the input Matrix
           and it's imaginary part.

           >>> from sympy import Matrix, I
           >>> A = Matrix([[1+2*I, 3], [4+7*I, 5]])
           >>> A.as_real_imag()
           (Matrix([
           [1, 3],
           [4, 5]]), Matrix([
           [2, 0],
           [7, 0]]))
           >>> from sympy.abc import x, y, z, w
           >>> B = Matrix([[x, y + x * I], [z + w * I, z]])
           >>> B.as_real_imag()
           (Matrix([
           [        re(x), re(y) - im(x)],
           [re(z) - im(w),         re(z)]]), Matrix([
           [        im(x), re(x) + im(y)],
           [re(w) + im(z),         im(z)]]))

        """
        from sympy.functions.elementary.complexes import re, im
        real_mat = self._new(self.rows, self.cols, lambda i, j: re(self[i, j]))
        im_mat = self._new(self.rows, self.cols, lambda i, j: im(self[i, j]))

        return (real_mat, im_mat)

    def as_real_imag(self):
        real_matrices = [re(matrix) for matrix in self.blocks]
        real_matrices = Matrix(self.blockshape[0], self.blockshape[1], real_matrices)

        im_matrices = [im(matrix) for matrix in self.blocks]
        im_matrices = Matrix(self.blockshape[0], self.blockshape[1], im_matrices)

        return (real_matrices, im_matrices)

 def _eval_rewrite_as_besseli(self, z):
     ot = Rational(1, 3)
     tt = Rational(2, 3)
     a = Pow(z, Rational(3, 2))
     if re(z).is_positive:
         return ot * sqrt(z) * (besseli(-ot, tt * a) - besseli(ot, tt * a))
     else:
         return ot * (Pow(a, ot) * besseli(-ot, tt * a) - z * Pow(a, -ot) * besseli(ot, tt * a))

    def as_real_imag(self):
        """Returns tuple containing (real , imaginary) part of sparse matrix"""
        from sympy.functions.elementary.complexes import re, im
        real_smat = self.copy()
        im_smat = self.copy()

        for key, value in self._smat.items():
            real_smat._smat[key] = re(value)
            im_smat._smat[key] = im(value)

        return (real_smat, im_smat)

 def _eval_evalf(self, prec):
     """ Careful! any evalf of polar numbers is flaky """
     from sympy import im, pi, re
     i = im(self.args[0])
     try:
         bad = (i <= -pi or i > pi)
     except TypeError:
         bad = True
     if bad:
         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_solve_sqrt_3():
    R = Symbol("R")
    eq = sqrt(2) * R * sqrt(1 / (R + 1)) + (R + 1) * (sqrt(2) * sqrt(1 / (R + 1)) - 1)
    sol = solveset_complex(eq, R)

    assert sol == FiniteSet(
        *[
            S(5) / 3 + 4 * sqrt(10) * cos(atan(3 * sqrt(111) / 251) / 3) / 3,
            -sqrt(10) * cos(atan(3 * sqrt(111) / 251) / 3) / 3
            + 40 * re(1 / ((-S(1) / 2 - sqrt(3) * I / 2) * (S(251) / 27 + sqrt(111) * I / 9) ** (S(1) / 3))) / 9
            + sqrt(30) * sin(atan(3 * sqrt(111) / 251) / 3) / 3
            + S(5) / 3
            + I
            * (
                -sqrt(30) * cos(atan(3 * sqrt(111) / 251) / 3) / 3
                - sqrt(10) * sin(atan(3 * sqrt(111) / 251) / 3) / 3
                + 40 * im(1 / ((-S(1) / 2 - sqrt(3) * I / 2) * (S(251) / 27 + sqrt(111) * I / 9) ** (S(1) / 3))) / 9
            ),
        ]
    )

    # the number of real roots will depend on the value of m: for m=1 there are 4
    # and for m=-1 there are none.
    eq = -sqrt((m - q) ** 2 + (-m / (2 * q) + S(1) / 2) ** 2) + sqrt(
        (-m ** 2 / 2 - sqrt(4 * m ** 4 - 4 * m ** 2 + 8 * m + 1) / 4 - S(1) / 4) ** 2
        + (m ** 2 / 2 - m - sqrt(4 * m ** 4 - 4 * m ** 2 + 8 * m + 1) / 4 - S(1) / 4) ** 2
    )
    unsolved_object = ConditionSet(
        q,
        Eq(
            (
                -2 * sqrt(4 * q ** 2 * (m - q) ** 2 + (-m + q) ** 2)
                + sqrt(
                    (-2 * m ** 2 - sqrt(4 * m ** 4 - 4 * m ** 2 + 8 * m + 1) - 1) ** 2
                    + (2 * m ** 2 - 4 * m - sqrt(4 * m ** 4 - 4 * m ** 2 + 8 * m + 1) - 1) ** 2
                )
                * Abs(q)
            )
            / Abs(q),
            0,
        ),
        S.Reals,
    )
    assert solveset_real(eq, q) == unsolved_object

    def eval(cls, n, a, x):
        # For negative n the polynomials vanish
        # See http://functions.wolfram.com/Polynomials/GegenbauerC3/03/01/03/0012/
        if n.is_negative:
            return S.Zero

        # Some special values for fixed a
        if a == S.Half:
            return legendre(n, x)
        elif a == S.One:
            return chebyshevu(n, x)
        elif a == S.NegativeOne:
            return S.Zero

        if not n.is_Number:
            # Handle this before the general sign extraction rule
            if x == S.NegativeOne:
                if (re(a) > S.Half) == True:
                    return S.ComplexInfinity
                else:
                    # No sec function available yet
                    #return (cos(S.Pi*(a+n)) * sec(S.Pi*a) * gamma(2*a+n) /
                    #            (gamma(2*a) * gamma(n+1)))
                    return None

            # Symbolic result C^a_n(x)
            # C^a_n(-x)  --->  (-1)**n * C^a_n(x)
            if x.could_extract_minus_sign():
                return S.NegativeOne**n * gegenbauer(n, a, -x)
            # We can evaluate for some special values of x
            if x == S.Zero:
                return (2**n * sqrt(S.Pi) * gamma(a + S.Half*n) /
                        (gamma((1 - n)/2) * gamma(n + 1) * gamma(a)) )
            if x == S.One:
                return gamma(2*a + n) / (gamma(2*a) * gamma(n + 1))
            elif x == S.Infinity:
                if n.is_positive:
                    return RisingFactorial(a, n) * S.Infinity
        else:
            # n is a given fixed integer, evaluate into polynomial
            return gegenbauer_poly(n, a, x)

 def _check(self, value):
     from sympy.functions import re, im
     super(ComplexType, self)._check(re(value))
     super(ComplexType, self)._check(im(value))

 def _cast_nocheck(self, value):
     from sympy.functions import re, im
     return (
         super(ComplexType, self)._cast_nocheck(re(value)) +
         super(ComplexType, self)._cast_nocheck(im(value))*1j
     )

 def _eval_rewrite_as_besselj(self, z):
     ot = Rational(1, 3)
     tt = Rational(2, 3)
     a = Pow(-z, Rational(3, 2))
     if re(z).is_negative:
         return ot * sqrt(-z) * (besselj(-ot, tt * a) + besselj(ot, tt * a))

 def _eval_Abs(self):   # Abs is never a polar number
     from sympy.functions.elementary.complexes import re
     return exp(re(self.args[0]))

 def as_real_imag(self):
     return (re(Matrix(self)), im(Matrix(self)))

 def _eval_rewrite_as_besselj(self, z, **kwargs):
     ot = Rational(1, 3)
     tt = Rational(2, 3)
     a = Pow(-z, Rational(3, 2))
     if re(z).is_negative:
         return sqrt(-z/3) * (besselj(-ot, tt*a) - besselj(ot, tt*a))

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.
    """
    from sympy.functions.elementary.complexes import re, im
    # 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.core.add import Add
        nint = int(to_int(nexpr, rnd))
        n, c, p, b = nexpr
        is_int = (p == 0)
        if not is_int:
            # if there are subs and they all contain integer re/im parts
            # then we can (hopefully) safely substitute them into the
            # expression
            s = options.get('subs', False)
            if s:
                doit = True
                from sympy.core.compatibility import as_int
                for v in s.values():
                    try:
                        as_int(v)
                    except ValueError:
                        try:
                            [as_int(i) for i in v.as_real_imag()]
                            continue
                        except (ValueError, AttributeError):
                            doit = False
                            break
                if doit:
                    expr = expr.subs(s)

            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(re(expr, evaluate=False), ire)
    if iim:
        im_, im_acc = calc_part(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