こんにちは｜こんばんは。カエルのアイコンで活動しております @kyamaz です。

はじめに

以前のエントリー『名前のついた素数コレクション〜Haskellを添えて〜』の「正則素数」を求めるプログラムで、関・ベルヌーイ数を使いました。

クンマーが、奇素数が正則であることと、素数 $p$ が $k = 2,4,6,\dots, p − 3$ における関・ベルヌーイ数の分子を割り切らないことが等価であることを示しました。その性質を利用して次のような正則素数を生成するプログラムを紹介しました。

正則素数を生成するプログラム

-- 素数列 p
p=2:3:5#p;n#x@(m:p:y)=[n|gcd m n<2]++(n+2)#last(x:[m*p:y|p^2-3<n]) 
b_n = undefined -- 関・ベルヌーイ数を求める計算（後述）
p'= [3,5,7,11]++[q| q<-drop 5 p, let bs = map (\n-> numerator . abs $ b_n!!n) [10,12..(fromEnum (q-3))], and $ map ((/=0) . (\(n :: Integer)-> mod n q)) bs]

関・ベルヌーイ数が詳しく調べられたのは、18世紀初頭に関孝和、ヤコブ・ベルヌーイにより独立で発見されました。また、諸説あるようですが、エイダ・ラブレスがチャールズ・バベッジの考案した初期の汎用計算機である解析機関についての著作の中に掲載した関・ベルヌーイ数を求めるためのプログラムのコードが、世界初のコンピュータプログラムと言われています｡¹
このように数学・計算科学とも歴史のある関・ベルヌーイ数ですが、本稿では関・ベルヌーイ数を求めるプログラムをみていきましょう。

関・ベルヌーイ数

関・ベルヌーイ数は、次の『べき乗和の公式（ファウルハーバーの公式）』の係数を表すための数として用いられました。

べき乗和の公式（ファウルハーバーの公式）
: $$\begin{flalign} &\sum_{z=1}^{n} z^{k} = \frac{1}{k+1}\sum_{j=0}^{k} {}_{k+1} C _{j} B _{j} n^{k+1-j} \tag{A}& \end{flalign}$$

関・ベルヌーイ数の定義は、次のようにされています。

関・ベルヌーイ数（定義）
: 関数 $\displaystyle \frac{x}{e^{x} -1}$ のマクローリン展開の展開計数として定義されます。
$$\begin{flalign}
& \frac{x}{e^{x} -1} = \sum_{n=0}^{\infty} B _{n} \frac{x^{n}}{n!} \tag{B}&
\end{flalign}$$

$n=1$の符号が逆転しますが、指数型母関数を $\displaystyle \frac{xe^{x}}{e^{x} -1}$ として、
$$\begin{flalign} \frac{xe^{x}}{e^{x} -1} &= \sum_{n=0}^{\infty} B _{n} \frac{x^{n}}{n!} \tag{B'} &\\
\end{flalign}$$ と表すこともあります。
そして、式(B)の両辺に$\displaystyle \frac{e^{x}-1}{x} = \sum _{n=0}^{\infty} \frac{1}{n!}x^{n-1}$をかけて、左辺を$1$にして計算すると、関・ベルヌーイ数の漸化式は次のように表せます。²

関・ベルヌーイ数（漸化式）
: $$\begin{flalign} & B_{0} = 1, B _{n} = - \frac{1}{n+1}\sum_{k=0}^{n-1} {}_{n+1} C _{k} B _{k} \tag{C}& \end{flalign}$$

関・ベルヌーイ数を求めるプログラム

関・ベルヌーイ数を求めるプログラムは、式(C)をもとにすれば計算できますが、それには組み合わせ $ {}_{n+1} C _{k} $ の値が必要となり、それなりに計算量が多くなります。組み合わせの値を使うのも工夫ができそうです³が、できれば組み合わせの値を使わずに計算できる方法を探してみましょう。
そこで、べき乗和の公式（式(A)）からはじめてみます。
式(A)の左辺を、

$$\begin{flalign}
& \sum_{z=1}^{n} z^{k} = 1^{k} + 2^{k} + \cdots + n^{k} := \mathfrak{s}_{k}(n)&
\end{flalign}$$とおいて、話しをすすめます。$k=0,1,2,3$のときを見てみると、

$$\begin{flalign}
& \mathfrak{s} _{0}(n) = 1^{0} + 2^{0} + \cdots + n^{0} = B _{0} \cdot n& \\
& \mathfrak{s} _{1}(n) = 1^{1} + 2^{1} + \cdots + n^{1} = \frac{B _{0}}{2} n^{2} + B _{1} \cdot n & \\
& \mathfrak{s} _{2}(n) = 1^{2} + 2^{2} + \cdots + n^{2} = \frac{B _{0}}{3} n^{3} + \frac{2B _{1}}{2} n^{2} + B _{2} \cdot n & \\
& \mathfrak{s} _{3}(n) = 1^{3} + 2^{3} + \cdots + n^{3} = \frac{B _{0}}{4} n^{4} + \frac{3B _{1}}{2 \cdot 3} n^{3} + \frac{3B _{2}}{2}n^{2} + B _{3} \cdot n & \\
\end{flalign}$$

関・ベルヌーイ数は、$ B _{0} = 1, B _{1} = 1/2, B _{2} = 1/6,B _{3} = 0 $ が入ります。ここで $k$ が１つ増える式を見てみると、
$$\begin{flalign}
& \mathfrak{s} _{k}(n) = \int (k \cdot \mathfrak{s} _{k-1}(n) + B _{k}) dn \tag{D}&
\end{flalign}$$ の関係があるようです。これを示すには次のようにします。
$ \mathfrak{s} _{k}(m) - \mathfrak{s} _{k}(m-1) = m^{k} $ を考えて、これを両辺 $n$ で微分すると、$ \mathfrak{s}^{\prime} _{k}(m) - \mathfrak{s}^{\prime} _{k}(m-1) = k \cdot m^{k-1} $ となります。$m=n,\dots,0$まで降順でこの関係式を使って書いてみます。

$ \mathfrak{s}^{\prime} _{k}(n) - \mathfrak{s}^{\prime} _{k}(n-1) = k \cdot n^{k-1} $
$ \mathfrak{s}^{\prime} _{k}(n-1) - \mathfrak{s}^{\prime} _{k}(n-2) = k \cdot (n-1)^{k-1} $
$:$
$:$
$ \mathfrak{s}^{\prime} _{k}(2) - \mathfrak{s}^{\prime} _{k}(1) = k \cdot 2^{k-1} $
$ \mathfrak{s}^{\prime} _{k}(1) - \mathfrak{s}^{\prime} _{k}(0) = k \cdot 1^{k-1} $

左辺と右辺をそれぞれ足し合わせると次のようになります。

$ \mathfrak{s}^{\prime} _{k}(n) - \mathfrak{s}^{\prime} _{k}(0) = k \cdot \mathfrak{s} _{k-1}(n)$

ここで、$\mathfrak{s}^{\prime} _{k}(0) = B _{k}$ であることを使うと、

$ \mathfrak{s}^{\prime} _{k}(n) = k \cdot \mathfrak{s} _{k-1}(n) + B _{k}$

が成り立ち、両辺を $n$ で積分すると式(D)が得られます。関・ベルヌーイ数の性質の１つを紹介したいために微積分に関する関係式を紹介しましたが、この延長にもうまい計算方法が探せなさそうです。

秋山・谷川三角形

秋山・谷川三角形とは、秋山茂樹、谷川好男両氏が提案した調和数列のある種の階差数列を反復計算することで関・ベルヌーイ数を効率的に計算するアルゴリズムです。

詳しくはせきゅーんさんのはてなブログに解説があります。

同じような解説をこちらに示すのも憚れる（おわりに、写経しています）ので、ここではプログラムの紹介に移りたいと思います。結果を知っての出力ですが、秋山・谷川三角形は以下のような数列の列の先頭を取ってくればよいことになります。

[1 % 1,1 % 2,1 % 3,1 % 4,1 % 5,1 % 6,1 % 7,1 % 8,1 % 9,1 % 10,1 % 11,1 % 12,1 % 13],
[1 % 2,1 % 3,1 % 4,1 % 5,1 % 6,1 % 7,1 % 8,1 % 9,1 % 10,1 % 11,1 % 12,1 % 13,1 % 14],
[1 % 6,1 % 6,3 % 20,2 % 15,5 % 42,3 % 28,7 % 72,4 % 45,9 % 110,5 % 66,11 % 156,6 % 91,13 % 210],
[0 % 1,1 % 30,1 % 20,2 % 35,5 % 84,5 % 84,7 % 120,28 % 495,3 % 55,15 % 286,55 % 1092,22 % 455,13 % 280],
[(-1) % 30,(-1) % 30,(-3) % 140,(-1) % 105,0 % 1,1 % 140,49 % 3960,8 % 495,27 % 1430,125 % 6006,121 % 5460,3 % 130,169 % 7140],
[0 % 1,(-1) % 42,(-1) % 28,(-4) % 105,(-1) % 28,(-29) % 924,(-7) % 264,(-28) % 1287,(-87) % 5005,(-27) % 2002,(-11) % 1092,(-11) % 1547,(-13) % 2856],
[1 % 42,1 % 42,1 % 140,(-1) % 105,(-5) % 231,(-9) % 308,(-343) % 10296,(-1576) % 45045,(-27) % 770,(-205) % 6006,(-605) % 18564,(-95) % 3094,(-3887) % 135660],
[0 % 1,1 % 30,1 % 20,8 % 165,5 % 132,295 % 12012,67 % 5720,4 % 6435,(-6) % 715,(-75) % 4862,(-55) % 2652,(-517) % 20995,(-7423) % 271320],
[(-1) % 30,(-1) % 30,1 % 220,7 % 165,200 % 3003,1543 % 20020,3997 % 51480,464 % 6435,1539 % 24310,775 % 14586,10769 % 251940,9643 % 293930,169 % 7140],
[0 % 1,(-5) % 66,(-5) % 44,(-44) % 455,(-629) % 12012,(-41) % 12012,133 % 3432,140 % 1989,1113 % 12155,9597 % 92378,38555 % 352716,3223 % 29393,21177 % 198968],[5 % 66,5 % 66,(-1017) % 20020,(-2663) % 15015,(-35) % 143,(-1013) % 4004,(-38759) % 175032,(-18536) % 109395,(-51219) % 461890,(-105155) % 1939938,(-1331) % 352716,3567 % 92378,2494609 % 34321980],
[0 % 1,691 % 2730,691 % 1820,368 % 1365,15 % 364,(-3515) % 18564,(-9653) % 26520,(-88508) % 188955,(-74976) % 146965,(-326115) % 646646,(-164455) % 352716,(-1381937) % 3380195,(-2132689) % 6240360],
[(-691) % 2730,(-691) % 2730,601 % 1820,1247 % 1365,5350 % 4641,32421 % 30940,1104901 % 1511640,441824 % 1322685,(-170073) % 3233230,(-738425) % 1939938,(-25619891) % 40562340,(-5441531) % 6760390,(-2836327) % 3120180]]

Haskellでは、無限数列の列でも先頭だけを取ってくる記述ができます。それが次のようなプログラムになります。

関・ベルヌーイ数を求めるプログラム

import Data.Ratio
b_n = map head . iterate (b' 1) . map (\i->1%(i+1)) $ enumFrom 0 where b' n ms = case ms of [_] -> []; (x:y:xs) -> (n%1)*(x-y) : b' (n+1) (y:xs)

おわりに

いかがでしたでしょうか。折角ですので、ここで金子昌信氏による証明⁴（せきゅーんさんのブログにも記載）を写経しておきます。（※長いので折り畳んでおきます。）

秋山・谷川三角形の解説

秋山・谷川三角形の数列から関・ベルヌーイ数を導出することを考えます。まず、秋山・谷川三角形の数列は $\{ b _{n,k} \} _{n,k \ge 0}$ とおいて、初期値は$\displaystyle \{ b _{0,k} \} = \{\frac{1}{1}, \frac{1}{2}, \frac{1}{3},\dots ,\frac{1}{k+1} \}$ とします。漸化式を以下で定義すると、$ b _{n+1,k} = (k+1) (b _{n,k} - b _{n,k+1}) $ となり、$b _{n,0}$ が関・ベルヌーイ数 $B _{n}$ であることを示します。
まず、母関数 $\displaystyle f _{n}(t) := \sum _{k=0}^{\infty} b _{n,k} t^{k}$ を考えます。
$ t = 0 $ のとき、$ f _{n}(0) = b _{n,0} $ となるので、$ t = 0 $ のときの値を計算していくことを目指します。

さて、$ n \ge 1 $ のとき、漸化式より
$$\begin{flalign}
f _{n}(t) &= \sum _{k=0}^{\infty} (k+1) (b _{n-1,k} - b _{n-1,k+1}) t^{k} &\\
&= \frac{d}{dt}\Big(\sum _{k=0}^{\infty} b _{n-1,k} t^{k+1} \Big) - \frac{d}{dt}\Big(\sum _{k=0}^{\infty} b _{n-1,k+1} t^{k+1} \Big) \\
&= \frac{d}{dt}\Big\{ t f _{n-1}(t) \Big\} - \frac{d}{dt}\Big\{ f _{n-1}(t) - b _{n-1,0} \Big\} \\
&= \frac{d}{dt}\Big\{ (t - 1) f _{n-1}(t) \Big\} \\
\end{flalign}$$
ここで、$ g _{n}(t) := (t - 1) f _{n}(t) $ とおくと、

$$\begin{flalign}
g _{n}(t) &= (t - 1) \frac{d}{dt}\Big\{ (t - 1) f _{n-1}(t) \Big\} &\\
&= (t - 1) \frac{d}{dt} g _{n-1}(t) \\
&= \Big( (t - 1) \frac{d}{dt} \Big) ^{n} g _{0}(t) \tag{★}\\
\end{flalign}$$

次に、第二種スターリング数 $\displaystyle \Big\{ \begin{matrix}
n\\ m \end{matrix} \Big\}$ の特性として下式を用います。

$$\begin{flalign}
\Big(x \frac{d}{dx} \Big)^{n} &= \sum _{m=0}^{n} \Big\{ \begin{matrix} n\\ m \end{matrix} \Big\} x^{m} \Big(\frac{d}{dx} \Big)^{m} &\\
\end{flalign}$$

この式に $ x \rightarrow (t - 1)$ とした式

$$\begin{flalign}
\Big((t-1) \frac{d}{dt} \Big)^{n} &= \sum _{m=0}^{n} \Big\{ \begin{matrix} n\\ m \end{matrix} \Big\} (t-1)^{m} \Big(\frac{d}{dt} \Big)^{m} &\\
\end{flalign}$$

を式(★)に適用して、

$$\begin{flalign}
g _{n}(t) &= \sum _{m=0}^{n} \Big\{ \begin{matrix} n\\ m \end{matrix} \Big\} (t-1)^{m} \Big(\frac{d}{dt} \Big)^{m} g _{0}(t)&\\
\end{flalign}$$

を得ます。この式の $ t = 0 $ のときで式を変形してみしましょう。

$$\begin{flalign}
{-} f _{n}(0) = - b _{n,0} &= \sum _{m=0}^{n} \Big\{ \begin{matrix} n\\ m \end{matrix} \Big\} (-1)^{m} \Big(\frac{d}{dt} \Big)^{m} g _{0}(t)&\\
&= \sum _{m=0}^{n} \Big\{ \begin{matrix} n\\ m \end{matrix} \Big\} (-1)^{m} m! (b _{0,m-1} - b _{0,m}) \\
&= \sum _{m=0}^{n-1} \Big\{ \begin{matrix} n\\ m+1 \end{matrix} \Big\} (-1)^{m+1} (m+1)! b _{0,m} - \sum _{m=0}^{n} \Big\{ \begin{matrix} n\\ m \end{matrix} \Big\} (-1)^{m} m! b _{0,m} \\
&= - \sum _{m=0}^{n} \Big( (m+1)\Big\{ \begin{matrix} n\\ m+1 \end{matrix} \Big\} + \Big\{ \begin{matrix} n\\ m \end{matrix} \Big\} \Big) (-1)^{m} m! b _{0,m}\\
&= - \sum _{m=0}^{n} \Big\{ \begin{matrix} n+1\\ m+1 \end{matrix} \Big\} (-1)^{m} m! b _{0,m} \\
&= - \sum _{m=0}^{n} \Big\{ \begin{matrix} n+1\\ m+1 \end{matrix} \Big\} (-1)^{m} \frac{m!}{m+1} \\
\end{flalign}$$

（ここでは、第二種スターリング数の特性 $\displaystyle \Big\{ \begin{matrix}
n+1\\ m+1 \end{matrix} \Big\} = (m+1)\Big\{ \begin{matrix} n\\ m+1 \end{matrix} \Big\} + \Big\{ \begin{matrix} n\\ m \end{matrix} \Big\}$ を使いました）

これが、関・ベルヌーイ数であることを示したいので、少し強引ですが式(B)（式(B')でも同様）の右辺の形を作って計算してみます。

$$\begin{flalign}
& \sum _{n=0}^{\infty}\Big( \sum _{m=0}^{n} \Big\{ \begin{matrix} n+1\\ m+1 \end{matrix} \Big\} (-1)^{m} \frac{m!}{m+1} \Big)\frac{x^{n}}{n!} & \\
&= \sum _{m=0}^{\infty} (-1)^{m} \frac{m!}{m+1} \sum _{n=m}^{\infty}\Big\{ \begin{matrix} n+1\\ m+1 \end{matrix} \Big\}\frac{x^{n}}{n!} \\
&= \sum _{m=0}^{\infty} (-1)^{m} \frac{m!}{m+1} \frac{e^{x}(e^{x}-1)^{m}}{m!} \\
&= e^{x} \sum _{m=0}^{\infty} \frac{(1- e^{x})^{m}}{m+1} = \frac{e^{x}}{1- e^{x}} \sum _{m=1}^{\infty} \frac{(1- e^{x})^{m}}{m} \\
&= \frac{e^{x}}{1- e^{x}} \big(- \log (1 - (1 - e^{x})) \big) = - \frac{x e^{x}}{1- e^{x}} \\
\end{flalign}$$

これは、式(B')の指数型母関数の形式と同じですので、目的の次式が示せたことになります。

$$\begin{flalign}
b _{n,0} &= \sum _{m=0}^{n} \Big\{ \begin{matrix} n+1\\ m+1 \end{matrix} \Big\} (-1)^{m} \frac{m!}{m+1} = B _{n} & \\
\end{flalign}$$

最後に、本稿を記載するために検証したHaskell環境を記しておきます。お手元の環境で検証する際に、動作が異なるときには参考になるかもしれません。

本稿の環境

本稿のために使用した環境は以下となります。

macOS: Sonoma 14.5 (chip: Apple M1)
GHCup: 0.1.30.0
GHC: 9.6.4

ご一読いただきまして有り難うございます。
（●）（●） Happy Hacking!
/""　＿＿""＼

@gretchiさんのエントリ『計算機オタクがコンピュータの歴史と現代CPUについて語る』より ↩
高校数学の美しい物語『ベルヌーイ数とゼータ関数』 ↩
Qiita）ベルヌーイ数を線形時間で求める ↩
Journal of Integer Sequences, Vol. 3 (2000), Article 00.2.9 ↩