bessel_Y is off by 3 ulps

[zimmermann6]

consider the following with Sage 6.0:

sage: R=RealField(113)
sage: a=R("8.935761195587725798762818805462843676e-01")
sage: b=bessel_Y(0,a)
sage: c=R(bessel_Y(0,RealField(200)(a)))
sage: (b-c)/c.ulp()
-3.00000000000000000000000000000000
sage: b
-7.44623881999333920107530266264974e-7
sage: c
-7.44623881999333920107530266264973e-7

Given that MPFR provides correct rounding for bessel_Y (mpfr_y0) this should not happen.

Depends on #17130

Component: basic arithmetic

Author: Jeroen Demeyer

Branch/Commit: bce4fb0

Reviewer: Paul Zimmermann

Issue created by migration from https://trac.sagemath.org/ticket/17122

[jdemeyer]

comment:1

The example is wrong, I get

sage: R = RealField(113)
sage: a = R("1.414213562373095048801688724209698177")
sage: b = bessel_Y(0,a)
sage: c = R(bessel_Y(0,RealField(200)(a)))
sage: (b-c)/c.ulp()
0.000000000000000000000000000000000
sage: b
0.344636931299712753154578621698097
sage: c
0.344636931299712753154578621698097

Anyway, numerical evaluation of Bessel functions is done using mpmath, not mpfr. I guess the reason is that mpfr can only compute Y(n,x) for integers n, while mpmath supports more general complex numbers for n.

Note that you can access the mpfr functions directly using

sage: R = RealField(113)
sage: a = R("1.414213562373095048801688724209698177")
sage: a.y0()
0.344636931299712753154578621698097

[zimmermann6]

Description changed:

--- 
+++ 
@@ -2,7 +2,7 @@
 
 ```
 sage: R=RealField(113)
-sage: a=R("1.414213562373095048801688724209698177")
+sage: a=R("8.935761195587725798762818805462843676e-01")
 sage: b=bessel_Y(0,a)
 sage: c=R(bessel_Y(0,RealField(200)(a)))
 sage: (b-c)/c.ulp()

[zimmermann6]

comment:2

sorry, the value of a was wrong, I changed it in the description.

Indeed MPFR can only handle integer n for Y(n,x), but it would better to call it in that case,
since it guarantees correct rounding (and thus numerical reproducibility).

[jdemeyer]

comment:3

mpfr is also much faster...

[jdemeyer]

comment:4

Is there an MPFR function which can check if an mpfr_t is an exact integer? Or which acts like mpfr_get_si but without rounding, just converting if the number already was an exact integer (and returning an error otherwise). Or alternatively, which acts like mpfr_get_z and returns a ternary value.

[zimmermann6]

comment:5

there is mpfr_integer_p, but you also want mpfr_fits_slong_p if you want to convert the (integer) argument to a long before giving it to say mpfr_yn.

[jdemeyer]

comment:6

I bumped into #17130 which I think should be fixed first (not by me).

[jdemeyer]

Dependencies: #17130

[jdemeyer]

Branch: u/jdemeyer/ticket/17122

[jdemeyer]

Author: Jeroen Demeyer

[jdemeyer]

Commit: bce4fb0