Download code

Jump to: navigation, search

Back to Bernoulli_numbers_(Python)

Download for Windows: zip

Download for UNIX: zip, tar.gz, tar.bz2

bernoulli_simple.py

 1 # The authors of this work have released all rights to it and placed it
 2 # in the public domain under the Creative Commons CC0 1.0 waiver
 3 # (http://creativecommons.org/publicdomain/zero/1.0/).
 4 # 
 5 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 6 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 7 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 8 # IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 9 # CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12 # 
13 # Retrieved from: http://en.literateprograms.org/Bernoulli_numbers_(Python)?oldid=5312
14 
15 from gmpy import mpq, bincoef as _bincoef
16 
17 def bincoef(n, k):
18     if k < 0:
19         return 0
20     else:
21         return _bincoef(n, k)
22 
23 class memo:
24     def __init__(self, f):
25         self.func = f
26         self.cache = {}
27         self.highest = 0
28     def __call__(self, m):
29         if m < self.highest:
30             return self.cache[m]
31         else:
32             x = self.func(m)
33             self.cache[m] = x
34             self.highest = m
35             return x
36 
37 @memo
38 def bernoulli(m):
39     if m == 0: return 1
40     if m == 1: return -mpq(1,2)
41     if m % 2: return 0
42     s = 0
43     for j in range(0, m):
44         s += bincoef(m+1, j) * bernoulli(j)
45     b = -mpq(s)/bincoef(m+1, m)
46     return b


hijacker
hijacker
hijacker
hijacker

bernoulli_faster_yet.py

 1 # The authors of this work have released all rights to it and placed it
 2 # in the public domain under the Creative Commons CC0 1.0 waiver
 3 # (http://creativecommons.org/publicdomain/zero/1.0/).
 4 # 
 5 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 6 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 7 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 8 # IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 9 # CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12 # 
13 # Retrieved from: http://en.literateprograms.org/Bernoulli_numbers_(Python)?oldid=5312
14 
15 from gmpy import mpq, bincoef as _bincoef
16 
17 def bincoef(n, k):
18     if k < 0:
19         return 0
20     else:
21         return _bincoef(n, k)
22 
23 class memo:
24     def __init__(self, f):
25         self.func = f
26         self.cache = {}
27         self.highest = 0
28     def __call__(self, m):
29         if m < self.highest:
30             return self.cache[m]
31         else:
32             x = self.func(m)
33             self.cache[m] = x
34             self.highest = m
35             return x
36 def product(low, high):
37     p = 1
38     for k in range(low, high+1):
39         p *= k
40     return p
41 
42 def A(m, mtop):
43     s = 0
44     a = bincoef(m+3, m-6)
45     for j in range(1, mtop+1):
46         s += a*bernoulli(m-6*j)
47         a *= product(m-6 - 6*j + 1, m-6*j)
48         a //= product(6*j+4, 6*j+9)
49     return s
50 @memo
51 def b0mod6(m):
52     return (mpq(m+3, 3) - A(m, m//6)) / bincoef(m+3, m)
53 
54 @memo
55 def b2mod6(m):
56     return (mpq(m+3, 3) - A(m, (m-2)//6)) / bincoef(m+3, m)
57 
58 @memo
59 def b4mod6(m):
60     return (-mpq(m+3, 6) - A(m, (m-4)//6)) / bincoef(m+3, m)
61 def bernoulli(m):
62     assert m >= 0
63     if m == 0: return 1
64     if m == 1: return -mpq(1,2)
65     if m % 6 == 0: return b0mod6(m)
66     if m % 6 == 2: return b2mod6(m)
67     if m % 6 == 4: return b4mod6(m)
68     return 0


hijacker
hijacker
hijacker
hijacker

bernoulli_faster.py

 1 # The authors of this work have released all rights to it and placed it
 2 # in the public domain under the Creative Commons CC0 1.0 waiver
 3 # (http://creativecommons.org/publicdomain/zero/1.0/).
 4 # 
 5 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 6 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 7 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 8 # IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 9 # CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12 # 
13 # Retrieved from: http://en.literateprograms.org/Bernoulli_numbers_(Python)?oldid=5312
14 
15 from gmpy import mpq, bincoef as _bincoef
16 
17 def bincoef(n, k):
18     if k < 0:
19         return 0
20     else:
21         return _bincoef(n, k)
22 
23 class memo:
24     def __init__(self, f):
25         self.func = f
26         self.cache = {}
27         self.highest = 0
28     def __call__(self, m):
29         if m < self.highest:
30             return self.cache[m]
31         else:
32             x = self.func(m)
33             self.cache[m] = x
34             self.highest = m
35             return x
36 def A(m, M):
37     s = 0
38     for j in range(1, M+1):
39         s += bincoef(m+3, m-6*j) * bernoulli(m-6*j)
40     return s
41 @memo
42 def b0mod6(m):
43     return (mpq(m+3, 3) - A(m, m//6)) / bincoef(m+3, m)
44 
45 @memo
46 def b2mod6(m):
47     return (mpq(m+3, 3) - A(m, (m-2)//6)) / bincoef(m+3, m)
48 
49 @memo
50 def b4mod6(m):
51     return (-mpq(m+3, 6) - A(m, (m-4)//6)) / bincoef(m+3, m)
52 def bernoulli(m):
53     assert m >= 0
54     if m == 0: return 1
55     if m == 1: return -mpq(1,2)
56     if m % 6 == 0: return b0mod6(m)
57     if m % 6 == 2: return b2mod6(m)
58     if m % 6 == 4: return b4mod6(m)
59     return 0


hijacker
hijacker
hijacker
hijacker

test.py

 1 # The authors of this work have released all rights to it and placed it
 2 # in the public domain under the Creative Commons CC0 1.0 waiver
 3 # (http://creativecommons.org/publicdomain/zero/1.0/).
 4 # 
 5 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 6 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 7 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 8 # IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 9 # CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12 # 
13 # Retrieved from: http://en.literateprograms.org/Bernoulli_numbers_(Python)?oldid=5312
14 
15 from bernoulli_simple import bernoulli as bernoulli1
16 from bernoulli_faster import bernoulli as bernoulli2
17 from bernoulli_faster_yet import bernoulli as bernoulli3
18 from time import clock
19 
20 for i in range(50):
21     assert bernoulli1(i) == bernoulli2(i) == bernoulli3(i)
22 
23 a = clock()
24 temp = bernoulli1(1000)
25 print "simple:", clock() - a, "seconds"
26 
27 a = clock()
28 temp = bernoulli2(1000)
29 print "faster:", clock() - a, "seconds"
30 
31 a = clock()
32 temp = bernoulli3(1000)
33 print "faster yet:", clock() - a, "seconds"
34 


hijacker
hijacker
hijacker
hijacker