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?

Answer 1

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.