SymPyで円周率がちゃんと１００万桁計算できたかな
2016年 11月 18日
せっかく早く収束する円周率の公式を使って計算をはじめたのだが、永久ループに入ってしまったようだ。
永久ループから抜ける判断をしているのは、次の個所だけである。
if abs(t)<10**(-len):
break
これから考えられるのは、この条件式がいつまで経っても成立しないことだ。
tは高精度な計算をしているので、たぶん問題ないだろう。
限界を決めるのに、１０のべき乗を利用しているのだが、もしかしてこれが問題だろうか。
それで、こんな実験をしてみた。
In [26]: 10**(-400)==0.0
Out[26]: True
In [27]: 10**(-300)==0.0
Out[27]: False10**(-400)は、0.0と等しくなってしまったようだ。
ということは、
if abs(t)<0.0:という絶対に成立しないことをしていたようだ。
高精度の計算を行うには、もっと微小な値をちゃんと計算し、比較しないといけない。
どうすればよいだろうか？
ということで、微小な下限値を、ｔと同様に、まず式で求め、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桁小さくした。
これで、１００１桁で実行してみた。
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時間。
一番最後の桁が３だけ違った。
まあいいか。
でもないかな。
１違うのは仕方がないにしても、３違うのは許容範囲を超える原因がありそうだ。
考えてみよう。
＃人気の記事
＃タグ
