Sharedwww / talks / bernoulli / misc.sageOpen in CoCalc
Author: William A. Stein
1
2
def bernoulli_pari_ish(m):
3
r"""
4
Returns the Bernoulli number $B_m$. Use PARI and MPFR when
5
possible for parts of the calculation.
6
"""
7
if m == 1:
8
return -1/2
9
if m % 2:
10
return 0
11
12
t0 = verbose('computing bernoulli number %s'%m)
13
14
tm = verbose('computing d...')
15
P = prime_range(m+2)
16
# IDIOTIC -- should compute by factoring m!
17
d = prod([p for p in P if m % (p-1) == 0])
18
verbose('got d=%s'%d, tm)
19
20
tm = verbose('computing prec')
21
t = log(d) + (m + 0.5) * log(m) - m*(1+log(2*pi)) + 1.712086
22
u = t/(log(2)*4)
23
prec0 = int(u.ceil()) + 3
24
prec = 4*(int(u.ceil()) + 12)
25
verbose('prec = %s; prec0 = %s'%(prec,prec0), tm)
26
27
#prec = 16*m
28
#verbose('prec = %s'%prec)
29
30
tm = verbose('bernoulli_python setup')
31
R = RealField(prec)
32
tm = verbose('computing pi...')
33
PI = R.pi()
34
verbose('got pi...', tm)
35
tm = verbose('converting to PARI')
36
PI = pari(PI)
37
verbose('done', tm)
38
39
tm = verbose('computing factorial...')
40
m_factorial = pari(factorial_gmp(m))
41
verbose('got factorial', tm)
42
43
tm = verbose('computing K ...')
44
K = (pari(2)*m_factorial)/((pari(2)*PI)^pari(m))
45
verbose('got K', tm)
46
47
48
tm = verbose('computing zeta...')
49
z = pari.new_with_prec(m, prec0).zeta()
50
verbose('got zeta...', tm)
51
52
tm = verbose('constructing bernoulli number')
53
dKz = pari(d)*K*z
54
a = long((dKz).round())
55
if m % 4 == 0:
56
a = -a
57
B = Rational(a)/Rational(d)
58
tm = verbose('done!', tm)
59
verbose('total time', t0)
60
return B
61
62
63
64
65
66
def bernoulli_invzeta(m):
67
r"""
68
Returns the Bernoulli number $B_m$. Use PARI and MPFR when
69
possible for parts of the calculation.
70
"""
71
if m == 1:
72
return -1/2
73
if m % 2:
74
return 0
75
76
t0 = verbose('computing bernoulli number %s'%m)
77
78
tm = verbose('computing d...')
79
P = prime_range(m+2)
80
#d = prod([p_minus_1 + 1 for p_minus_1 in divisors(m) if is_prime(p_minus_1+1)])
81
d = 6*prod([2*d + 1 for d in divisors(m//2)[1:] if is_prime(2*d+1)])
82
tm = verbose('got d=%s'%d, tm)
83
d0 = prod([p for p in P if m % (p-1) == 0])
84
assert d==d0, 'd=%s, d0=%s'%(d,d0)
85
verbose('got d=%s old way'%d, tm)
86
87
tm = verbose('computing prec')
88
89
t = log(d) + (m + 0.5) * log(m) - m*(1+log(2*pi)) + 1.712086
90
u = t/(log(2)*4)
91
prec0 = int(u.ceil()) + 3
92
93
prec = 4*(int(u.ceil()) + 12)
94
verbose('prec = %s; prec0 = %s'%(prec,prec0), tm)
95
96
#prec = 16*m
97
#verbose('prec = %s'%prec)
98
99
tm = verbose('bernoulli_python setup')
100
R = RealField(prec)
101
tm = verbose('computing pi...')
102
PI = R.pi()
103
verbose('got pi...', tm)
104
105
tm = verbose('computing factorial...')
106
m_factorial = factorial_gmp(m)
107
#m_factorial = R.factorial(m)
108
verbose('got factorial', tm)
109
110
tm = verbose('computing K ...')
111
K = (R(2*m_factorial))/((2*PI)^m)
112
verbose('got K', tm)
113
114
115
tm = verbose('computing inv zeta...')
116
invz = inverse_zeta(m, t, prec)
117
verbose('got zeta...', tm)
118
119
tm = verbose('constructing bernoulli number')
120
dKz = d*K/invz
121
#print 'frac = ', str(dKz.frac())[:30]
122
a = long((dKz).round())
123
if m % 4 == 0:
124
a = -a
125
B = Rational(a)/Rational(d)
126
tm = verbose('done!', tm)
127
verbose('total time', t0)
128
return B
129
130
def inverse_zeta(n, t, prec):
131
D = exp((t - log(n-1)) / (n-1))
132
lim = 1 + long(D.ceil())
133
prec = prec + 1
134
R = RealField(int((prec)))
135
print lim, prec
136
n = R(n)
137
z = 1 - 1/R(2)^n
138
for p in prime_range(3, lim+1):
139
#print p
140
h = z/(R(p)^R(n))
141
z = z - h
142
return z
143