Opened 6 years ago
Closed 5 years ago
#17122 closed defect (fixed)
bessel_Y is off by 3 ulps
Reported by:  zimmerma  Owned by:  

Priority:  major  Milestone:  sage6.4 
Component:  basic arithmetic  Keywords:  
Cc:  Merged in:  
Authors:  Jeroen Demeyer  Reviewers:  Paul Zimmermann 
Report Upstream:  N/A  Work issues:  
Branch:  bce4fb0 (Commits)  Commit:  bce4fb0c1665025c98ec15d89c1345a0a67e925c 
Dependencies:  #17130  Stopgaps: 
Description (last modified by )
consider the following with Sage 6.0:
sage: R=RealField(113) sage: a=R("8.935761195587725798762818805462843676e01") sage: b=bessel_Y(0,a) sage: c=R(bessel_Y(0,RealField(200)(a))) sage: (bc)/c.ulp() 3.00000000000000000000000000000000 sage: b 7.44623881999333920107530266264974e7 sage: c 7.44623881999333920107530266264973e7
Given that MPFR provides correct rounding for bessel_Y (mpfr_y0) this should not happen.
Change History (12)
comment:1 Changed 6 years ago by
comment:2 Changed 6 years ago by
 Description modified (diff)
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).
comment:3 Changed 6 years ago by
mpfr is also much faster...
comment:4 Changed 6 years ago by
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.
comment:5 Changed 6 years ago by
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
.
comment:6 Changed 6 years ago by
 Dependencies set to #17130
I bumped into #17130 which I think should be fixed first (not by me).
comment:7 Changed 6 years ago by
 Branch set to u/jdemeyer/ticket/17122
 Created changed from 10/09/14 14:58:37 to 10/09/14 14:58:37
 Modified changed from 10/10/14 09:26:23 to 10/10/14 09:26:23
comment:8 Changed 6 years ago by
 Commit set to bce4fb0c1665025c98ec15d89c1345a0a67e925c
 Status changed from new to needs_review
New commits:
6d10729  Simplify numerical evaluation of BuiltinFunctions

b6e1ed4  Merge remotetracking 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

comment:9 Changed 6 years ago by
I wish I could review that patch, but unfortunately since the change to git I don't know any more how to do...
Paul
comment:10 Changed 6 years ago by
I think it's sufficient to look at http://git.sagemath.org/sage.git/commit/?id=bce4fb0c1665025c98ec15d89c1345a0a67e925c (the other commits come from #17130).
comment:11 Changed 5 years ago by
 Reviewers set to Paul Zimmermann
 Status changed from needs_review to positive_review
this commit looks correct to me. Positive review.
Paul
comment:12 Changed 5 years ago by
 Branch changed from u/jdemeyer/ticket/17122 to bce4fb0c1665025c98ec15d89c1345a0a67e925c
 Resolution set to fixed
 Status changed from positive_review to closed
The example is wrong, I get
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 forn
.Note that you can access the mpfr functions directly using