def q_zeta(m, N):
return 1/prod([1-1/(p^m) for p in primes(N)])
def r_zeta(m, prec):
w = RealField(prec)(m)
return w.zeta()
def bernoulli1(m):
r"""
Returns the Bernoulli number $B_m$; do every detail directly.
"""
if m == 1:
return -1/2
if m % 2:
return 0
prec = 16*m
verbose('prec = %s'%prec)
tm = verbose('bernoulli_python setup')
R = RealField(prec)
tm = verbose('computing pi...', tm)
pi = R.pi()
tm = verbose('computing factorial...', tm)
m_factorial = factorial(m)
tm = verbose('computing pi pow...', tm)
K = 2*m_factorial/((2*pi)^m)
tm = verbose('computing P...', tm)
P = prime_range(m+2)
d = prod([p for p in P if m % (p-1) == 0])
tm = verbose('computing N...', tm)
dK = d*K
n = str(dK).find('.') + 1
N = float(10)^(n/(m-1))
tm = verbose('N = %s;\ncomputing product...'%N, tm)
assert N < m
z = 1
z_prev = z
Rl = RealField()
for p in P:
if p > N:
break
z /= 1 - 1/(R(p)^R(m))
diff = abs(z - z_prev).log()
print p, Rl(diff)
if diff < -prec:
break
z_prev =z
print z
a = long((dK*z).ceil())
if m % 4 == 0:
a = -a
return Rational(a)/Rational(d)
def bernoulli2(m):
r"""
Returns the Bernoulli number $B_m$. Use PARI and MPFR when
possible for parts of the calculation.
OBSERVATIONS:
* Pi:
MPFR computes pi *much* more quickly than PARI:
sage: time s = pari.new_with_bits_prec (pi, 10^5)
CPU times: user 2.48 s, sys: 0.01 s, total: 2.49 s
Wall time: 2.67
sage: time p = RealField(10^5)(pi)
CPU times: user 0.20 s, sys: 0.00 s, total: 0.20 s
* ZETA:
PARI computes zeta *much* more quickly than MPFR
sage: time m = zeta(RealField(10000)(1000))
CPU times: user 5.12 s, sys: 0.00 s, total: 5.13 s
Wall time: 5.22
sage: time n = pari.new_with_bits_prec(1000, 10000).zeta()
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 0.01
sage: time n=pari.new_with_bits_prec(1000, 20000).zeta()
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 0.06
sage: time m=zeta(RealField(20000)(1000))
CPU times: user 139.82 s, sys: 0.65 s, total: 140.46 s
Wall time: 164.22
* Factorials: My Pyrex compiled wrapping of GMP is best:
sage: time n = factorial(100000)
CPU times: user 6.93 s, sys: 0.00 s, total: 6.93 s
Wall time: 6.95
MAGMA:
sage: time n = magma('Factorial(100000)')
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 9.02
GP interpreter
sage: time gp('m=prod(n=1,100000,n);')
Wall time: 16.00
PARI C library:
sage: time n=pari('prod(n=1,100000,n)')
CPU times: user 17.64 s, sys: 0.22 s, total: 17.85 s
Wall time: 18.81
"""
if m == 1:
return -1/2
if m % 2:
return 0
t0 = verbose('computing bernoulli number %s'%m)
tm = verbose('computing d...')
P = prime_range(m+2)
d = prod([p for p in P if m % (p-1) == 0])
verbose('got d=%s'%d, tm)
tm = verbose('computing prec')
t = log(d) + (m + 0.5) * log(m) - m*(1+log(2*pi)) + 1.712086
u = t/(log(2)*4)
prec0 = int(u.ceil()) + 3
prec = 4*(int(u.ceil()) + 12)
verbose('prec = %s; prec0 = %s'%(prec,prec0), tm)
tm = verbose('bernoulli_python setup')
R = RealField(prec)
tm = verbose('computing pi...')
PI = R.pi()
verbose('got pi...', tm)
tm = verbose('computing factorial...')
m_factorial = factorial_gmp(m)
verbose('got factorial', tm)
tm = verbose('computing K ...')
K = (R(2*m_factorial))/((2*PI)^m)
verbose('got K', tm)
tm = verbose('computing zeta... (prec=%s)'%prec)
w = pari.new_with_prec(m, prec0).zeta()
verbose('got zeta...', tm)
tm = verbose("converting...")
z = R(w)
verbose('converted', tm)
tm = verbose('constructing bernoulli number')
dKz = d*K*z
a = long((dKz).round())
if m % 4 == 0:
a = -a
B = Rational(a)/Rational(d)
tm = verbose('done!', tm)
verbose('total time', t0)
return B
def ndigits(m):
"""
Compute the number of digits of the numerator of the m-th
Bernoulli number without actually computing the m-th Bernoulli
number.
"""
P = prime_range(m+2)
d = prod([p for p in P if m % (p-1) == 0])
x = log(2) + sum(math.log(n) for n in range(1, m+1)) \
- m*log(2) - m*log(pi) + log(d)
y = x/log(10)
return ZZ(int(y.ceil()))