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]

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

[jdemeyer]

New commits:

6d10729	Simplify numerical evaluation of BuiltinFunctions
b6e1ed4	Merge remote-tracking branches 'origin/u/jdemeyer/ticket/17131' and 'origin/u/jdemeyer/ticket/17133' into ticket/17130
382f97a	Merge branch 'u/jdemeyer/ticket/17130' of trac.sagemath.org:sage into 6.5beta1
7265989	17130: reviewer's patch: fix typo
c47dbd4	Fix Trac #17328 again in a better way
a486db2	Call the factorial() method of an Integer
bce4fb0	Use mpfr for Bessel functions if possible

[zimmermann6]

comment:9

I wish I could review that patch, but unfortunately since the change to git I don't know any more how
to do...

Paul

[jdemeyer]

comment:10

I think it's sufficient to look at sagemath/sagetrac-mirror@bce4fb0 (the other commits come from #17130).

[zimmermann6]

Reviewer: Paul Zimmermann

[zimmermann6]

comment:11

this commit looks correct to me. Positive review.

Paul

[vbraun]

Changed branch from u/jdemeyer/ticket/17122 to bce4fb0