ネイピア数eをPython（decimal）で100桁計算してみよう

どうも、木村（@kimu3_slime）です。

ネイピア数（オイラー数）\(e = 2.718\dots\)は、円周率\(\pi\)に並び、数学において重要な定数です。

指数関数\({e^x}\)は微分しても形が変わらない、自然対数\(\log_e x\)の底に利用されるだけでなく、複素関数に関するオイラーの公式、統計学における正規分布（ガウス分布）など、幅広く使われています。

今回は、\(e\)をPythonによるプログラムで、100桁以上計算する方法を紹介します。

（数学定数に関してはmath,numpyモジュールがありますが、今回はそれらを用いず1から計算します）

floatによる素直な計算
decimalで桁を増やして計算
こちらもおすすめ

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のプログラムを書いてみましょう。

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

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))」を実行します。

2.7182818011463845
2.7182818284590455

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\)を計算してみましょう。

import decimal

decimal.getcontext().prec = 101 # 扱える桁数を指定

# 関数定義式は上で紹介したものと同じ、省略

print(series1(100,decimal.Decimal(1)))

import decimal

decimal.getcontext().prec = 101 # 扱える桁数を指定

# 関数定義式は上で紹介したものと同じ、省略

print(series1(100,decimal.Decimal(1)))

2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274

1	2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274

「10000 digits of e」と照らし合わせても表示されている桁はすべて一致し、小数点以下100桁求められたことになります。

僕が最初に計算しようとして、間違えたことがあります。関数定義式の方にdecimal.Decimalを使い、「print(series1(100,1))」を実行してしまったことでした。

これだと結果が20桁くらいまでは正確に求められるのですが、それ以降はいくら\(n\)を増やしても精度が増さないのです。問題は、関数に「1」を代入していることです。これはfloatの1として扱われるので、計算のうちに誤差が生じてしまうのです。decimalの1を代入すべきだった、というわけですね。

ちなみに、級数ではなく極限の定義でネイピア数を計算すると

import decimal
decimal.getcontext().prec = 101

def napier(n):
    return (1 + decimal.Decimal(1/n))**n

print(napier(decimal.Decimal(10**20)))

import decimal

decimal.getcontext().prec = 101

def napier(n):

return (1 + decimal.Decimal(1/n))**n

print(napier(decimal.Decimal(10**20)))

2.7182818284590452353466960622103672715805702442603339687181342161830472120109028383467661917294983144

1	2.7182818284590452353466960622103672715805702442603339687181342161830472120109028383467661917294983144

\(n= 10 ^{20}\)であっても、小数20桁ぐらいまでしか一致しません。級数の方が小さい値（100）の代入で精度良い値が出てきますね。

ちなみに、当たり前ですが、級数による定義だから収束が早いというわけではありません。円周率を求めるなら、ライプニッツの公式は級数として計算できますが、収束が遅いです。ガウス＝ルジャンドルのアルゴリズムと呼ばれる方法が、単純かつ早いようですね。

僕は最近Pythonを勉強し始めたばかりなのですが、16桁くらいしかないfloatの計算だけでは、ネイピア数の100桁も求められないじゃないかと不満がありました。今回、decimalの扱いを知り、精度の良い計算が素早くできたので満足です。

計算方法によって、真の値への収束の速さ、結果の精度が大きく変わってくることは、理論数学と違い、コンピュータならではのことして意識しなければならないなと思いました。

木村すらいむ（@kimu3_slime）でした。ではでは。

Pythonからはじめる数学入門

posted with amazlet at 19.09.10

Amit Saha
オライリージャパン
売り上げランキング: 62,016

Amazon.co.jpで詳細を見る

こちらもおすすめ

Pythonで統計量関数（平均、中央値、分散、相関係数）を作り、可視化しよう

2次方程式をプログラムで解くときに気をつける「誤差」とは？

論理に関するド・モルガンの法則を真偽値の計算（プログラミング）で確かめる