Skip to content

mpf_bernoulli sometimes results sometimes non-normalized  #928

Closed
#939
Closed
mpf_bernoulli sometimes results sometimes non-normalized #928
#939
Labels
bugan unexpected problem or unintended behavior
Milestone
1.4

Description

@skirpichev
skirpichev
opened on Mar 16, 2025

For example:

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

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

mpmath/mpmath/libmp/gammazeta.py

Lines 418 to 450 in 9ea5250

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]

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugan unexpected problem or unintended behavior

    Type

    No type

    Fields

    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions