副題は「プログラミング問題集」。共立出版の月刊誌「bit」で毎月出題されていた。 その抜粋である。
プログラムを JavaScript で作ろうとしたができなかった。残念だが数式だけで説明する。 出題は「カタラン数」。カタランの定数 (Catalan's constant) C は次のように級数の和で定義される。 この値を十進法で小数 10 桁以上正しく求めることである。
`C=sum_(k=0)^oo (-1)^k/(2k+1)^2 = 1 - 1 / 3^2 + 1 / 5^2 - 1/ 7^2 + cdots`
答は C = 0.91596 55941 772192 … である。
電卓＋手計算で解答者が提案しているのは、既知の公式
`pi^2/48 = sum_(k=1)^oo (-1)^(k-1) / (2k)^2`
とカタランの定数の定義式を組み合わせて得られる式
`C = pi^2/48 + sum_(k=1)^oo (-1)^(k-1) (1/(2k-1)^2 + 1/ (2k)^2)`
を使うものである。上の式の第２項の級数は、オイラー変換を利用して解く、というものである。
以下は内容というより編集のことになる。
この問題と解答が第４部 組合せにあるが、出題者＝解答者である一松信が、これは，プログラミングよりも数値計算の問題ですから
と述べている通り、第６部 数値計算にもっていくべきだったろう。
カタランの名前は組合せ数学においてカタラン数の名前で知られているが、カタラン数とカタランの定数との間には直接の関係はない。 カタランの定数はある種の組合せ関数の漸近的評価に使われる。 一方、カタラン数は別のある種の `n` 種類の組み合わせを表す自然数として登場し、通常 n 番目のカタラン数 Cnとして使われる。 Cn は次の式で表される。
`C_n = ((2n)!) / ((n+1)!n!)`
こちらは perl でプログラミングしようとしてできなかった問題である。
自然数の逆数の和
`sum_(k=1)^oo 1/k ` は発散する。言い換えれば、任意の正整数 `N` に対して、 `sum_(k=1)^M 1/k ge N` となるような正整数 `M` が存在する。 そこで関数 `f(N)` を、
`f(N) = min_M (sum_(k=1)^M 1/k ge N)`
のように定義する。たとえば、`f(1) = 1, f(2) = 4, f(3) = 11, cdots` である。
12までの正整数 `N` に対して，`f(N)` を求めよ。できれば，さらに大きな整数 `N` に対しても，`f(N)` を求めよ。
素数の逆数の和
素数を小さいほうから順に`p_1, p_2, p_3, cdots` とする。`p_1=2, p_2 =3, p_3 =5, cdots` である。 無限級数`sum_(i=1)^oo 1 / p_i `も，発散することが知られている。そこで，関数 `g(N)` を
`g(N) = min_M (sum_(k=1)^M 1/p_i ge N)`
のように定義する。
`g(2)`と`g(3)`を求めよ。できれば，さらに大きな整数 `N` に対しても，`g(N)` を求めよ。注意1. 計算の精度に注意すること.
注意2. 計算結果だけでなく，計算方法(プログラムでもよい)とその計算方法で正しい答が得られる根拠(証明，あるいは説明)を明示すること.
本問は、分数の和の正確な計算をいかに行うかが問題である。この計算は誤差評価と密接に絡んでくる。 perl のプログラムについては、数列の和の計算というページに書いた。 要は、うまくいかなかった。
うまくいかなかった理由は、誤差である。結果として得られた答が正しくとも、
わずかな誤差が累積すると、単純に計算機で累積した和による `f(N)` では誤る可能性がある。
つまり、`f(13)=248397` という結果は一致しているが、
出題者のいうその計算方法で正しい答が得られる根拠(証明，あるいは説明)を明示すること.
には答えているとは言い難い。
出題者が紹介した解答では、自然数の逆数の和に関する次の公式
`sum_(k=1)^M 1/k = gamma + log M + 1/(2M) - 1/(12M(M+1)) + 1/(12M(M+1)(M+2)) - cdots`
を利用している。ここで `gamma` はオイラーの定数であり、` 0.5772156649 cdots` という数値である。
なお、上記の公式は通常次の漸近展開で表される。
`sum_(k=1)^M 1/k ~= gamma + log M + 1/(2M) - sum_(k=1)^oo (B_(2k))/(2kM^(2k))`
ここで、`B_k` はベルヌーイ数である。フルヴィッツのゼータ関数に関する結果により、次のように誤差項が評価できる：
`sum_(k=1)^M 1/k ~= gamma + log M + 1/(2M) - 1/(12M^2) + 1 / (120M^4) - epsilon`
ここで `0 lt epsilon lt 1 / (252M^6)`
左辺はほぼ `N` であるから、`1//2M` 以下の項が 1 に比べて十分に小さいとすれば
`N ~= gamma + log M`
が成り立つ。すなわち、
`N - gamma ~= log M`
`M ~= exp(N-gamma)`
ナノピコ教室では上記の `M` の近似式を受けて
`:. f(N+1) ~~ e * f(N)`
とあるが、なぜこうなるのか、理解できていない。困った。具体例で考えよう。
`f(N)` の正確な定義は上記のとおりであるが、近似ということであれば次の式を書き下していいだろう。
`f(1) = 1 ↔ 1 ~~ 1`
`f(2) = 4 ↔ 2 ~~ 1 + 1/2 + 1/3 + 1/4`
`f(3) = 11 ↔ 3 ~~ 1 + 1/2 + 1/3 + 1/4 + 1/5 + 1/6 + 1/7 + 1/8 + 1/9 + 1/10 + 1/11`
`f(2)` は `f(1)` の４倍である。これはちょっと違うが、`f(3)` は `f(2)` の 2.75 倍である。 `e = 2.718...` だから、いい線をいっている。
一方、M の近似式であるが、次々にあてはめてみよう。
`1 ~~ exp(1-gamma)`
`4 ~~ exp(2-gamma)`
`11 ~~ exp(3-gamma)`
ふむふむ。ということは、`M ~~ exp(N-gamma)` と `M' ~~ exp(N+1-gamma)` を見比べれば、 `M' ~~ M * e` がわかる。そして考えてみれば、`f(N) = M` なのである。ということで、 `f(N+1) ~~ e * f(N)` がわかった。次に、ナノピコ教室では新たな式が出てきている。
`:.f(13)~~e*f(12)~~248396`
これはどういうことだろうか。左の `~~` はわかるが、右の `~~` はどのように求めたのだろうか。
`f(12)` を実際に数列の和を求めた結果として使ったのだろうか。それとも、
`f(13) ~~ e*f(12) ~~ e^2 * f(11) ~~ cdots ~~ e^10 * f(3)`
のように求めたものだろうか。後者の方法で JavaScript を使って計算してみると、f(3) = 11 だから、
これでは足りない。おそらく、f(N) を作っていったときに、
`f(N+1) ~~ ceil(e * f(N))`
ということを N = 1, 2, 3, 4, ... , と順に確認しなければならないのだろう。地道に計算してみた。
f(1)=1, f(2)=4 は目の子でわかる。f(3)=11 からは地道に計算する必要があるだろう。
まず、`f(3) ~~ e * f(2) ~~ 2.7128 * 4 ~~ 10.8512`
よって f(3) は 10 か 11 だろう。計算してあたりをつける。 以下、`H~(M) = gamma + log M + 1/(2M) ~~ sum_(k=1)^M 1/k ` とする。
本来はベルヌーイ数が出てくるところの和の評価を行わないといけないが、省略する。 ということで、f(3) = 11である。
次は f(4) だ。`f(4) ~= e * 11 ~= ` だから、M = 29 か M = 30 を調べればいいだろう。
あれ、M = 30 でも足りない。M = 31 ではどうだろうか。
そんなことを調べていると、整数列を列挙しているページ https://oeis.org を発見した。 数列 1, 4, 11 で検索すると、https://oeis.org/A004080 がこの数列であることが確認できた。 なお、この A004080 では、f(0) = 0 としている。さて、このページを見てみると、f(4) = 31 となっている。 ただ、ベルヌーイ数の箇所を切り捨てているこということは過大評価になっているということだ。 ということは今までの結果が無になってしまう。こわいことだ。 そこで誤差項を正確に評価した式を H(M) として、改めて計算し直した。
よかった。誤差項も評価した結果、正しい結果が得られている。
さて、現在では計算機が発達しているから、手計算で一つずつ f(1), f(2), f(3), ... を求めていく必要はないと思う。
数式の表示は、MathJaxを用いている。 以前は MathML を用いていた。
| 書名 | ナノピコ教室 |
| 編者 | 駒木 悠二、有澤 誠 |
| 発行日 | |
| 発行元 | 共立出版 |
| 定 価 | 円（本体） |
| サイズ | ??版 |
| ISBN | ??? |
| その他 | 現在手元にはない |
まりんきょ学問所 ＞ コンピュータの部屋 ＞ プログラミング ＞ 駒木 悠二、有澤 誠（編）：ナノピコ教室
