-
Notifications
You must be signed in to change notification settings - Fork 1
/
b.py
71 lines (57 loc) · 2 KB
/
b.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
# The implementation of the zeta-function algorithm for computing Bernoulli numbers quickly and highly accurate.
#
# Computational details:
# 1. Sieve of Eratosthenes for the prime numbers.
# 2. Computation of the zeta-function over primes.
# 3. The expression of the zeta-function and the Bernoulli numbers.
# 4. The Clausen and von Staudt theorem for computing denominator.
#
# Notice that b(1) = +0.5 which could be essential for the Euler-Maclaurin formula.
#
# Ref. Kevin J. McGown, Computing Bernoulli Numbers Quickly, 2005
from decimal import *
import math
pi = Decimal('3.141592653589793238462643383279502884197169399375105820974944592307816406286')
def factorial(n):return reduce(lambda a,b:a*b,[1]+range(1,n+1))
# Algorithm for finding all prime numbers up to the given limit known as Sieve of Eratosthenes
def prime_sieve(n):
prime_list = []
res = []
for i in range(2, n+1):
if i not in prime_list:
res.append(i)
for j in range(i*i, n+1, i):
prime_list.append(j)
return res
# The general procedure for the McGown algorithm
def b(n):
if n == 0:
return 1
elif n == 1:
return 0.5
elif (n-1)%2 == 0:
return 0
else:
K = 2*factorial(n)*1/Decimal((2*pi)**(n))
primes = prime_sieve(n+1)
d = Decimal(1)
for p in primes:
if n%(p-1)==0:
d *= p
N = math.ceil((K*d)**(Decimal(1.0/(n-1))))
z = Decimal(1)
for p in primes:
if p <= N:
z *= 1/(1-1/(Decimal(p)**n))
a = (-1)**(n/2+1)*Decimal(d*K*z)
a = a.to_integral_exact(rounding=ROUND_HALF_EVEN)
if a == 0:
a = (-1)**(n/2+1)
print(a) # Print numerator
print(d) # Print denominator
return a/d
n = input()
# Set precision. It allows accurate cumputing up to the 154th Bernoulli number.
if n > 28:
getcontext().prec = n
print(b(n)) # Print n-th Bernoulli number