なくるなつき備忘録

Pythonで再起を使って円周率を求める

Python

Pythonで再起関数の練習がてら円周率を求めるプログラムを書いてみたのでその記録です。一応、1000桁目までを求めることを目標にしていました。

このプログラムを作成するにあたり躓いた点は次の２点です。

  • 再起回数に上限がある
  • float型で扱える桁数に上限がある

この２点さえ解決すれば、あとは数学なので気楽にやっていきましょう。

再起回数の上限を変更

まずPythonはデフォルトで再起回数が1000回までに制限されているのでその上限を変更します。変更方法は簡単で標準モジュールであるsysモジュールをimportし、setrecursionlimit関数の引数に希望の値をint型で指定するだけです。

import sys
sys.setrecursionlimit(1000000)

私のPCではこんなに再起したら普通にスタックオーバーフローするはずなのでこれ以上に設定する必要はないでしょう。たぶん。心配な方は安心できるまで大きな値にしてもらって大丈夫です。

大きな桁数の計算を行えるようにする

次にfloat型(浮動小数点型)では扱える桁数に上限がありますのでこれに対応する必要があります。

float型の扱える桁数に上限があるとはどういうことでしょうか？
例えば以下のようなコードを実行してみます。

print(1e-322)
print(1e-323)
print(1e-324)
print(1e-325)
print(1e-326)

このプログラムの結果はPython 3.11時点で以下のようになります。

1e-322
1e-323
0.0
0.0
0.0

本来は1e-324などのように出力されてほしいのですが0.0になってしまっていますね。このようにあまりに小さい桁の値はfloat型の仕様上、無視されてしまうのです。

またそのほかにもfloat型には問題があり、

print(1.1 + 2.2)

の結果が3.3にならないことなんかは有名な話だと思います。

そのため数値の計算には別の型を用いました。
それがdecimalモジュールにて定義されているDecimal型です。

decimalモジュールは指定した桁数の中で正確な10進法による計算を行うことができる機能を持っているため少なくとも1000桁目までは正確に計算がしたい本記事のような場合にはぴったりです。

早速decimalモジュールをプログラムの最初でimportします。decimalモジュールもsysモジュール同様、標準モジュールなのでPythonをインストールした際に勝手に使えるようになっているはずです。

import decimal

次にdemicalモジュールで計算を行う際の精度がデフォルトでは28桁なので、環境設定を変更し、とりあえず2000桁の精度で計算を行うようにします。

decimalモジュールの計算精度は演算コンテキストにより設定されているのでdecimal.setcontext()を用いて変更していきます。

setcontext関数とはdecimal.Context型の演算コンテキストを引数に取り演算の環境設定を変更するものです。

演算コンテキストにいくつかある変数の中にprecという名前で計算の精度(precision)を決めるものがありますので、decimal.Contextクラスを生成する際に第１引数(もしくはprecとして名前付き引数)にint型の整数を指定します。

これにより指定した精度が設定された演算コンテキストが生成されます。

newContext = decimal.Context(2000)

#以下のようにしても大丈夫
newContext = decimal.Context(prec=2000)

計算精度が2000桁の演算コンテキストを生成できたので次はこの演算コンテキストをsetContext関数で適用させます。

decimal.setcontext(newContext)

一旦ここまでのプログラムをまとめてみます。

import decimal
import sys

sys.setrecursionlimit(1000000)
#再起回数の上限を変更
newContext = decimal.Context(2000)
#計算精度が2000桁のコンテキストを生成
decimal.setcontext(newContext)
#生成したコンテキストを適用

実際に円周率を求めてみる

ここまでできたらいよいよ本題の円周率を求めるプログラムを書いていきます。

今回は再起関数の練習なので逆正接関数の連分数展開を用いて円周率を求めていきました。

(詳しい導出方法などは円周率の連分数展開についてという資料に記載がありましたので興味のある方は参照してください。)

arctan(x)=x1+(1x)23+(2x)25+(3x)27+⋯arctan(x)=\cfrac{x}{1+\cfrac{(1x)^2}{3+\cfrac{(2x)^2}{5+\cfrac{(3x)^2}{7+\cdots}}}}

逆正接関数を連分数展開すると上のようになるのでx=1を代入すると

arctan(1)=π4=11+123+225+327+⋯arctan(1)=\frac{\pi}{4}=\cfrac{1}{1+\cfrac{1^2}{3+\cfrac{2^2}{5+\cfrac{3^2}{7+\cdots}}}}

このようになります。再起によりn回分数を降りるようにすると、このプログラムにより求められる円周率は次のように表されます。

π∼4×(11++123++225++327++⋯+n2(2n+1))\pi \sim 4 \times (\frac{1}{1+}+\frac{1^2}{3+}+\frac{2^2}{5+}+\frac{3^2}{7+}+\cdots+\frac{n^2}{(2n+1)})

何回の再起で1000桁まで求められるのか楽しみですね。

では実際にプログラムを書いてみます。とりあえずnは100くらいでいいでしょうか。

n=100

def calcpi(x):
    if x==0:
        return 1 / ( 1 + calcpi(x+1) )
    if x==n:
        return x**2 / ( 2*x + 1 )
    return  x**2 / ( 2*x + 1 + calcpi(x+1) )

print(calcpi(decimal.Decimal(0)) * 4)

果たして結果は…

3.141592653589793238462643383279502884197169399375105820974944592307816406286233829351709883016261375854029526481260733370440096368269497777002332851415819979261449104957113743197403338602767594670169430151859820167008082781953685687043154990277804384442670798939091628015302821937369049481762577894439496072159337856885674788507872434666769067711359016597999386988790362006810651046880689066634903236667911152996847425220525542157005152701779888293197763231363356123346485213278883090128533282247267846852458391056279202497857096331228754032463123472258579496318749402806242711342573324831804072615896050863359563228564657692058885909850006071055508716711269357945940925348984536638953440070463490197992567014674924373103517599798324519435552571750528184270390526342138291745026610107950864197022392850570496314322240272692801613870354373365953874497329709937389973267053399347877983314491297436458829543239205694779422158173594232557087798385963286295129174415699850199338156807526597744493651111579...

小数第76位まで正確でした。小数第1000位には程遠いです。

その後どんどんnの値を増やして実験してみた結果、n=1307で小数第1000位まで正しい円周率になりました。意外と早くて助かりました。

3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989...

またその後nの値をどんどん増やしてみましたが私のPCではn=7408が限界で小数第5670位までしか正確に求めることができませんでした。

どうやらこの式では約1300回計算毎に1000桁正確に求められるようです。

最後に完成形のコードを載せて終わりにしたいと思います。

import decimal
import sys

sys.setrecursionlimit(1000000)
#再起回数の上限を変更
newContext = decimal.Context(2000)
#計算精度が2000桁のコンテキストを生成
decimal.setcontext(newContext)
#生成したコンテキストを適用

n=1307

def calcpi(x):
    if x==0:
        return 1 / ( 1 + calcpi(x+1) )
    if x==n:
        return x**2 / ( 2*x + 1 )
    return  x**2 / ( 2*x + 1 + calcpi(x+1) )

print(calcpi(decimal.Decimal(0))*4)

おまけ

1行で書くとしたらこうでしょうか。

from decimal import *;import sys;sys.setrecursionlimit(9999);getcontext().prec*=99;n=1308;c=lambda x:4/(1+c(x+1))if x==0 else x**2/(2*x+1)if x==n else x**2/(2*x+1+c(x+1));print(c(Decimal(0)))