|
|---|
|
これは日々の作業を通して学んだことや毎日の生活で気づいたことをを記録しておく備忘録である。
HTML ファイル生成日時: 2024/09/11 17:53:24.921 (台灣標準時)
インド人の研究者によって、円周率を求めることに使える新しい公式が見つけ られたそうでござる。
この公式を使って、円周率の近似値を計算するプログラムを書いてみたでござ る。
#!/usr/pkg/bin/python3.12 # # Time-stamp: <2024/08/02 13:06:52 (UT+8) daisuke> # # importing math module import math # a function to calculate approximate value of pi def new_pi (l, n): # initialising value of pi pi = 4.0 # calculating n terms for i in range (1, n+1): a = 1.0 / (i + l) - 4.0 / (2.0 * i + 1.0) b = (2.0 * i + 1.0)**2 / (4.0 * (i + l) ) - i c = i - 1.0 d = math.factorial (i) pi += a * math.gamma (b + c) / math.gamma (b) / d # returning calculated approximate value of pi return (pi) # calculation of approximate value of pi pi = new_pi (10.0, 100) # printing result print (f'approximate value of pi') print (f'50 digits: 3.14159265358979323846264338327950288419716939937510') print (f'new_pi (10,100): {pi}')
実行結果は以下の通りでござる。
% ./pi_2024.py approximate value of pi 50 digits: 3.14159265358979323846264338327950288419716939937510 new_pi (10,100): 3.1415926535897927
λ=10 で、 100 項も計算すれば、 15 桁くらいまで合っているようでござる。
|
|
|---|
|
|
上の方法だと、倍精度の範囲内でしか円周率の近似値を求められないでござる。 更に、桁数を増やすにはどうすればよいのか考えてみたでござる。 Python に 標準で含まれる Decimal モジュールを使えばよいかと思ったのでござるが、 ガンマ関数の計算をする部分を自分で書かねばならぬので、面倒そうでござっ た。調べてみると、 mpmath というパッ ケージがあることがわかったでござる。このパッケージにはガンマ関数を計算 してくれる関数も用意されているようなので、これを使えば簡便に円周率の近 似値を求められそうでござる。以下のようなコードを準備してみたでござる。
#!/usr/pkg/bin/python3.12 # # Time-stamp: <2024/08/02 14:50:52 (UT+8) daisuke> # # importing mpmath module import mpmath # settings for mpmath module mpmath.mp.dps = 100 mpmath.mp.pretty = True # a function to calculate approximate value of pi def new_pi (l, n): # initialising value of pi mp_1 = mpmath.mpf (1) mp_2 = mpmath.mpf (2) mp_4 = mpmath.mpf (4) mp_l = mpmath.mpf (l) mp_pi = mp_4 # calculating n terms for i in range (1, n+1): mp_i = mpmath.mpf (i) a = mp_1 / (mp_i + mp_l) - mp_4 / (mp_2 * mp_i + mp_1) b = (mp_2 * mp_i + mp_1)**mp_2 / (mp_4 * (mp_i + mp_l) ) - mp_i c = mp_i - mp_1 d = mpmath.factorial (i) mp_pi += a * mpmath.gamma (b + c) / mpmath.gamma (b) / d # returning calculated approximate value of pi return (mp_pi) # calculation of approximate value of pi l = 10 n = 10000000 pi = new_pi (l, n) # printing result print (f'approximate value of pi:') print (f'100 digits:') print (f'3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679') print (f'new_pi ({l},{n}):') print (f'{pi}')
実行結果は、以下の通りでござる。
% ./pi_2024_100.py approximate value of pi: 100 digits: 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679 new_pi (10,10000000): 3.141592653589793238462643383279502884197169399375105820974944592307816406286206269139035079106063948
1000 万項の計算をすると、 78 桁くらいまでいけるようでござる。
|
|---|
|