どうも、木村(@kimu3_slime)です。
ネイピア数(オイラー数)\(e = 2.718\dots\)は、円周率\(\pi\)に並び、数学において重要な定数です。
指数関数\({e^x}\)は微分しても形が変わらない、自然対数\(\log_e x\)の底に利用されるだけでなく、複素関数に関するオイラーの公式、統計学における正規分布(ガウス分布)など、幅広く使われています。
今回は、\(e\)をPythonによるプログラムで、100桁以上計算する方法を紹介します。
(数学定数に関してはmath,numpyモジュールがありますが、今回はそれらを用いず1から計算します)
floatによる素直な計算
\(e\)には複数の定義があります。高校数学で扱うものは
\[ \begin{aligned}e= \lim _{n\to \infty} (1+ \frac{1}{n})^n\end{aligned} \]
です。これでも計算はできますが、収束が遅く、計算する桁\(n\)が大きすぎる割に精度ある値が求まりません。
収束がより早く計算しやすいのが、\(e^x\)の冪級数展開(テイラー展開)です。
\[ \begin{aligned}e^x = \sum _{n=1} ^\infty \frac{x^n}{n!}\end{aligned} \]
参考:なぜテイラー展開を学ぶ? 単振り子を例にわかりやすく解説、テイラー展開の展開式の覚え方、導き方、証明
これに\(x=1\)を代入すれば、\(e\)が求まります。
素直にPythonのプログラムを書いてみましょう。
1 2 3 4 5 6 7 8 9 10 11 12 | def factorial(n): # 階乗関数の定義 f=1 for i in range(n): f = f*(i+1) return f def series1(n,x): # ネイピア数の計算 series = 1 + x for i in range(2,n+1): series = series + (x**i)/factorial(i) # 1項ずつ加えていく print(series) return series |
さて、ここで「print(series1(10,1)),print(series1(100,1))」を実行します。
1 2 | 2.7182818011463845 2.7182818284590455 |
「10000 digits of e」と照らし合わせると、前者は小数点以下7桁、後者は15桁一致しています。
しかし、もうちょっと桁を求めたいものです。Pythonの小数を含む数値計算(浮動小数点演算 float)は、デフォルトで17桁までの数値しか扱えません。
参考:15. 浮動小数点演算、その問題と制限 – Pythonドキュメント
decimalで桁を増やして計算
より多くの桁を計算するために役立つのが、計算精度を任意に指定できるdecimalモジュールです。
例えば、「decimal.Decimal(1.000)」とすれば0が省略されることなく、有効数字を保ったまま数値を扱えます。
参考:decimal — 十進固定及び浮動小数点数の算術演算 – Pythonドキュメント
これを使って\(e\)を計算してみましょう。
1 2 3 4 5 6 7 | import decimal decimal.getcontext().prec = 101 # 扱える桁数を指定 # 関数定義式は上で紹介したものと同じ、省略 print(series1(100,decimal.Decimal(1))) |
1 | 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274 |
「10000 digits of e」と照らし合わせても表示されている桁はすべて一致し、小数点以下100桁求められたことになります。
僕が最初に計算しようとして、間違えたことがあります。関数定義式の方にdecimal.Decimalを使い、「print(series1(100,1))」を実行してしまったことでした。
これだと結果が20桁くらいまでは正確に求められるのですが、それ以降はいくら\(n\)を増やしても精度が増さないのです。問題は、関数に「1」を代入していることです。これはfloatの1として扱われるので、計算のうちに誤差が生じてしまうのです。decimalの1を代入すべきだった、というわけですね。
ちなみに、級数ではなく極限の定義でネイピア数を計算すると
1 2 3 4 5 6 7 | import decimal decimal.getcontext().prec = 101 def napier(n): return (1 + decimal.Decimal(1/n))**n print(napier(decimal.Decimal(10**20))) |
1 | 2.7182818284590452353466960622103672715805702442603339687181342161830472120109028383467661917294983144 |
\(n= 10 ^{20}\)であっても、小数20桁ぐらいまでしか一致しません。級数の方が小さい値(100)の代入で精度良い値が出てきますね。
ちなみに、当たり前ですが、級数による定義だから収束が早いというわけではありません。円周率を求めるなら、ライプニッツの公式は級数として計算できますが、収束が遅いです。ガウス=ルジャンドルのアルゴリズムと呼ばれる方法が、単純かつ早いようですね。
僕は最近Pythonを勉強し始めたばかりなのですが、16桁くらいしかないfloatの計算だけでは、ネイピア数の100桁も求められないじゃないかと不満がありました。今回、decimalの扱いを知り、精度の良い計算が素早くできたので満足です。
計算方法によって、真の値への収束の速さ、結果の精度が大きく変わってくることは、理論数学と違い、コンピュータならではのことして意識しなければならないなと思いました。
木村すらいむ(@kimu3_slime)でした。ではでは。
オライリージャパン
売り上げランキング: 62,016
こちらもおすすめ
Pythonで統計量関数(平均、中央値、分散、相関係数)を作り、可視化しよう
論理に関するド・モルガンの法則を真偽値の計算(プログラミング)で確かめる