mpf_bernoulli sometimes results sometimes non-normalized 

For example:
```pycon
>>> from mpmath import libmp
>>> libmp.mpf_bernoulli(10, 336, rnd='d')
(0, mpz(5970000938847648365496824257597517243193899889464461616355802030946321480529880395668828220437108430425246665198001), -385, 382)
>>> libmp.mpf_bernoulli(10, 336, rnd='u')
(0, mpz(84838815991583492824145540930622889322789602167979149049242011400143631484097005506851420199682087719), -339, 336)
```
while
```pycon
>>> from mpmath import libmp
>>> libmp.mpf_bernoulli(10, 336, rnd='u')
(0, mpz(5970000938847648365496824257597517243193899889464461616355802030946321480529880395668828220437108430425246665198001), -385, 382)
```

It seems that main loop assumes default `rnd=round_fast` value:
https://github.com/mpmath/mpmath/blob/9ea52503dd085ff8a15e17b1fc4f2d0cc36c2bb1/mpmath/libmp/gammazeta.py#L418-L450

	while m <= n:
	#print m
	case = m % 6
	# Accurately estimate size of B_m so we can use
	# fixed point math without using too much precision
	szbm = bernoulli_size(m)
	s = MPZ(0)
	sexp = max(0, szbm) - wp
	if m < 6:
	a = MPZ_ZERO
	else:
	a = bin1
	for j in range(1, m//6+1):
	usign, uman, uexp, ubc = u = numbers[m-6*j]
	if usign:
	uman = -uman
	s += lshift(a*uman, uexp-sexp)
	# Update inner binomial coefficient
	j6 = 6*j
	a = ((m-5-j6)(m-4-j6)(m-3-j6)(m-2-j6)(m-1-j6)(m-j6))
	a //= ((4+j6)(5+j6)(6+j6)(7+j6)(8+j6)*(9+j6))
	if case == 0: b = mpf_rdiv_int(m+3, f3, wp)
	if case == 2: b = mpf_rdiv_int(m+3, f3, wp)
	if case == 4: b = mpf_rdiv_int(-m-3, f6, wp)
	s = from_man_exp(s, sexp, wp)
	b = mpf_div(mpf_sub(b, s, wp), from_int(bin), wp)
	numbers[m] = b
	m += 2
	# Update outer binomial coefficient
	bin = bin * ((m+2)(m+3)) // (m(m-1))
	if m > 6:
	bin1 = bin1 * ((2+m)(3+m)) // ((m-7)(m-6))
	state[:] = [m, bin, bin1]

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

mpf_bernoulli sometimes results sometimes non-normalized #928

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

mpf_bernoulli sometimes results sometimes non-normalized #928

Description

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions