Sharedwww / talks / bernoulli / bern.sageOpen in CoCalc
Author: William A. Stein
1
def q_zeta(m, N):
2
return 1/prod([1-1/(p^m) for p in primes(N)])
3
4
def r_zeta(m, prec):
5
w = RealField(prec)(m)
6
return w.zeta()
7
8
def bernoulli1(m):
9
r"""
10
Returns the Bernoulli number $B_m$; do every detail directly.
11
"""
12
if m == 1:
13
return -1/2
14
if m % 2:
15
return 0
16
17
prec = 16*m
18
verbose('prec = %s'%prec)
19
tm = verbose('bernoulli_python setup')
20
R = RealField(prec)
21
tm = verbose('computing pi...', tm)
22
pi = R.pi()
23
tm = verbose('computing factorial...', tm)
24
m_factorial = factorial(m)
25
tm = verbose('computing pi pow...', tm)
26
K = 2*m_factorial/((2*pi)^m)
27
tm = verbose('computing P...', tm)
28
P = prime_range(m+2)
29
30
# IDIOTIC -- should compute by factoring m!
31
d = prod([p for p in P if m % (p-1) == 0])
32
33
tm = verbose('computing N...', tm)
34
dK = d*K
35
n = str(dK).find('.') + 1
36
N = float(10)^(n/(m-1))
37
tm = verbose('N = %s;\ncomputing product...'%N, tm)
38
assert N < m
39
z = 1
40
z_prev = z
41
Rl = RealField()
42
for p in P:
43
if p > N:
44
break
45
z /= 1 - 1/(R(p)^R(m))
46
diff = abs(z - z_prev).log()
47
print p, Rl(diff)
48
if diff < -prec:
49
break
50
z_prev =z
51
52
print z
53
a = long((dK*z).ceil())
54
if m % 4 == 0:
55
a = -a
56
return Rational(a)/Rational(d)
57
#return K
58
59
def bernoulli2(m):
60
r"""
61
Returns the Bernoulli number $B_m$. Use PARI and MPFR when
62
possible for parts of the calculation.
63
64
OBSERVATIONS:
65
* Pi:
66
MPFR computes pi *much* more quickly than PARI:
67
sage: time s = pari.new_with_bits_prec (pi, 10^5)
68
CPU times: user 2.48 s, sys: 0.01 s, total: 2.49 s
69
Wall time: 2.67
70
sage: time p = RealField(10^5)(pi)
71
CPU times: user 0.20 s, sys: 0.00 s, total: 0.20 s
72
73
* ZETA:
74
PARI computes zeta *much* more quickly than MPFR
75
sage: time m = zeta(RealField(10000)(1000))
76
CPU times: user 5.12 s, sys: 0.00 s, total: 5.13 s
77
Wall time: 5.22
78
sage: time n = pari.new_with_bits_prec(1000, 10000).zeta()
79
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
80
Wall time: 0.01
81
sage: time n=pari.new_with_bits_prec(1000, 20000).zeta()
82
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
83
Wall time: 0.06
84
sage: time m=zeta(RealField(20000)(1000))
85
CPU times: user 139.82 s, sys: 0.65 s, total: 140.46 s
86
Wall time: 164.22
87
88
* Factorials: My Pyrex compiled wrapping of GMP is best:
89
sage: time n = factorial(100000)
90
CPU times: user 6.93 s, sys: 0.00 s, total: 6.93 s
91
Wall time: 6.95
92
93
MAGMA:
94
sage: time n = magma('Factorial(100000)')
95
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
96
Wall time: 9.02
97
98
GP interpreter
99
sage: time gp('m=prod(n=1,100000,n);')
100
Wall time: 16.00
101
102
PARI C library:
103
sage: time n=pari('prod(n=1,100000,n)')
104
CPU times: user 17.64 s, sys: 0.22 s, total: 17.85 s
105
Wall time: 18.81
106
"""
107
if m == 1:
108
return -1/2
109
if m % 2:
110
return 0
111
112
t0 = verbose('computing bernoulli number %s'%m)
113
114
tm = verbose('computing d...')
115
P = prime_range(m+2)
116
# IDIOTIC -- should compute by factoring m!
117
d = prod([p for p in P if m % (p-1) == 0])
118
verbose('got d=%s'%d, tm)
119
120
tm = verbose('computing prec')
121
t = log(d) + (m + 0.5) * log(m) - m*(1+log(2*pi)) + 1.712086
122
u = t/(log(2)*4)
123
prec0 = int(u.ceil()) + 3
124
prec = 4*(int(u.ceil()) + 12)
125
verbose('prec = %s; prec0 = %s'%(prec,prec0), tm)
126
127
#prec = 16*m
128
#verbose('prec = %s'%prec)
129
130
tm = verbose('bernoulli_python setup')
131
R = RealField(prec)
132
tm = verbose('computing pi...')
133
PI = R.pi()
134
verbose('got pi...', tm)
135
136
tm = verbose('computing factorial...')
137
m_factorial = factorial_gmp(m)
138
#m_factorial = R.factorial(m)
139
verbose('got factorial', tm)
140
141
tm = verbose('computing K ...')
142
K = (R(2*m_factorial))/((2*PI)^m)
143
verbose('got K', tm)
144
145
146
tm = verbose('computing zeta... (prec=%s)'%prec)
147
#w = pari.new_with_bits_prec(m, prec).zeta()
148
w = pari.new_with_prec(m, prec0).zeta()
149
verbose('got zeta...', tm)
150
#z = R(w)
151
tm = verbose("converting...")
152
z = R(w)
153
verbose('converted', tm)
154
155
tm = verbose('constructing bernoulli number')
156
dKz = d*K*z
157
#print 'frac = ', str(dKz.frac())[:30]
158
a = long((dKz).round())
159
if m % 4 == 0:
160
a = -a
161
B = Rational(a)/Rational(d)
162
tm = verbose('done!', tm)
163
verbose('total time', t0)
164
return B
165
166
167
168
def ndigits(m):
169
"""
170
Compute the number of digits of the numerator of the m-th
171
Bernoulli number without actually computing the m-th Bernoulli
172
number.
173
"""
174
P = prime_range(m+2)
175
d = prod([p for p in P if m % (p-1) == 0])
176
x = log(2) + sum(math.log(n) for n in range(1, m+1)) \
177
- m*log(2) - m*log(pi) + log(d)
178
y = x/log(10)
179
return ZZ(int(y.ceil()))
180
181
182