「なかなか収束しない級数の数値計算例」を自分もやってみた
id:Hyperion64 さんのブログの、
なかなか収束しない級数の数値計算例 - 完全無欠で荒唐無稽な夢
の記事が興味深かったので、ちょっと Ruby でやってみた。
級数
の極限を求めようというのである。こういうコードでやってみた。Ruby は、任意の精度で 10 進表現された浮動小数点数を扱うことができる。また、lambda のクロージャとしての性質を使っている。
require 'bigdecimal' def sum1 a = BigDecimal.new("0"); n = 0 lambda do n += 1 a += 1 / BigDecimal.new(n) end end def sum2 x = sum1 b = BigDecimal.new("0"); s = -1 lambda do s *= -1 b += s * (1 / x.call) end end y = sum2 999999.times do y.call end puts y.call #=>0.59159271505818740835529550854823676E0 puts y.call #=>0.661072248003938175801125710276193031E0
これで第100万項と第100万1項が出る(実行時間は 8.6 秒くらい)。
もう一桁増やしてみると、第1000万項、第1000万1項はそれぞれ
0.596383954822163458431155651252054767E0
0.656281010474278423000546156574263931E0
でなかなか収束しない。
では、さらに二桁増やしてみる。結果が出るのは明日の朝になるだろう。(AM2:35)
→古い PC をつけっぱなしにして計算させたのだけれど、何故かコマンドプロンプトが消えていた…。(AM8:49)
よく考えたら、こんな冗長なコードにする必要はなかった。
require 'bigdecimal' sum = BigDecimal.new("0") dnm = BigDecimal.new("0") s = 1 for i in 1..100_0000 dnm += 1 / BigDecimal.new(i) sum += s * (1 / dnm) s *= -1 end puts sum #=>0.59159271505818740835529550854823676E0 dnm += 1 / BigDecimal.new(i + 1) puts sum += s * (1 / dnm) #=>0.661072248003938175801125710276193031E0
これでいいな。
これで、10億項目、10億1項目がそれぞれ
0.602858834973980440783949444112709997E0
0.649806130500742268711055098195187070799508852E0
となる。古い PC で 12443 秒(約 3.46時間)かかった。これでは時間がかかりすぎて、収束の値などというところまでいかない。根本的に別の方法でやらないと。C言語でやるという手もあるが、C言語の浮動小数点演算は何桁まであるのか。
なお、上の二数の平均はおおよそ 0.626332482 ほどで、上のリンク先記事の収束値とだいたい同じ値が出る。
C言語でやってみる。
#include <stdio.h> int main(void) { double sum = 0, dnm = 0; int s = 1, i; for (i = 1; i <= 1000000; i++) { dnm += 1 / (double)i; sum += s * (1 / dnm); s *= -1; } printf("%17.15f\n", sum); //=>0.591592715058186 dnm += 1 / ((double)i + 1); sum += s * (1 / dnm); printf("%17.15f\n", sum); //=>0.661072248003945 return 0; }
Ruby での場合と比べてみると、既に最後の方の桁が信用できないことがわかる。
10億項目、10億1項目を計算すると、それぞれ
0.602858834975832
0.649806130502593
をわずか 23 秒で計算するが、これも Ruby の場合と比べると、最後 4桁が信用できない。このあたりは計算上の分母がまったく不正確なので、殆ど無意味だと思う。
参考としてさらに 2桁増やしたいが、ループ変数が int型の範囲を超えてしまう。なので、Visual C++ の __int64型(非標準)を使う(参照)。1000億項目、1000億1項目はそれぞれ
0.607031675099378
0.645633290380083
と求められるが、上に出した平均値ほどの意味しかあるまい。ちなみにこれで 1085 秒(約 18 分)かかった。
いずれにせよ、これ以上のスピード(Ruby)あるいは精度(C言語)を得るには、別の考え方をするしかない。元記事はどうやってあれだけの桁数を計算しているのか、自分には見当もつかない。教えてもらいたいものだが。