Inaccurate Bernoulli numbers from Flint

[Ralf Stephan]

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

[Leif Leonhardy]

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...

[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.

[Vincent Delecroix]

New commits:

​98b9df5

fix Bernoulli numbers

[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 :-).

[Travis Scrimshaw]

LGTM.

[Vincent Delecroix]

thanks