はじめに｜レポート
私が毎月開催している桜井進の算数・数学教室（Zoomオンラインセミナー）。
11/21桜井進の魔法の算数教室、12/12桜井進のPython・UNIX・Math教室のテーマはそれぞれ「小学生からわかる超入門関・ベルヌーイ数」「超入門・Pythonで関・ベルヌーイ数」でした。
小学生1年生から大人まで多数参加いただき盛況に終わりました。
目次
関・ベルヌーイ数
ゼータ関数を語る時、避けて通れないのが関・ベルヌーイ数です。
ゼータの負の自然数に対する公式を見れば明らかですが、ゼータは関・ベルヌーイ数でできています。
高校数学に登場するべき乗和の公式の係数にも関・ベルヌーイ数は現れます。べき乗和の公式を一般化したものが関・ベルヌーイの公式に他なりません。
クラスではまず、関・ベルヌーイ数Bnを使ってみる練習を行います。
関・ベルヌーイ数自体の算出は少し高度なので後回しにします。
最初のステップとして、関・ベルヌーイ数の数表の数値を用いて、べき乗和の公式を計算してみます。
その際に二項係数（コンビネーション）も必要になるので、あわせてパスカルの三角形の解説も行います。
次のステップとして、関・ベルヌーイ数の算出に進んでいきます。
1712年『活用算法』に登場した関・ベルヌーイ数
1712年、わが国において関・ベルヌーイ数は世界で初めて登場しました。それが関孝和の死後に弟子たちによって編纂された『活用算法』の朶積術（だせきじゅつ）です。朶とは垂れているものという意味です。
ヤコブ・ベルヌーイの『推測術』が著されたのは、その1年後の1713年です。隔絶した世界の2地点で、同時期に同じことが考えられていたことに驚くばかりです。
ベルヌーイ数Bnとして世界には知られるようになりましたが、事実は関孝和の方が出版が早いので、本来は関数（せきすう）Snと呼ばれてもよいのですが、1年差を考慮して関・ベルヌーイ数と呼ぶべきでしょう。
ちなみにウィキペディアのベルヌーイ数を見てみると、拙書『天才たちが愛した美しい数式』（PHP研究所、2008）で、ベルヌーイ数は関・ベルヌーイ数と呼ぶべきであると書かれてあることが紹介されています。
関・ベルヌーイ数の漸化式
関・ベルヌーイ数は次のような式で定義されます。これはヤコブ・ベルヌーイによる定義です。
これをPythonコードに翻訳して関・ベルヌーイ数を算出することを目指します。
このシグマの式からどのようにBnが求まっていくかというと、まずn=1としてB0からB1が求まります。次にn=2としてB0、B1からB2が求まります。
次にn=3としてB0、B1、B2からB3が求まります。
漸化式をPthonコードに翻訳
シグマの式をBn=…と式変形することで、Pythonコードに翻訳することを考えます。
次の図はその詳細です。一番式の式がBn=…であることがわかります。
この式を用いることで関数B(n)を再帰的定義します。その詳細が次の図です。
コード1「seki-bernoulli_1.py」
>>> # 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))
このプログラムファイルは次からダウンロードできます
https://drive.google.com/file/d/1oI8u6PutBn5zow9g6Q0ifpf7U8X2GXl0/view?usp=sharing
実行結果
python seki-bernoulli_1.py
として実行させると、
B(n)のnの上限がたずねられます。
20としてみたのが次の結果です。
SymPyライブラリをインポートしてsympy.bernoulli(n)とすればn番目の関・ベルヌーイ数を返してくれます。
そこでB(n)に対してsympy.bernoulli(n)も表示するようにしました。
B(1)= 1/2 であるのに対してsympy.bernoulli(2)= -1/2です。
これは、ヤコブ・ベルヌーイによるオリジナルな関・ベルヌーイ数の定義がB(1)= 1/2であるのに対して、現在の関・ベルヌーイ数は指数関数を用いた関数をマクローリン展開したときの展開係数として定義されます。
それがsympy.bernoulli()です。
歴史的にこの2つの定義が混在しています。一番目の関・ベルヌーイ数だけの違いですが、関・ベルヌーイ数を扱う場合には注意が必要です。
このプログラムによる関・ベルヌーイ数の算出はn=22が限界です。再帰的定義が原因です。
このプログラムはヤコブ・ベルヌーイによる定義を忠実に翻訳してみたものです。
コード2「seki-bernoulli_2.py」
>>> # seki-bernoulli_2.py
>>>
>>> print('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))
このプログラムファイルは次からダウンロードできます
https://drive.google.com/file/d/1PhkgXMSbvc5BoMuKgLgK6FLiJNjcvRg-/view?usp=sharing
関・ベルヌーイ数とゼータ関数の関係から関・ベルヌーイ数を求めるプログラムです。
その関係式はいたってシンプルです。
Bm = - m * zeta(1-m)
ゼータ関数はmpmathライブラリにある関数zeta()を用いています。
実行結果
B(1)= 1/2であることから、関係式Bm = - m * zeta(1-m)の関・ベルヌーイ数はヤコブ・ベルヌーイの定義によるものだとわかります。
そもそも関・ベルヌーイ数を計算する簡単な公式はありません。これまでに関・ベルヌーイ数を算出するアルゴリズムはいくつも考案されています。
興味を持った方はそれらのコーディングにチャレンジしてみてはいかがでしょうか。
