TIM Labs

SymPyで円周率がちゃんと100万桁計算できたかな

| コメント(0) | トラックバック(0)
せっかく早く収束する円周率の公式を使って計算をはじめたのだが、永久ループに入ってしまったようだ。
永久ループから抜ける判断をしているのは、次の個所だけである。
        if abs(t)<10**(-len):
            break

これから考えられるのは、この条件式がいつまで経っても成立しないことだ。 tは高精度な計算をしているので、たぶん問題ないだろう。 限界を決めるのに、10のべき乗を利用しているのだが、もしかしてこれが問題だろうか。 それで、こんな実験をしてみた。
In [26]: 10**(-400)==0.0
Out[26]: True

In [27]: 10**(-300)==0.0
Out[27]: False
10**(-400)は、0.0と等しくなってしまったようだ。 ということは、
        if abs(t)<0.0:
という絶対に成立しないことをしていたようだ。 高精度の計算を行うには、もっと微小な値をちゃんと計算し、比較しないといけない。 どうすればよいだろうか?

ということで、微小な下限値を、tと同様に、まず式で求め、evalfで桁数を指定して評価するように変更してみた。
def FabriceBellard(len):
    var('a:z')
    term = (-1)**n/2**(10*n+6) * (-2**5/(4*n+1) - 1/(4*n+3) + 2**8/(10*n+1) 
        - 2**6/(10*n+3) - 2**2/(10*n+5) - 2**2/(10*n+7) + 1/(10*n+9))
    eps  = 10**(-n).evalf(2,subs={n:len+2})
    sigma = 0
    m = 0
    while True:
        t = term.evalf(len+2,subs={n:m})
        if abs(t)<eps:
            break
        sigma += t
        m += 1
    return sigma.evalf(len)
epsの値だが、余裕をもって2桁小さくした。
これで、1001桁で実行してみた。
In [51]: len = 1001

In [51]: str(FabriceBellard(len))[-50:]
Out[51]: '18577805321712268066130019278766111959092164201989'

In [52]: str(sympy.pi.evalf(len))[-50:]
Out[52]: '18577805321712268066130019278766111959092164201989'
大丈夫なようだ。
時間計測を含めて、ちょっと100万桁求めてみた。

In [68]: len = 1000001

In [69]: st = time.time()

In [70]: p = FabriceBellard(len)

In [71]: time.time()-st
Out[71]: 43364.538419008255

In [72]: str(p)[-50:]
Out[72]: '56787961303311646283996346460422090106105779458148'

In [73]: str(sympy.pi.evalf(len))[-50:]
Out[73]: '56787961303311646283996346460422090106105779458151'

43364秒かかった。約12時間。 一番最後の桁が3だけ違った。
まあいいか。
でもないかな。
1違うのは仕方がないにしても、3違うのは許容範囲を超える原因がありそうだ。

考えてみよう。

トラックバック(0)

トラックバックURL: http://labs.timedia.co.jp/mt/mt-tb.cgi/531

コメントする

このブログ記事について

このページは、fujiが2016年11月18日 00:00に書いたブログ記事です。

ひとつ前のブログ記事は「SymPyで高性能な円周率の公式を使ったら計算が止まらない」です。

次のブログ記事は「書評:『行列プログラマー』」です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。