Opened 6 years ago

Closed 3 years ago

#21439 closed defect (fixed)

Inaccurate Bernoulli numbers from Flint

Reported by: Ralf Stephan Owned by:
Priority: major Milestone: sage-8.8
Component: numerical Keywords:
Cc: Leif Leonhardy Merged in:
Authors: Vincent Delecroix Reviewers: Travis Scrimshaw
Report Upstream: N/A Work issues:
Branch: 98b9df5 (Commits, GitHub, GitLab) Commit: 98b9df59f4ce32d74591f2a172d335842379b4bd
Dependencies: Stopgaps:

Status badges

Description (last modified by Vincent Delecroix)

tldr: Bernoulli(250408) has a million digit denominator whose value differs by 1428094259 from the correct value.

One has

sage: %time bflint = bernoulli(250408, algorithm='flint')
CPU times: user 9.86 s, sys: 49.9 ms, total: 9.91 s
Wall time: 9.95 s
sage: %time barb = bernoulli(250408, algorithm='arb')
CPU times: user 10.2 s, sys: 23.9 ms, total: 10.2 s
Wall time: 10.3 s
sage: %time bpari = bernoulli(250408, algorithm='pari')
CPU times: user 11.6 s, sys: 36.9 ms, total: 11.7 s
Wall time: 11.7 s
sage: %time bbernmm = bernoulli(250408, algorithm='bernmm')
CPU times: user 13.3 s, sys: 0 ns, total: 13.3 s
Wall time: 13.3 s

And pari/arb/bernmm agrees on the answer (and disagree with flint)

sage: bbernmm == barb == bpari
True
sage: barb == bflint
False

Fredrik Johansson advises to switch to Arb for computing Bernoulli numbers because of https://github.com/wbhart/flint2/issues/288

This ticket modifies the default so that flint is used in small range, arb in intermediate range and bernmm in large range. The choice algorithm=flint with large input generates a warning. With the branch applied

sage: bernoulli(250408) == bernoulli(250408, algorithm='bernmm')
True
sage: b = bernoulli(250408, algorithm='flint')
.../sage/arith/misc.py:368: UserWarning: flint is known to not be accurate for large Bernoulli numbers

Change History (11)

comment:1 Changed 6 years ago by Ralf Stephan

Description: modified (diff)

comment:2 Changed 6 years ago by Ralf Stephan

Component: number theorynumerical

comment:3 Changed 6 years ago by Leif Leonhardy

Cc: Leif Leonhardy added

We may temporarily just do what the reporter on sage-devel did, i.e., change the "algorithm" used (perhaps depending on the argument):

Since Sage uses flint as the default algorithm for n<=300000 the error
affects Sage and can be observed, e.g., by comparing
bernoulli(250408,'flint') and bernoulli(250408,'bernmm').

(Just noticed algorithm="arb" is also already supported; there are various backends.)

Or apply the upstream fix...

comment:4 Changed 6 years ago by wbhart

Note the flint fix does not guarantee anything. It should fix the problem for all practical sizes, but there are no proven bounds. Arb guarantees the correct result up to bugs.

comment:5 Changed 3 years ago by Vincent Delecroix

Authors: Vincent Delecroix
Description: modified (diff)
Milestone: sage-7.4sage-8.8

comment:6 Changed 3 years ago by Vincent Delecroix

Description: modified (diff)

comment:7 Changed 3 years ago by Vincent Delecroix

Branch: u/vdelecroix/21439
Commit: 98b9df59f4ce32d74591f2a172d335842379b4bd
Status: newneeds_review

New commits:

98b9df5fix Bernoulli numbers

comment:8 in reply to:  4 Changed 3 years ago by Vincent Delecroix

Replying to wbhart:

Note the flint fix does not guarantee anything. It should fix the problem for all practical sizes, but there are no proven bounds. Arb guarantees the correct result up to bugs.

All values n <= 20000 currently used in the thresholds have been tested one by one. So that there is a guarantee (until we upgrade flint :-).

comment:9 Changed 3 years ago by Travis Scrimshaw

Reviewers: Travis Scrimshaw
Status: needs_reviewpositive_review

LGTM.

comment:10 Changed 3 years ago by Vincent Delecroix

thanks

comment:11 Changed 3 years ago by Volker Braun

Branch: u/vdelecroix/2143998b9df59f4ce32d74591f2a172d335842379b4bd
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.