Open in CoCalc
1def q_zeta(m, N):
2    return 1/prod([1-1/(p^m) for p in primes(N)])
3
4def r_zeta(m, prec):
5    w = RealField(prec)(m)
6    return w.zeta()
7
8def 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
59def 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
168def 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