Opened 5 years ago

Closed 2 years ago

#21439 closed defect (fixed)

Inaccurate Bernoulli numbers from Flint

Reported by: rws Owned by:
Priority: major Milestone: sage-8.8
Component: numerical Keywords:
Cc: leif 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 vdelecroix)

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 5 years ago by rws

  • Description modified (diff)

comment:2 Changed 5 years ago by rws

  • Component changed from number theory to numerical

comment:3 Changed 5 years ago by leif

  • Cc leif 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 follow-up: Changed 5 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 2 years ago by vdelecroix

  • Authors set to Vincent Delecroix
  • Description modified (diff)
  • Milestone changed from sage-7.4 to sage-8.8

comment:6 Changed 2 years ago by vdelecroix

  • Description modified (diff)

comment:7 Changed 2 years ago by vdelecroix

  • Branch set to u/vdelecroix/21439
  • Commit set to 98b9df59f4ce32d74591f2a172d335842379b4bd
  • Status changed from new to needs_review

New commits:

98b9df5fix Bernoulli numbers

comment:8 in reply to: ↑ 4 Changed 2 years ago by vdelecroix

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 2 years ago by tscrim

  • Reviewers set to Travis Scrimshaw
  • Status changed from needs_review to positive_review

LGTM.

comment:10 Changed 2 years ago by vdelecroix

thanks

comment:11 Changed 2 years ago by vbraun

  • Branch changed from u/vdelecroix/21439 to 98b9df59f4ce32d74591f2a172d335842379b4bd
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.