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: sage-6.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 zimmerma)

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.

Change History (12)

comment:1 Changed 6 years ago by jdemeyer

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

comment:2 Changed 6 years ago by zimmerma

  • 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 jdemeyer

mpfr is also much faster...

comment:4 Changed 6 years ago by jdemeyer

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 zimmerma

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 jdemeyer

  • 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 jdemeyer

  • 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 jdemeyer

  • Authors set to Jeroen Demeyer
  • Commit set to bce4fb0c1665025c98ec15d89c1345a0a67e925c
  • Status changed from new to needs_review

New commits:

6d10729Simplify numerical evaluation of BuiltinFunctions
b6e1ed4Merge remote-tracking branches 'origin/u/jdemeyer/ticket/17131' and 'origin/u/jdemeyer/ticket/17133' into ticket/17130
382f97aMerge branch 'u/jdemeyer/ticket/17130' of trac.sagemath.org:sage into 6.5beta1
726598917130: reviewer's patch: fix typo
c47dbd4Fix Trac #17328 again in a better way
a486db2Call the factorial() method of an Integer
bce4fb0Use mpfr for Bessel functions if possible

comment:9 Changed 6 years ago by zimmerma

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 jdemeyer

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 zimmerma

  • 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 vbraun

  • Branch changed from u/jdemeyer/ticket/17122 to bce4fb0c1665025c98ec15d89c1345a0a67e925c
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.