CoCalc Public Fileswww / talks / bernoulli / misc.sage
Author: William A. Stein
1
2def 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
66def 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
130def 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