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