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