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


2016年 11月 18日

せっかく早く収束する円周率の公式を使って計算をはじめたのだが、永久ループに入ってしまったようだ。
永久ループから抜ける判断をしているのは、次の個所だけである。

        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違うのは許容範囲を超える原因がありそうだ。
考えてみよう。