def _eval_evalf(self, prec):
        # The default code is insufficient for polar arguments.
        # mpmath provides an optional argument "r", which evaluates
        # G(z**(1/r)). I am not sure what its intended use is, but we hijack it
        # here in the following way: to evaluate at a number z of |argument|
        # less than (say) n*pi, we put r=1/n, compute z' = root(z, n)
        # (carefully so as not to loose the branch information), and evaluate
        # G(z'**(1/r)) = G(z'**n) = G(z).
        from sympy.functions import exp_polar, ceiling
        from sympy import Expr
        import mpmath
        z = self.argument
        znum = self.argument._eval_evalf(prec)
        if znum.has(exp_polar):
            znum, branch = znum.as_coeff_mul(exp_polar)
            if len(branch) != 1:
                return
            branch = branch[0].args[0]/I
        else:
            branch = S(0)
        n = ceiling(abs(branch/S.Pi)) + 1
        znum = znum**(S(1)/n)*exp(I*branch / n)

        # Convert all args to mpf or mpc
        try:
            [z, r, ap, bq] = [arg._to_mpmath(prec)
                    for arg in [znum, 1/n, self.args[0], self.args[1]]]
        except ValueError:
            return

        with mpmath.workprec(prec):
            v = mpmath.meijerg(ap, bq, z, r)

        return Expr._from_mpmath(v, prec)

 def _eval_evalf(self, prec):
     from mpmath import mp, workprec
     from sympy import Expr
     z = self.args[0]._to_mpmath(prec)
     with workprec(prec):
         res = mp.airybi(z, derivative=1)
     return Expr._from_mpmath(res, prec)

 def _eval_evalf(self, prec):
     z = self.argument._to_mpmath(prec)
     ap = [a._to_mpmath(prec) for a in self.ap]
     bp = [b._to_mpmath(prec) for b in self.bq]
     with mpmath.workprec(prec):
         res = mpmath.hyper(ap, bp, z, eliminate=False)
     return Expr._from_mpmath(res, prec)

    def apply(self, z, evaluation):
        "%(name)s[z__]"

        args = z.get_sequence()

        if len(args) != self.nargs:
            return

        # if no arguments are inexact attempt to use sympy
        if len([True for x in args if Expression("InexactNumberQ", x).evaluate(evaluation).is_true()]) == 0:
            expr = Expression(self.get_name(), *args).to_sympy()
            result = from_sympy(expr)
            # evaluate leaves to convert e.g. Plus[2, I] -> Complex[2, 1]
            result = result.evaluate_leaves(evaluation)
        else:
            prec = min_prec(*args)
            with mpmath.workprec(prec):
                mpmath_args = [sympy2mpmath(x.to_sympy()) for x in args]
                if None in mpmath_args:
                    return
                try:
                    result = self.eval(*mpmath_args)
                    result = from_sympy(mpmath2sympy(result, prec))
                except ValueError, exc:
                    text = str(exc)
                    if text == "gamma function pole":
                        return Symbol("ComplexInfinity")
                    else:
                        raise
                except ZeroDivisionError:
                    return
                except SpecialValueError, exc:
                    return Symbol(exc.name)

 def apply_N(self, prec, evaluation):
     'N[E, prec_]'
     
     prec = get_precision(prec, evaluation)
     if prec is not None:
         with workprec(prec):
             return Real(mpmath2gmpy(mpmath.e))

    def apply(self, z, evaluation):
        '%(name)s[z__]'

        args = z.get_sequence()

        if len(args) != self.nargs:
            return

        # if no arguments are inexact attempt to use sympy
        if all(not x.is_inexact() for x in args):
            result = Expression(self.get_name(), *args).to_sympy()
            result = self.prepare_mathics(result)
            result = from_sympy(result)
            # evaluate leaves to convert e.g. Plus[2, I] -> Complex[2, 1]
            result = result.evaluate_leaves(evaluation)
        else:
            prec = min_prec(*args)
            with mpmath.workprec(prec):
                mpmath_args = [sympy2mpmath(x.to_sympy()) for x in args]
                if None in mpmath_args:
                    return
                try:
                    result = self.eval(*mpmath_args)
                    result = from_sympy(mpmath2sympy(result, prec))
                except ValueError, exc:
                    text = str(exc)
                    if text == 'gamma function pole':
                        return Symbol('ComplexInfinity')
                    else:
                        raise
                except ZeroDivisionError:
                    return
                except SpecialValueError, exc:
                    return Symbol(exc.name)

    def _eval_evalf(self, prec):
        m = self.args[0]

        if m.is_Integer and m.is_nonnegative:
            m = m._to_mpmath(prec)
            with workprec(prec):
                res = mp.eulernum(m)
            return Expr._from_mpmath(res, prec)

 def _eval_evalf(self, prec):
     from mpmath import mp, workprec
     from sympy import Expr
     a = self.args[0]._to_mpmath(prec)
     z = self.args[1]._to_mpmath(prec)
     with workprec(prec):
         res = mp.gammainc(a, z, mp.inf)
     return Expr._from_mpmath(res, prec)

def real_power(x, y):
    x = g_mpf(x)
    y = g_mpf(y)
    prec = min(x.getprec(), y.getprec())
    with workprec(prec):
        x = gmpy2mpmath(x)
        y = gmpy2mpmath(y)
        return mpmath2gmpy(x ** y)

    def _eval_evalf(self, prec):
        """Evaluate this complex root to the given precision."""
        with workprec(prec):
            g = self.poly.gen
            if not g.is_Symbol:
                d = Dummy('x')
                func = lambdify(d, self.expr.subs({g: d}), "mpmath")
            else:
                func = lambdify(g, self.expr, "mpmath")

            try:
                interval = self.interval
            except DomainError:
                return super()._eval_evalf(prec)

            while True:
                if self.is_extended_real:
                    a = mpf(str(interval.a))
                    b = mpf(str(interval.b))
                    if a == b:
                        root = a
                        break
                    x0 = mpf(str(interval.center))
                else:
                    ax = mpf(str(interval.ax))
                    bx = mpf(str(interval.bx))
                    ay = mpf(str(interval.ay))
                    by = mpf(str(interval.by))
                    x0 = mpc(*map(str, interval.center))
                    if ax == bx and ay == by:
                        root = x0
                        break

                try:
                    root = findroot(func, x0)
                    # If the (real or complex) root is not in the 'interval',
                    # then keep refining the interval. This happens if findroot
                    # accidentally finds a different root outside of this
                    # interval because our initial estimate 'x0' was not close
                    # enough. It is also possible that the secant method will
                    # get trapped by a max/min in the interval; the root
                    # verification by findroot will raise a ValueError in this
                    # case and the interval will then be tightened -- and
                    # eventually the root will be found.
                    if self.is_extended_real:
                        if (a <= root <= b):
                            break
                    elif (ax <= root.real <= bx and ay <= root.imag <= by
                          and (interval.ay > 0 or interval.by < 0)):
                        break
                except (ValueError, UnboundLocalError):
                    pass
                self.refine()
                interval = self.interval

        return ((Float._new(root.real._mpf_, prec) if not self.is_imaginary else 0) +
                I*Float._new(root.imag._mpf_, prec))

    def _eval_evalf(self, prec):
        """Evaluate this complex root to the given precision. """
        with workprec(prec):
            g = self.poly.gen
            if not g.is_Symbol:
                d = Dummy('x')
                func = lambdify(d, self.expr.subs(g, d))
            else:
                func = lambdify(g, self.expr)

            interval = self._get_interval()
            if not self.is_real:
                # For complex intervals, we need to keep refining until the
                # imaginary interval is disjunct with other roots, that is,
                # until both ends get refined.
                ay = interval.ay
                by = interval.by
                while interval.ay == ay or interval.by == by:
                    interval = interval.refine()

            while True:
                if self.is_real:
                    x0 = mpf(str(interval.center))
                else:
                    x0 = mpc(*map(str, interval.center))
                try:
                    root = findroot(func, x0, verify=False)
                    # If the (real or complex) root is not in the 'interval',
                    # then keep refining the interval. This happens if findroot
                    # accidentally finds a different root outside of this
                    # interval because our initial estimate 'x0' was not close
                    # enough.
                    if self.is_real:
                        a = mpf(str(interval.a))
                        b = mpf(str(interval.b))
                        if a == b:
                            root = a
                            break
                        if not (a < root < b):
                            raise ValueError("Root not in the interval.")
                    else:
                        ax = mpf(str(interval.ax))
                        bx = mpf(str(interval.bx))
                        ay = mpf(str(interval.ay))
                        by = mpf(str(interval.by))
                        if ax == bx and ay == by:
                            root = ax + S.ImaginaryUnit*by
                            break
                        if not (ax < root.real < bx and ay < root.imag < by):
                            raise ValueError("Root not in the interval.")
                except ValueError:
                    interval = interval.refine()
                    continue
                else:
                    break

        return Float._new(root.real._mpf_, prec) + I*Float._new(root.imag._mpf_, prec)

    def apply_N(self, k, precision, evaluation):
        'N[AiryBiZero[k_Integer], precision_]'

        prec = get_precision(precision, evaluation)
        k_int = k.get_int_value()

        with mpmath.workprec(prec):
            result = mpmath2sympy(mpmath.airybizero(k_int), prec)
        return from_sympy(result)

 def _eval_evalf(self, prec):
     from mpmath import mp, workprec
     from sympy import Expr
     if all(x.is_number for x in self.args):
         a = self.args[0]._to_mpmath(prec)
         z = self.args[1]._to_mpmath(prec)
         with workprec(prec):
             res = mp.gammainc(a, 0, z)
         return Expr._from_mpmath(res, prec)
     else:
         return self

def gmpy2mpmath(value):
    if isinstance(value, mpcomplex):
        return value.to_mpmath()
    else:
        if get_type(value) != 'f':
            value = g_mpf(value)
        with workprec(value.getprec()):
            value = str(g_mpf(value))
            if value and value[0] == '-':
                return -mp_mpf(value[1:])
            else:
                return mp_mpf(value)

 def apply_inexact(self, n, k, evaluation):
     'Binomial[n_?InexactNumberQ, k_?NumberQ]'
     
     with workprec(min_prec(n, k)):
         n = gmpy2mpmath(n.value)
         k = gmpy2mpmath(k.value)
         result = mpmath.binomial(n, k)
         try:
             result = mpmath2gmpy(result)
         except SpecialValueError, exc:
             return Symbol(exc.name)
         number = Number.from_mp(result)
         return number

 def _eval_evalf(self, prec):
     # Note: works without this function by just calling
     #       mpmath for Legendre polynomials. But using
     #       the dedicated function directly is cleaner.
     from mpmath import mp, workprec
     from sympy import Expr
     n = self.args[0]._to_mpmath(prec)
     m = self.args[1]._to_mpmath(prec)
     theta = self.args[2]._to_mpmath(prec)
     phi = self.args[3]._to_mpmath(prec)
     with workprec(prec):
         res = mp.spherharm(n, m, theta, phi)
     return Expr._from_mpmath(res, prec)

    def apply(self, z, evaluation):
        '%(name)s[z__]'

        args = z.numerify(evaluation).get_sequence()
        mpmath_function = self.get_mpmath_function(args)
        result = None

        # if no arguments are inexact attempt to use sympy
        if all(not x.is_inexact() for x in args):
            result = Expression(self.get_name(), *args).to_sympy()
            result = self.prepare_mathics(result)
            result = from_sympy(result)
            # evaluate leaves to convert e.g. Plus[2, I] -> Complex[2, 1]
            return result.evaluate_leaves(evaluation)
        elif mpmath_function is None:
            return

        if not all(isinstance(arg, Number) for arg in args):
            return

        if any(arg.is_machine_precision() for arg in args):
            # if any argument has machine precision then the entire calculation
            # is done with machine precision.
            float_args = [arg.round().get_float_value(permit_complex=True) for arg in args]
            if None in float_args:
                return

            result = self.call_mpmath(mpmath_function, float_args)
            if isinstance(result, (mpmath.mpc, mpmath.mpf)):
                if mpmath.isinf(result) and isinstance(result, mpmath.mpc):
                    result = Symbol('ComplexInfinity')
                elif mpmath.isinf(result) and result > 0:
                    result = Expression('DirectedInfinity', Integer(1))
                elif mpmath.isinf(result) and result < 0:
                    result = Expression('DirectedInfinity', Integer(-1))
                elif mpmath.isnan(result):
                    result = Symbol('Indeterminate')
                else:
                    result = Number.from_mpmath(result)
        else:
            prec = min_prec(*args)
            d = dps(prec)
            args = [Expression('N', arg, Integer(d)).evaluate(evaluation) for arg in args]
            with mpmath.workprec(prec):
                mpmath_args = [x.to_mpmath() for x in args]
                if None in mpmath_args:
                    return
                result = self.call_mpmath(mpmath_function, mpmath_args)
                if isinstance(result, (mpmath.mpc, mpmath.mpf)):
                    result = Number.from_mpmath(result, d)
        return result

 def apply_inexact(self, z, evaluation):
     '%(name)s[z_Real|z_Complex?InexactNumberQ]'
     
     with workprec(z.get_precision()):
         z = gmpy2mpmath(z.value)
         try:
             result = self.eval(z)
         except ValueError, exc:
             text = str(exc)
             if text == 'gamma function pole':
                 return Symbol('ComplexInfinity')
             else:
                 raise
         except ZeroDivisionError:
             return

    def apply_N(self, k, precision, evaluation):
        'N[AiryBiZero[k_Integer], precision_]'

        try:
            d = get_precision(precision, evaluation)
        except PrecisionValueError:
            return

        if d is None:
            p = machine_precision
        else:
            p = _prec(d)

        k_int = k.get_int_value()

        with mpmath.workprec(p):
            result = mpmath.airybizero(k_int)
            return Number.from_mpmath(result, d)

 def apply_inexact(self, z, evaluation):
     '%(name)s[z_Real|z_Complex?InexactNumberQ]'
     
     prec = z.get_precision()
     with mpmath.workprec(prec):
         z = sympy2mpmath(z.to_sympy())
         if z is None:
             return
         try:
             result = self.eval(z)
             result = mpmath2sympy(result, prec)
         except ValueError, exc:
             text = str(exc)
             if text == 'gamma function pole':
                 return Symbol('ComplexInfinity')
             else:
                 raise
         except ZeroDivisionError:
             return