■12月12日(日)11:00-12:00 桜井進のPython･UNIX・Math教室（応用コース）
超入門・Pythonで関・ベルヌーイ数
- My code
- Akiyama-Tanigawa algorithm
- B(n) is just sum of k^n formula linear term coefficient.
- Zeta function algorithm
# 1. My code
# seki-bernoulli_1.py
from fractions import Fraction
# 二項係数 Combination nCk
def comb(n, k):
prod = 1
for t in range(min(k, n-k)):
prod = prod * (n-t)//(t+1)
return prod
# Seki-Bernoulli number Bn
def B(n):
if n == 0:
return 1
else:
ss = 0
for k in range(0, n):
ss = ss + ((-1)**k) * comb(n+1, k)*B(k)
return Fraction((-1)**(n+1), n+1) * ss
n = int(input('B(n)のnの上限 >>> '))
import sympy
for i in range(n + 1):
print(f"B({i})= {B(i)}".ljust(20, " "),f"sympy.bernoulli({i})=",sympy.bernoulli(i))
実行結果
$ python seki-bernoulli_1.py
B(n)のnの上限 >>> 20
B(0)= 1 sympy.bernoulli(0)= 1
B(1)= 1/2 sympy.bernoulli(1)= -1/2
B(2)= 1/6 sympy.bernoulli(2)= 1/6
B(3)= 0 sympy.bernoulli(3)= 0
B(4)= -1/30 sympy.bernoulli(4)= -1/30
B(5)= 0 sympy.bernoulli(5)= 0
B(6)= 1/42 sympy.bernoulli(6)= 1/42
B(7)= 0 sympy.bernoulli(7)= 0
B(8)= -1/30 sympy.bernoulli(8)= -1/30
B(9)= 0 sympy.bernoulli(9)= 0
B(10)= 5/66 sympy.bernoulli(10)= 5/66
B(11)= 0 sympy.bernoulli(11)= 0
B(12)= -691/2730 sympy.bernoulli(12)= -691/2730
B(13)= 0 sympy.bernoulli(13)= 0
B(14)= 7/6 sympy.bernoulli(14)= 7/6
B(15)= 0 sympy.bernoulli(15)= 0
B(16)= -3617/510 sympy.bernoulli(16)= -3617/510
B(17)= 0 sympy.bernoulli(17)= 0
B(18)= 43867/798 sympy.bernoulli(18)= 43867/798
B(19)= 0 sympy.bernoulli(19)= 0
B(20)= -174611/330 sympy.bernoulli(20)= -174611/330
案の定、関数B(n)は再帰的定義をしているため、B(30)の計算はできません。
結果をsympy.bernoulli()と比較してみます。
このアルゴリズムはヤコブ・ベルヌーイによるべき乗和による導入によるもの
したがって、B(1)=\frac{1}{2}
sympy.bernoulli(1)=-1/2 は、現在主流である母関数\frac{x}{e^x-1}による関・ベルヌーイ数の定義
次はゼータの負の整数に対する公式によるもの
関数mpmath.zeta()を用いるだけのコード
# 4 zeta fuction algorithm
# seki-bernoulli_4.py
print('4 zeta fuction algorithm')
n = int(input('B(n)のnの上限 >>> '))
from mpmath import zeta
def B(m):
b = -m * zeta(1-m)
return b
from sympy import *
for i in range(1, n + 1):
b = B(i)
c = Rational(b)
print(f"B({i})= {c}".ljust(20, " "),f"sympy.bernoulli({i})=",bernoulli(i))
実行結果
$ python seki-bernoulli_4.py
4 zeta fuction algorithm
B(n)のnの上限 >>> 50
B(1)= 1/2 sympy.bernoulli(1)= -1/2
B(2)= 1/6 sympy.bernoulli(2)= 1/6
B(3)= 0 sympy.bernoulli(3)= 0
B(4)= -1/30 sympy.bernoulli(4)= -1/30
B(5)= 0 sympy.bernoulli(5)= 0
B(6)= 1/42 sympy.bernoulli(6)= 1/42
B(7)= 0 sympy.bernoulli(7)= 0
B(8)= -1/30 sympy.bernoulli(8)= -1/30
B(9)= 0 sympy.bernoulli(9)= 0
B(10)= 5/66 sympy.bernoulli(10)= 5/66
B(11)= 0 sympy.bernoulli(11)= 0
B(12)= -691/2730 sympy.bernoulli(12)= -691/2730
B(13)= 0 sympy.bernoulli(13)= 0
B(14)= 7/6 sympy.bernoulli(14)= 7/6
B(15)= 0 sympy.bernoulli(15)= 0
B(16)= -3617/510 sympy.bernoulli(16)= -3617/510
B(17)= 0 sympy.bernoulli(17)= 0
B(18)= 43867/798 sympy.bernoulli(18)= 43867/798
B(19)= 0 sympy.bernoulli(19)= 0
B(20)= -174611/330 sympy.bernoulli(20)= -174611/330
B(21)= 0 sympy.bernoulli(21)= 0
B(22)= 854513/138 sympy.bernoulli(22)= 854513/138
B(23)= 0 sympy.bernoulli(23)= 0
B(24)= -236364091/2730 sympy.bernoulli(24)= -236364091/2730
B(25)= 0 sympy.bernoulli(25)= 0
B(26)= 8553103/6 sympy.bernoulli(26)= 8553103/6
B(27)= 0 sympy.bernoulli(27)= 0
B(28)= -23749461029/870 sympy.bernoulli(28)= -23749461029/870
B(29)= 0 sympy.bernoulli(29)= 0
B(30)= 8615841276005/14322 sympy.bernoulli(30)= 8615841276005/14322
B(31)= 0 sympy.bernoulli(31)= 0
B(32)= -7709321041217/510 sympy.bernoulli(32)= -7709321041217/510
B(33)= 0 sympy.bernoulli(33)= 0
B(34)= 2577687858367/6 sympy.bernoulli(34)= 2577687858367/6
B(35)= 0 sympy.bernoulli(35)= 0
B(36)= -26315271553053477373/1919190 sympy.bernoulli(36)= -26315271553053477373/1919190
B(37)= 0 sympy.bernoulli(37)= 0
B(38)= 2929993913841559/6 sympy.bernoulli(38)= 2929993913841559/6
B(39)= 0 sympy.bernoulli(39)= 0
B(40)= -261082718496449122051/13530 sympy.bernoulli(40)= -261082718496449122051/13530
B(41)= 0 sympy.bernoulli(41)= 0
B(42)= 1520097643918070802691/1806 sympy.bernoulli(42)= 1520097643918070802691/1806
B(43)= 0 sympy.bernoulli(43)= 0
B(44)= -27833269579301024235023/690 sympy.bernoulli(44)= -27833269579301024235023/690
B(45)= 0 sympy.bernoulli(45)= 0
B(46)= 596451111593912163277961/282 sympy.bernoulli(46)= 596451111593912163277961/282
B(47)= 0 sympy.bernoulli(47)= 0
B(48)= -5609403368997817686249127547/46410 sympy.bernoulli(48)= -5609403368997817686249127547/46410
B(49)= 0 sympy.bernoulli(49)= 0
B(50)= 495057205241079648212477525/66 sympy.bernoulli(50)= 495057205241079648212477525/66