Compute Bernoulli numbers with Python recursive program

Question

I'm trying to solve a problem about Bernoulli numbers using Python. The aim is to output the numerator and the denominator of the $n$-th Bernoulli number. I use the conventions and the generic formula given in this source.

Here is my code. I use the auxiliary function aux_bernoulli to compute Bernoulli numbers using recursivity.

from fractions import Fraction
from math import factorial

def aux_bernoulli(n):
    if n == 0:
        return 1
    elif n == 1: # convention
        return -0.5
    elif (n-1)%2==0: # B(n)=0 when n is odd
        return 0
    else:
        somme = 0
        for k in range(n):
            somme += (factorial(n)/(factorial(n+1-k)*factorial(k))) * aux_bernoulli(k)
        return -somme

def bernoulli(n):
    ber = aux_bernoulli(n)
    print(ber) # for debugging purposes
    numerator, denominator = Fraction(ber).numerator, Fraction(ber).denominator
    return numerator, denominator

This code is giving me wrong values that are very close to the right ones and I can't understand figure out why. Here are some examples:

bernoulli(4)
bernoulli(6)
bernoulli(8)

Output:

-0.03333333333333338
(-600479950316067, 18014398509481984)

0.023809523809524058
(214457125112883, 9007199254740992)

-0.033333333333335075
(-1200959900632195, 36028797018963968)

Correct values according to this source:

-0.033333
(-1, 30)

0.0280952
(1/42)

-0.033333
(-1, 30)

Does anyone know what's wrong with my approach?

I suggest using `from math import comb` and then `comb(n, k) / (n+1-k)` rather than `(factorial(n)/(factorial(n+1-k)*factorial(k)))` — Stef, Jan 31 '23 at 12:36
If you do the computations with `float` but then convert to `Fraction`, I suggest using `.limit_denominator()` to simplify the fractions. — Stef, Jan 31 '23 at 12:38
Interesting reading: [Is floating-point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) and [What Every Computer Scientist Should Know About Floating-Point Arithmetic](http://download.oracle.com/docs/cd/E19957-01/806-3568/ncg_goldberg.html) — Stef, Jan 31 '23 at 12:43
Thanks a lot for the advice and the resources! Seems like it gives the desired output now. — Barbara Gendron, Jan 31 '23 at 12:48
Also note that your recursive function is somewhat inefficient because it recomputes the same values over and over again. The recursive calls are branching exponentially. For instance, `bernouilli(50)` will make recursive calls to all the `bernouilli(k)` for `k` between 0 and 49; and `bernouilli(49)` will again make recursive calls to all the `bernouilli(k)` for k between 0 and 48, even though those have already been calculated by the recursive calls from `bernouilli(50)`; etc — Stef, Jan 31 '23 at 12:48
To avoid these recomputations, I suggest filling an array with the computed values so that you don't have to recompute them. — Stef, Jan 31 '23 at 12:49
If you want to return a `Fraction`, I suggest doing all the computations with `Fraction` and never using `float`. This will avoid all approximations. [Try it online!](https://tio.run/##pVDBToQwFLz3K15iTFoXlLK6MSZ49MgPbIgpS3Eb4LUpJdGvx5ZdAc3qxTn1zZu@mYz5cEeN20djx7G2uoNOuCOozmjr4KC7kkxsbcXBKY391@rlTBBSyRrE8P5aSot6UG2raM6eCHioGnJ4ziAFgZV/XvtXlgE/bQOsdIPF@RpN2LQqIYP9THIWLYqYR5CyYpLV2gKCQrAC3yRNI8g3nC3XvT9eMJ0cboUxEiu6smazQra9/EX/jQ2IoR86GqqiGEHD4AbKfVPAHVDc8NgTIWazxMSVUcBpOjdREmKsQkd/VMpDvKt1KRFc6GVF@Hm3npMLH7bJXwovuE//d@LBh9ixYhw/AQ) — Stef, Jan 31 '23 at 12:57

score 0 · Accepted Answer · answered Jan 31 '23 at 18:00

Combining @Stef's various suggestions (multiple +1s), I came up with the following simplification:

from math import comb
from fractions import Fraction
from functools import lru_cache

@lru_cache
def bernoulli(n):
    if n == 0:
        return Fraction(1)

    if n == 1: # convention
        return Fraction(-1/2)

    somme = Fraction(0)

    if n % 2:  # B(n) = 0 when n is odd
        return somme

    for k in range(n):
        somme += bernoulli(k) * comb(n, k) / (n + 1 - k)

    return -somme

print(bernoulli(60).as_integer_ratio())

It's easy to mess up the result by moving between Fraction and float.

Compute Bernoulli numbers with Python recursive program

1 Answers1