「なかなか収束しない級数の数値計算例」を自分もやってみた

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言語)を得るには、別の考え方をするしかない。元記事はどうやってあれだけの桁数を計算しているのか、自分には見当もつかない。教えてもらいたいものだが。


参考:
library bigdecimal (Ruby 2.2.0)