def bessel_Y(nu,z,algorithm="maxima", prec=53):
    r"""
    Implements the "Y-Bessel function", or "Bessel function of the 2nd
    kind", with index (or "order") nu and argument z.
    
    .. note::

       Currently only prec=53 is supported.
    
    Defn::
    
                    cos(pi n)*bessel_J(nu, z) - bessel_J(-nu, z)
                   -------------------------------------------------
                                     sin(nu*pi)
    
    if nu is not an integer and by taking a limit otherwise.
    
    Sometimes bessel_Y(n,z) is denoted Y_n(z) in the literature.
    
    This is computed using Maxima by default.
    
    EXAMPLES::
    
        sage: bessel_Y(2,1.1,"scipy")
        -1.4314714939...
        sage: bessel_Y(2,1.1)   
        -1.4314714939590...
        sage: bessel_Y(3.001,2.1) 
        -1.0299574976424...

    TESTS::

        sage: bessel_Y(2,1.1, algorithm="pari")
        Traceback (most recent call last):
        ...
        NotImplementedError: The Y-Bessel function is only implemented for the maxima and scipy algorithms
    """
    if algorithm=="scipy":
        if prec != 53:
            raise ValueError, "for the scipy algorithm the precision must be 53"
        import scipy.special
        ans = str(scipy.special.yv(float(nu),complex(real(z),imag(z))))
        ans = ans.replace("(","")
        ans = ans.replace(")","")
        ans = ans.replace("j","*I")
        ans = sage_eval(ans)
        return real(ans) if z in RR else ans
    elif algorithm == "maxima":
        if prec != 53:
            raise ValueError, "for the maxima algorithm the precision must be 53"
        return RR(maxima.eval("bessel_y(%s,%s)"%(float(nu),float(z))))
    elif algorithm == "pari":
        raise NotImplementedError, "The Y-Bessel function is only implemented for the maxima and scipy algorithms"
    else:
        raise ValueError, "unknown algorithm '%s'"%algorithm

    def _evalf_(self, x, **kwargs):
        """
        EXAMPLES::

            sage: from sage.functions.airy import airy_ai_simple
            sage: airy_ai_simple(0.0)
            0.355028053887817
            sage: airy_ai_simple(1.0 * I)
            0.331493305432141 - 0.317449858968444*I

        We can use several methods for numerical evaluation::

            sage: airy_ai_simple(3).n(algorithm='mpmath')
            0.00659113935746072
            sage: airy_ai_simple(3).n(algorithm='mpmath', prec=100)
            0.0065911393574607191442574484080
            sage: airy_ai_simple(3).n(algorithm='scipy')  # rel tol 1e-10
            0.006591139357460719
            sage: airy_ai_simple(I).n(algorithm='scipy')  # rel tol 1e-10
            0.33149330543214117 - 0.3174498589684438*I
            
        TESTS::

            sage: parent(airy_ai_simple(3).n(algorithm='scipy'))                                          
            Real Field with 53 bits of precision
            sage: airy_ai_simple(3).n(algorithm='scipy', prec=200)
            Traceback (most recent call last):
            ...
            NotImplementedError: airy_ai not implemented for precision > 53
        """
        algorithm = kwargs.get('algorithm', 'mpmath') or 'mpmath'
        parent = kwargs.get('parent')
        if algorithm == 'scipy':
            if hasattr(parent, 'prec') and parent.prec() > 53:
                raise NotImplementedError("%s not implemented for precision > 53"%self.name())
            from sage.rings.all import RR, CC
            from sage.functions.other import real,imag
            from scipy.special import airy as airy
            if x in RR:
                y = airy(real(x))[0]
                if parent is None:
                    return RR(y)
            else:
                y = airy(complex(real(x),imag(x)))[0]
                if parent is None:
                    return CC(y)
            return parent(y)
        elif algorithm == 'mpmath':
            import mpmath
            from sage.libs.mpmath import utils as mpmath_utils
            return mpmath_utils.call(mpmath.airyai, x, parent=parent)
        else:
            raise ValueError("unknown algorithm '%s'" % algorithm)

    def sexp_0(self,t):
        #convergence radius 1
        t = self._in_prec(t)
        sexp = self.sexp_0
        b = self.b

        #development point -1 convergence radius 1
        if real(t)>1:
            return b**(sexp(t-1))
        if real(t)<0:
            #sage bug, log(z,b) does not work for complex z
            return log(sexp(t+1))/log(b)
	return self.sexp_0_raw(t)

    def show(self, boundary=True, **options):
        r"""
        Plot ``self``.

        EXAMPLES:

        First some lines::

            sage: PD = HyperbolicPlane().PD()
            sage: PD.get_geodesic(0, 1).show()
            Graphics object consisting of 2 graphics primitives
            sage: PD.get_geodesic(0, 0.3+0.8*I).show()
            Graphics object consisting of 2 graphics primitives

        Then some generic geodesics::

            sage: PD.get_geodesic(-0.5, 0.3+0.4*I).show()
            Graphics object consisting of 2 graphics primitives
            sage: PD.get_geodesic(-1, exp(3*I*pi/7)).show(linestyle="dashed", color="red")
            Graphics object consisting of 2 graphics primitives
            sage: PD.get_geodesic(exp(2*I*pi/11), exp(1*I*pi/11)).show(thickness=6, color="orange")
            Graphics object consisting of 2 graphics primitives
        """
        opts = {'axes': False, 'aspect_ratio': 1}
        opts.update(self.graphics_options())
        opts.update(options)
        end_1, end_2 = [CC(k.coordinates()) for k in self.endpoints()]
        bd_1, bd_2 = [CC(k.coordinates()) for k in self.ideal_endpoints()]
        # Check to see if it's a line
        if abs(bd_1 + bd_2) < EPSILON:
            pic = line([end_1, end_2], **opts)
        else:
            # If we are here, we know it's not a line
            # So we compute the center and radius of the circle
            invdet = RR.one() / (real(bd_1)*imag(bd_2) - real(bd_2)*imag(bd_1))
            centerx = (imag(bd_2) - imag(bd_1)) * invdet
            centery = (real(bd_1) - real(bd_2)) * invdet
            center = centerx + I * centery
            radius = RR(abs(bd_1 - center))
            # Now we calculate the angles for the arc
            theta1 = CC(end_1 - center).arg()
            theta2 = CC(end_2 - center).arg()
            theta1, theta2 = sorted([theta1, theta2])
            # Make sure the sector is inside the disk
            if theta2 - theta1 > pi:
                theta1 += 2 * pi
            pic = arc((centerx, centery), radius,
                      sector=(theta1, theta2), **opts)
        if boundary:
            pic += self._model.get_background_graphic()
        return pic

    def isometry_from_fixed_points(self, repel, attract):
        r"""
        Given two fixed points ``repel`` and ``attract`` as complex
        numbers return a hyperbolic isometry with ``repel`` as repelling
        fixed point and ``attract`` as attracting fixed point.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: UHP.isometry_from_fixed_points(2 + I, 3 + I)
            Traceback (most recent call last):
            ...
            ValueError: fixed points of hyperbolic elements must be ideal

            sage: UHP.isometry_from_fixed_points(2, 0)
            Isometry in UHP
            [  -1    0]
            [-1/3 -1/3]

        TESTS::

            sage: UHP = HyperbolicPlane().UHP()
            sage: UHP.isometry_from_fixed_points(0, 4)
            Isometry in UHP
            [  -1    0]
            [-1/5 -1/5]
            sage: UHP.isometry_from_fixed_points(UHP.get_point(0), UHP.get_point(4))
            Isometry in UHP
            [  -1    0]
            [-1/5 -1/5]
        """
        if isinstance(repel, HyperbolicPoint):
            repel = repel._coordinates
        if isinstance(attract, HyperbolicPoint):
            attract = attract._coordinates

        if imag(repel) + imag(attract) > EPSILON:
            raise ValueError("fixed points of hyperbolic elements must be ideal")
        repel = real(repel)
        attract = real(attract)
        if repel == infinity:
            A = self._moebius_sending([infinity, attract, attract + 1],
                                     [infinity, attract, attract + 2])
        elif attract == infinity:
            A = self._moebius_sending([repel, infinity, repel + 1],
                                     [repel, infinity, repel + 2])
        else:
            A = self._moebius_sending([repel, attract, infinity],
                                     [repel, attract, max(repel, attract) + 1])
        return self.get_isometry(A)

    def sexp_1(self,t):
        t = self._in_prec(t)
        sexp = self.sexp_1
        b = self.b
        IM = self.IM
        N = self.N

        #development point 0 convergence radius 2
        if real(t)>1:
            return b**(sexp(t-1))
        if real(t)<0:
            #sage bug, log(z,b) does not work for complex z
            return log(sexp(t+1))/log(b)
	return self.sexp_1_raw(t)

    def show(self, boundary=True, **options):
        r"""
        Plot ``self``.

        EXAMPLES::

            sage: HyperbolicPlane().UHP().get_geodesic(0, 1).show()
            Graphics object consisting of 2 graphics primitives
        """
        opts = {'axes': False, 'aspect_ratio': 1}
        opts.update(self.graphics_options())
        opts.update(options)
        end_1, end_2 = [CC(k.coordinates()) for k in self.endpoints()]
        bd_1, bd_2 = [CC(k.coordinates()) for k in self.ideal_endpoints()]
        if (abs(real(end_1) - real(end_2)) < EPSILON) \
                or CC(infinity) in [end_1, end_2]: #on same vertical line
            # If one of the endpoints is infinity, we replace it with a
            # large finite  point
            if end_1 == CC(infinity):
                end_1 = (real(end_2), (imag(end_2) + 10))
                end_2 = (real(end_2), imag(end_2))
            elif end_2 == CC(infinity):
                end_2 = (real(end_1), (imag(end_1) + 10))
                end_1 = (real(end_1), imag(end_1))
            from sage.plot.line import line
            pic = line((end_1, end_2), **opts)
            if boundary:
                cent = min(bd_1, bd_2)
                bd_dict = {'bd_min': cent - 3, 'bd_max': cent + 3}
                bd_pic = self._model.get_background_graphic(**bd_dict)
                pic = bd_pic + pic
                return pic
        else:
            center = (bd_1 + bd_2)/2 # Circle center
            radius = abs(bd_1 - bd_2)/2
            theta1 = CC(end_1 - center).arg()
            theta2 = CC(end_2 - center).arg()
            if abs(theta1 - theta2) < EPSILON:
                theta2 += pi
            [theta1, theta2] = sorted([theta1, theta2])
            from sage.calculus.var import var
            from sage.plot.plot import parametric_plot
            x = var('x')
            pic = parametric_plot((radius*cos(x) + real(center),
                                   radius*sin(x) + imag(center)),
                                  (x, theta1, theta2), **opts)
            if boundary:
                # We want to draw a segment of the real line.  The
                # computations below compute the projection of the
                # geodesic to the real line, and then draw a little
                # to the left and right of the projection.
                shadow_1, shadow_2 = [real(k) for k in [end_1, end_2]]
                midpoint = (shadow_1 + shadow_2)/2
                length = abs(shadow_1 - shadow_2)
                bd_dict = {'bd_min': midpoint - length, 'bd_max': midpoint +
                           length}
                bd_pic = self._model.get_background_graphic(**bd_dict)
                pic = bd_pic + pic
            return pic

    def image_coordinates(self, x):
        """
        Return the image of the coordinates of the hyperbolic point ``x``
        under ``self``.

        EXAMPLES::

            sage: PD = HyperbolicPlane().PD()
            sage: KM = HyperbolicPlane().KM()
            sage: phi = KM.coerce_map_from(PD)
            sage: phi.image_coordinates(0.5+0.5*I)
            (0.666666666666667, 0.666666666666667)
        """
        return (2*real(x)/(Integer(1) + real(x)**2 + imag(x)**2),
                2*imag(x)/(Integer(1) + real(x)**2 + imag(x)**2))

    def reflection_involution(self):
        r"""
        Return the isometry of the involution fixing the geodesic ``self``.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: g1 = UHP.get_geodesic(0, 1)
            sage: g1.reflection_involution()
            Isometry in UHP
            [ 1  0]
            [ 2 -1]
            sage: UHP.get_geodesic(I, 2*I).reflection_involution()
            Isometry in UHP
            [ 1  0]
            [ 0 -1]
        """
        x, y = [real(k.coordinates()) for k in self.ideal_endpoints()]
        if x == infinity:
            M = matrix([[1, -2*y], [0, -1]])
        elif y == infinity:
            M = matrix([[1, -2*x], [0, -1]])
        else:
            M = matrix([[(x+y)/(y-x), -2*x*y/(y-x)], [2/(y-x), -(x+y)/(y-x)]])
        return self._model.get_isometry(M)

    def super(self,x):
        """
        Development point is x0-1
        """
        if isinstance(x,float) and self.iprec != None:
           x = RealField(self.iprec)(x)
            
        super = self.super
        super_raw = self.super_raw
        b = self.b
        c = self.c

        xt = x - c
        if real(xt)<-0.5:
            return log(super(x+1))/log(b)
        if real(xt)>0.5:
            return b**(super(x-1))
        return super_raw(x)

    def _dist_points(self, p1, p2):
        r"""
        Compute the distance between two points in the Upper Half Plane
        using the hyperbolic metric.

        INPUT:

        - ``p1``, ``p2`` -- the coordinates of the points

        EXAMPLES::

           sage: HyperbolicPlane().UHP()._dist_points(4.0*I, I)
           1.38629436111989
        """
        num = (real(p2) - real(p1))**2 + (imag(p2) - imag(p1))**2
        denom = 2 * imag(p1) * imag(p2)
        if denom == 0:
            return infinity
        return arccosh(1 + num/denom)

    def midpoint(self):  # UHP
        r"""
        Return the (hyperbolic) midpoint of ``self`` if it exists.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: g = UHP.random_geodesic()
            sage: m = g.midpoint()
            sage: d1 = UHP.dist(m, g.start())
            sage: d2 = UHP.dist(m, g.end())
            sage: bool(abs(d1 - d2) < 10**-9)
            True

        Infinite geodesics have no midpoint::

            sage: UHP.get_geodesic(0, 2).midpoint()
            Traceback (most recent call last):
            ...
            ValueError: the length must be finite
        """
        if self.length() == infinity:
            raise ValueError("the length must be finite")

        start = self._start.coordinates()
        end = self._end.coordinates()
        d = self._model._dist_points(start, end) / 2
        S = self.complete()._to_std_geod(start)
        T = matrix([[exp(d), 0], [0, 1]])
        M = S.inverse() * T * S
        if ((real(start - end) < EPSILON)
                or (abs(real(start - end)) < EPSILON
                    and imag(start - end) < EPSILON)):
            end_p = start
        else:
            end_p = end
        return self._model.get_point(mobius_transform(M, end_p))

    def symmetry_involution(self):
        r"""
        Return the involutory isometry fixing the given point.

        EXAMPLES::

            sage: HyperbolicPlane().UHP().get_point(3 + 2*I).symmetry_involution()
            Isometry in UHP
            [  3/2 -13/2]
            [  1/2  -3/2]
        """
        p = self._coordinates
        x, y = real(p), imag(p)
        if y > 0:
            M = matrix([[x/y, -(x**2/y) - y], [1/y, -(x/y)]])
            return self.parent().get_isometry(M)
        raise ValueError("cannot determine the isometry of a boundary point")

def SL2R_to_SO21(A):
    r"""
    Given a matrix in `SL(2, \RR)` return its irreducible representation in
    `O(2,1)`.

    Note that this is not the only homomorphism, but it is the only one
    that works in the context of the implemented 2D hyperbolic geometry
    models.

    EXAMPLES::

        sage: from sage.geometry.hyperbolic_space.hyperbolic_coercion import SL2R_to_SO21
        sage: A = SL2R_to_SO21(identity_matrix(2))
        sage: J = matrix([[1,0,0],[0,1,0],[0,0,-1]]) #Lorentzian Gram matrix
        sage: norm(A.transpose()*J*A - J) < 10**-4
        True
    """
    a, b, c, d = (A/A.det().sqrt()).list()

    # Kill ~0 imaginary parts
    components = [
        a*d + b*c, a*c - b*d, a*c + b*d, a*b - c*d,
        Integer(1)/Integer(2)*a**2 - Integer(1)/Integer(2)*b**2 -
                Integer(1)/Integer(2)*c**2 + Integer(1)/Integer(2)*d**2,
        Integer(1)/Integer(2)*a**2 + Integer(1)/Integer(2)*b**2 -
                Integer(1)/Integer(2)*c**2 - Integer(1)/Integer(2)*d**2,
        a*b + c*d, Integer(1)/Integer(2)*a**2 -
                Integer(1)/Integer(2)*b**2 + Integer(1)/Integer(2)*c**2 -
        Integer(1)/Integer(2)*d**2, Integer(1)/Integer(2)*a**2 +
                Integer(1)/Integer(2)*b**2 + Integer(1)/Integer(2)*c**2 +
        Integer(1)/Integer(2)*d**2
    ]
    B = matrix(3, [real(c) for c in components])

    #B = B.apply_map(attrcall('real'))
    if A.det() > 0:
        return B
    else:
        # Orientation-reversing isometries swap the nappes of
        #  the lightcone.  This fixes that issue.
        return -B

    def show(self, boundary=True, **options):
        r"""
        Plot ``self``.

        EXAMPLES::

            sage: HyperbolicPlane().UHP().get_point(I).show()
            Graphics object consisting of 2 graphics primitives
            sage: HyperbolicPlane().UHP().get_point(0).show()
            Graphics object consisting of 2 graphics primitives
            sage: HyperbolicPlane().UHP().get_point(infinity).show()
            Traceback (most recent call last):
            ...
            NotImplementedError: can't draw the point infinity
        """
        p = self.coordinates()
        if p == infinity:
            raise NotImplementedError("can't draw the point infinity")
        opts = {'axes': False, 'aspect_ratio': 1}
        opts.update(self.graphics_options())
        opts.update(options)
        from sage.misc.functional import numerical_approx
        p = numerical_approx(p + 0 * I)
        from sage.plot.point import point
        if self._bdry:
            pic = point((p, 0), **opts)
            if boundary:
                bd_pic = self.parent().get_background_graphic(bd_min=p - 1,
                                                              bd_max=p + 1)
                pic = bd_pic + pic
        else:
            pic = point(p, **opts)
            if boundary:
                cent = real(p)
                bd_pic = self.parent().get_background_graphic(bd_min=cent - 1,
                                                              bd_max=cent + 1)
                pic = bd_pic + pic
        return pic

    def ideal_endpoints(self):
        r"""
        Determine the ideal (boundary) endpoints of the complete
        hyperbolic geodesic corresponding to ``self``.

        OUTPUT:

        - a list of 2 boundary points

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: UHP.get_geodesic(I, 2*I).ideal_endpoints()
            [Boundary point in UHP 0,
             Boundary point in UHP +Infinity]
            sage: UHP.get_geodesic(1 + I, 2 + 4*I).ideal_endpoints()
            [Boundary point in UHP -sqrt(65) + 9,
             Boundary point in UHP sqrt(65) + 9]
        """
        start = self._start.coordinates()
        end = self._end.coordinates()
        [x1, x2] = [real(k) for k in [start, end]]
        [y1, y2] = [imag(k) for k in [start, end]]
        M = self._model
        # infinity is the first endpoint, so the other ideal endpoint
        # is just the real part of the second coordinate
        if start == infinity:
            return [M.get_point(start), M.get_point(x2)]
        # Same idea as above
        if end == infinity:
            return [M.get_point(x1), M.get_point(end)]
        # We could also have a vertical line with two interior points
        if x1 == x2:
            return [M.get_point(x1), M.get_point(infinity)]
        # Otherwise, we have a semicircular arc in the UHP
        c = ((x1+x2)*(x2-x1) + (y1+y2)*(y2-y1)) / (2*(x2-x1))
        r = sqrt((c - x1)**2 + y1**2)
        return [M.get_point(c - r), M.get_point(c + r)]

 def cmp_ir(self,z):
     """ 
     returns -1 for left, 0 for in, and 1 for right from initial region
     cut line is on the north ray from L.
     """
     L = self.L
     x0 = self.x0
     if x0 > 0.5:
         if real(z) > real(L) and abs(z) < abs(L):
             return 0
         if real(z) < real(L):
             return -1
         if real(z) > real(L):
             return 1
     else:
         if imag(z) > imag(L):
             if real(z) > real(L):
                 return 1
             if real(z) < real(L):
                 return -1
         if real(z) < real(L) and real(z) > log(real(L)) + log(sqrt(1+tan(imag(z))**2)):
             return 0
         if real(z) > real(L):
             return 1
         if real(z) < real(L):
             return -1

def bessel_K(nu,z,algorithm="pari",prec=53):
    r"""
    Implements the "K-Bessel function", or "modified Bessel function,
    2nd kind", with index (or "order") nu and argument z. Defn::
    
                    pi*(bessel_I(-nu, z) - bessel_I(nu, z))
                   ----------------------------------------
                                2*sin(pi*nu)
            
    
    if nu is not an integer and by taking a limit otherwise.
    
    Sometimes bessel_K(nu,z) is denoted K_nu(z) in the literature. In
    PARI, nu can be complex and z must be real and positive.
    
    EXAMPLES::
    
        sage: bessel_K(3,2,"scipy")
        0.64738539094...
        sage: bessel_K(3,2)
        0.64738539094...
        sage: bessel_K(1,1)
        0.60190723019...
        sage: bessel_K(1,1,"pari",10)
        0.60
        sage: bessel_K(1,1,"pari",100)
        0.60190723019723457473754000154

    TESTS::

        sage: bessel_K(2,1.1, algorithm="maxima")
        Traceback (most recent call last):
        ...
        NotImplementedError: The K-Bessel function is only implemented for the pari and scipy algorithms

        Check whether the return value is real whenever the argument is real (#10251)::

        sage: bessel_K(5, 1.5, algorithm='scipy') in RR
        True

    """
    if algorithm=="scipy":
        if prec != 53:
            raise ValueError, "for the scipy algorithm the precision must be 53"
        import scipy.special
        ans = str(scipy.special.kv(float(nu),float(z)))
        ans = ans.replace("(","")
        ans = ans.replace(")","")
        ans = ans.replace("j","*I")
        ans = sage_eval(ans)
        return real(ans) if z in RR else ans
    elif algorithm == 'pari':
        from sage.libs.pari.all import pari
        try:
            R = RealField(prec)
            nu = R(nu)
            z = R(z)
        except TypeError:
            C = ComplexField(prec)
            nu = C(nu)
            z = C(z)
            K = C
        K = z.parent()
        return K(pari(nu).besselk(z, precision=prec))
    elif algorithm == 'maxima':
        raise NotImplementedError, "The K-Bessel function is only implemented for the pari and scipy algorithms"
    else:
        raise ValueError, "unknown algorithm '%s'"%algorithm

def bessel_J(nu,z,algorithm="pari",prec=53):
    r"""
    Return value of the "J-Bessel function", or "Bessel function, 1st
    kind", with index (or "order") nu and argument z.
    
    ::
    
            Defn:
            Maxima:
                             inf
                            ====          - nu - 2 k  nu + 2 k
                            \     (-1)^k 2           z
                             >    -------------------------
                            /        k! Gamma(nu + k + 1)
                            ====
                            k = 0
        
            PARI:
            
                             inf
                            ====          - 2k    2k
                            \     (-1)^k 2      z    Gamma(nu + 1)
                             >    ----------------------------
                            /         k! Gamma(nu + k + 1)
                            ====
                            k = 0
            
    
    Sometimes bessel_J(nu,z) is denoted J_nu(z) in the literature.
    
    .. warning::

       Inaccurate for small values of z.
    
    EXAMPLES::
    
        sage: bessel_J(2,1.1)
        0.136564153956658
        sage: bessel_J(0,1.1)
        0.719622018527511
        sage: bessel_J(0,1) 
        0.765197686557967
        sage: bessel_J(0,0)
        1.00000000000000
        sage: bessel_J(0.1,0.1)
        0.777264368097005
    
    We check consistency of PARI and Maxima::
    
        sage: n(bessel_J(3,10,"maxima"))
        0.0583793793051...
        sage: n(bessel_J(3,10,"pari"))  
        0.0583793793051868
        sage: bessel_J(3,10,"scipy")
        0.0583793793052...

    Check whether the return value is real whenever the argument is real (#10251)::                                                                                                                                                           
        sage: bessel_J(5, 1.5, algorithm='scipy') in RR                                                                      
        True
    """
    
    if algorithm=="pari":
        from sage.libs.pari.all import pari
        try:
            R = RealField(prec)
            nu = R(nu)
            z = R(z)
        except TypeError:
            C = ComplexField(prec)
            nu = C(nu)
            z = C(z)
            K = C
        if nu == 0:
            nu = ZZ(0)
        K = z.parent()
        return K(pari(nu).besselj(z, precision=prec))
    elif algorithm=="scipy":
        if prec != 53:
            raise ValueError, "for the scipy algorithm the precision must be 53"
        import scipy.special
        ans = str(scipy.special.jv(float(nu),complex(real(z),imag(z))))
        ans = ans.replace("(","")
        ans = ans.replace(")","")
        ans = ans.replace("j","*I")
        ans = sage_eval(ans)
        return real(ans) if z in RR else ans
    elif algorithm == "maxima":
        if prec != 53:
            raise ValueError, "for the maxima algorithm the precision must be 53"
        return maxima_function("bessel_j")(nu, z)
    else:
        raise ValueError, "unknown algorithm '%s'"%algorithm

def bessel_I(nu,z,algorithm = "pari",prec=53):
    r"""
    Implements the "I-Bessel function", or "modified Bessel function,
    1st kind", with index (or "order") nu and argument z.
    
    INPUT:
    
    
    -  ``nu`` - a real (or complex, for pari) number
    
    -  ``z`` - a real (positive) algorithm - "pari" or
       "maxima" or "scipy" prec - real precision (for PARI only)
    
    
    DEFINITION::
    
            Maxima:
                             inf
                            ====   - nu - 2 k  nu + 2 k
                            \     2          z
                             >    -------------------
                            /     k! Gamma(nu + k + 1)
                            ====
                            k = 0
        
            PARI:
            
                             inf
                            ====   - 2 k  2 k
                            \     2      z    Gamma(nu + 1)
                             >    -----------------------
                            /       k! Gamma(nu + k + 1)
                            ====
                            k = 0
        
            
    
    Sometimes ``bessel_I(nu,z)`` is denoted
    ``I_nu(z)`` in the literature.
    
    .. warning::

       In Maxima (the manual says) i0 is deprecated but
       ``bessel_i(0,*)`` is broken. (Was fixed in recent CVS patch
       though.)
    
    EXAMPLES::
    
        sage: bessel_I(1,1,"pari",500)
        0.565159103992485027207696027609863307328899621621092009480294489479255640964371134092664997766814410064677886055526302676857637684917179812041131208121
        sage: bessel_I(1,1)
        0.565159103992485
        sage: bessel_I(2,1.1,"maxima")  
        0.16708949925104...
        sage: bessel_I(0,1.1,"maxima") 
        1.32616018371265...
        sage: bessel_I(0,1,"maxima")   
        1.2660658777520...
        sage: bessel_I(1,1,"scipy")
        0.565159103992...

    Check whether the return value is real whenever the argument is real (#10251)::
    
        sage: bessel_I(5, 1.5, algorithm='scipy') in RR
        True
        
    """
    if algorithm=="pari":
        from sage.libs.pari.all import pari
        try:
            R = RealField(prec)
            nu = R(nu)
            z = R(z)
        except TypeError:
            C = ComplexField(prec)
            nu = C(nu)
            z = C(z)
            K = C
        K = z.parent()
        return K(pari(nu).besseli(z, precision=prec))
    elif algorithm=="scipy":
        if prec != 53:
            raise ValueError, "for the scipy algorithm the precision must be 53"
        import scipy.special
        ans = str(scipy.special.iv(float(nu),complex(real(z),imag(z))))
        ans = ans.replace("(","")
        ans = ans.replace(")","")
        ans = ans.replace("j","*I")
        ans = sage_eval(ans)
        return real(ans) if z in RR else ans # Return real value when arg is real
    elif algorithm == "maxima":
        if prec != 53:
            raise ValueError, "for the maxima algorithm the precision must be 53"
        return sage_eval(maxima.eval("bessel_i(%s,%s)"%(float(nu),float(z))))
    else:
        raise ValueError, "unknown algorithm '%s'"%algorithm