Use mpfr for Bessel functions if possible

Author: jdemeyer

src/sage/functions/bessel.py

@@ -308,10 +308,32 @@ def _evalf_(self, n, x, parent=None, algorithm=None):
         EXAMPLES::
 
             sage: bessel_J(0.0, 1.0)
-            0.765197686557966
+            0.765197686557967
             sage: bessel_J(0, 1).n(digits=20)
             0.76519768655796655145
+            sage: bessel_J(0.5, 1.5)
+            0.649838074753747
+
+        Check for correct rounding (:trac:`17122`)::
+
+            sage: R = RealField(113)
+            sage: a = R("8.935761195587725798762818805462843676e-01")
+            sage: aa = RealField(200)(a)
+            sage: for n in [-10..10]:
+            ....:     b = bessel_J(R(n), a)
+            ....:     bb = R(bessel_J(n, aa))
+            ....:     if b != bb:
+            ....:         print n, b-bb
         """
+        if parent is not None:
+            x = parent(x)
+
+        try:
+            return x.jn(Integer(n))
+        except Exception:
+            pass
+
+        n, x = get_coercion_model().canonical_coercion(n, x)
         import mpmath
         return mpmath_utils.call(mpmath.besselj, n, x, parent=parent)
 
@@ -457,11 +479,33 @@ def _evalf_(self, n, x, parent=None, algorithm=None):
         """
         EXAMPLES::
 
+            sage: bessel_Y(0.5, 1.5)
+            -0.0460831658930974
             sage: bessel_Y(1.0+2*I, 3.0+4*I)
             0.699410324467538 + 0.228917940896421*I
             sage: bessel_Y(0, 1).n(256)
             0.08825696421567695798292676602351516282781752309067554671104384761199978932351
+
+        Check for correct rounding (:trac:`17122`)::
+
+            sage: R = RealField(113)
+            sage: a = R("8.935761195587725798762818805462843676e-01")
+            sage: aa = RealField(200)(a)
+            sage: for n in [-10..10]:
+            ....:     b = bessel_Y(R(n), a)
+            ....:     bb = R(bessel_Y(n, aa))
+            ....:     if b != bb:
+            ....:         print n, b-bb
         """
+        if parent is not None:
+            x = parent(x)
+
+        try:
+            return x.yn(Integer(n))
+        except Exception:
+            pass
+
+        n, x = get_coercion_model().canonical_coercion(n, x)
         import mpmath
         return mpmath_utils.call(mpmath.bessely, n, x, parent=parent)
 

src/sage/libs/mpmath/ext_impl.pyx

@@ -25,10 +25,8 @@ from cpython.float cimport *
 from cpython.complex cimport *
 from cpython.number cimport *
 
-cdef extern from "math.h":
-    cdef double fpow "pow" (double, double)
-    cdef double fsqrt "sqrt" (double)
-    cdef double frexp "frexp" (double, int*)
+from libc.math cimport sqrt as fsqrt
+from libc.math cimport frexp
 
 from sage.libs.gmp.all cimport *
 from sage.libs.mpfr cimport *
@@ -1147,6 +1145,7 @@ cdef cy_exp_basecase(mpz_t y, mpz_t x, int prec):
     mpz_fdiv_q_2exp(x2, x2, prec)
     mpz_set(a, x2)
     while mpz_sgn(a):
+        sig_check()
         mpz_fdiv_q_ui(a, a, k)
         mpz_add(s0, s0, a)
         k += 1
@@ -1160,6 +1159,7 @@ cdef cy_exp_basecase(mpz_t y, mpz_t x, int prec):
     mpz_add(s0, s0, s1)
     u = r
     while r:
+        sig_check()
         mpz_mul(s0, s0, s0)
         mpz_fdiv_q_2exp(s0, s0, prec)
         r -= 1
@@ -2037,6 +2037,7 @@ cdef MPF_hypsum(MPF *a, MPF *b, int p, int q, param_types, str ztype, coeffs, z,
     else:
         mpz_set_tuple_fixed(ZRE, z, wp)
     for i in range(0,p):
+        sig_check()
         if param_types[i] == 'Z':
             mpz_init(AINT[aint])
             mpz_set_integer(AINT[aint], coeffs[i])
@@ -2060,6 +2061,7 @@ cdef MPF_hypsum(MPF *a, MPF *b, int p, int q, param_types, str ztype, coeffs, z,
         else:
             raise ValueError
     for i in range(p,p+q):
+        sig_check()
         if param_types[i] == 'Z':
             mpz_init(BINT[bint])
             mpz_set_integer(BINT[bint], coeffs[i])
@@ -2113,22 +2115,28 @@ cdef MPF_hypsum(MPF *a, MPF *b, int p, int q, param_types, str ztype, coeffs, z,
 
         # Multiply real factors
         for k in range(0, cancellable_real):
+            sig_check()
             mpz_mul(PRE, PRE, AREAL[k])
             mpz_fdiv_q(PRE, PRE, BREAL[k])
         for k in range(cancellable_real, areal):
+            sig_check()
             mpz_mul(PRE, PRE, AREAL[k])
             mpz_fdiv_q_2exp(PRE, PRE, wp)
         for k in range(cancellable_real, breal):
+            sig_check()
             mpz_mul_2exp(PRE, PRE, wp)
             mpz_fdiv_q(PRE, PRE, BREAL[k])
         if have_complex:
             for k in range(0, cancellable_real):
+                sig_check()
                 mpz_mul(PIM, PIM, AREAL[k])
                 mpz_fdiv_q(PIM, PIM, BREAL[k])
             for k in range(cancellable_real, areal):
+                sig_check()
                 mpz_mul(PIM, PIM, AREAL[k])
                 mpz_fdiv_q_2exp(PIM, PIM, wp)
             for k in range(cancellable_real, breal):
+                sig_check()
                 mpz_mul_2exp(PIM, PIM, wp)
                 mpz_fdiv_q(PIM, PIM, BREAL[k])
 
@@ -2162,6 +2170,7 @@ cdef MPF_hypsum(MPF *a, MPF *b, int p, int q, param_types, str ztype, coeffs, z,
                 mpz_fdiv_q(PIM, PIM, DIV)
 
             for i in range(acomplex):
+                sig_check()
                 mpz_mul(TRE, PRE, ACRE[i])
                 mpz_submul(TRE, PIM, ACIM[i])
                 mpz_mul(TIM, PIM, ACRE[i])
@@ -2170,6 +2179,7 @@ cdef MPF_hypsum(MPF *a, MPF *b, int p, int q, param_types, str ztype, coeffs, z,
                 mpz_fdiv_q_2exp(PIM, TIM, wp)
 
             for i in range(bcomplex):
+                sig_check()
                 mpz_mul(URE, BCRE[i], BCRE[i])
                 mpz_addmul(URE, BCIM[i], BCIM[i])
                 mpz_mul(TRE, PRE, BCRE[i])

src/sage/rings/real_mpfr.pyx

@@ -3769,7 +3769,7 @@ cdef class RealNumber(sage.structure.element.RingElement):
 
     def is_integer(self):
         """
-        Return ``True`` if this number is a integer
+        Return ``True`` if this number is a integer.
 
         EXAMPLES::
         
@@ -3778,7 +3778,7 @@ cdef class RealNumber(sage.structure.element.RingElement):
             sage: RR(0.1).is_integer()
             False
         """
-        return self in ZZ
+        return mpfr_integer_p(self.value) != 0
 
     def __nonzero__(self):
         """