コラム 人と星とともにある数学 数学

Pythonで関・ベルヌーイ数|Pythonで数学を学ぼう! 第16回

はじめに|レポート

私が毎月開催している桜井進の算数・数学教室(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)の関・ベルヌーイ数はヤコブ・ベルヌーイの定義によるものだとわかります。

そもそも関・ベルヌーイ数を計算する簡単な公式はありません。これまでに関・ベルヌーイ数を算出するアルゴリズムはいくつも考案されています。

興味を持った方はそれらのコーディングにチャレンジしてみてはいかがでしょうか。

  • この記事を書いた人
  • 最新記事

桜井進(さくらいすすむ)様

1968年山形県生まれ。 サイエンスナビゲーター®。株式会社sakurAi Science Factory 代表取締役CEO。 (略歴) 東京工業大学理学部数学科卒、同大学大学院院社会理工学研究科博士課程中退。 東京理科大学大学院非常勤講師。 理数教育研究所Rimse「算数・数学の自由研究」中央審査委員。 高校数学教科書「数学活用」(啓林館)著者。 公益財団法人 中央教育研究所 理事。 国土地理院研究評価委員会委員。 2000年にサイエンスナビゲーターを名乗り、数学の驚きと感動を伝える講演活動をスタート。東京工業大学世界文明センターフェローを経て現在に至る。 子どもから大人までを対象とした講演会は年間70回以上。 全国で反響を呼び、テレビ・新聞・雑誌など様々なメディアに出演。 著書に『感動する!数学』『わくわく数の世界の大冒険』『面白くて眠れなくなる数学』など50冊以上。 サイエンスナビゲーターは株式会社sakurAi Science Factoryの登録商標です。

-コラム, 人と星とともにある数学, 数学
-

© 2022 株式会社インフォマティクス