在Python中计算e
怎么表示下面这个方程式最简单呢?
我想说明一下,我的问题是想要一些代码来计算这个方程的答案。
这里有两个问题:
这个求和会导致一个无限循环,根本无法得到答案。
我希望能得到一个长一点、详细一点的答案(大概40位数字左右)。
4 个回答
为了得到和我之前关于近似 pi
的回答类似的结果,你可以参考以下内容:
from functools import wraps
def memoize(f):
"""Store and retrive previous results of the decorated function f."""
cache = {}
@wraps(f)
def func(*args):
if args not in cache:
cache[args] = f(*args)
return cache[args]
return func
@memoize
def fact(n):
"""Recursively calculates n!."""
if n <= 1:
return 1
return n * fact(n - 1)
def inverse_fact_n(start_n=0):
"""Generator producing the infinite series 1/n!."""
numerator = 1.0
denominator = start_n
while True:
yield numerator / fact(denominator)
denominator += 1
def approximate_e(steps=None, tolerance=None):
"""Calculate an approximation of e from summation of 1/n!."""
if steps is None and tolerance is None:
raise ValueError("Must supply one of steps or tolerance.")
series = inverse_fact_n()
if steps is not None: # stepwise method
return sum(next(series) for _ in range(steps))
output = 0 # tolerance method
term = next(series)
while abs(term) > tolerance:
output += term
term = next(series)
return output
if __name__ == "__main__":
from math import e
print("math.e:\t\t{0:.20f}.".format(e))
stepwise = approximate_e(steps=100)
print("Stepwise:\t{0:.20f}.".format(stepwise))
tolerated = approximate_e(tolerance=0.0000000001)
print("Tolerated:\t{0:.20f}.".format(tolerated))
这个 approximate_e
函数让你可以选择:
- 一个步骤数(“我希望它花这么长时间”);或者
- 一个期望的精度(“我希望它这么准确”)。
这里面有一些相对高级的 Python 内容(比如记忆化装饰器函数和生成器函数来生成序列),但你可以只关注主要的函数,在这里 next(series)
会给你求和的下一个项。
这样我就得到了输出:
math.e: 2.71828182845904509080.
Stepwise: 2.71828182845904553488.
Tolerated: 2.71828182844675936281.
最有效的方法是利用指数函数的特性。
exp(x)=(exp(x/N))^N
这样,你可以计算 x=exp(2^(-n)),这个计算的精度比最终结果需要的多出 2n 位,然后通过将结果平方 n 次来计算 e。
对于小数字 x,当你在计算 exp(x) 的时候,如果把级数截断在 m-1 次方的项,误差会小于下一个 m 次方项的两倍。
总结一下,要以 d 位的精度计算 e,你需要选择一个适中的 n,并选择 m,使得
2^(1-mn)/m! is smaller than 2^(-d-2n)
这个 m 的确定也可以动态进行(可以参考用户22698的回答,使用 Decimal)。
from decimal import Decimal, getcontext
def eulernumber(d):
dd=d
n=4
while dd > 1:
dd /= 8
n += 1
getcontext().prec = d+n
x = Decimal(1)/Decimal(1 << n)
eps = Decimal(1)/Decimal(1 << (1 + (10*d)/3 ))
term = x
expsum = Decimal(1) + x
m = 2
while term > eps:
term *= x / Decimal(m)
m += 1
expsum += term
for k in range(n):
expsum *= expsum
getcontext().prec = d
expsum += Decimal(0)
return expsum
if __name__ == "__main__":
for k in range(1,6):
print(k,eulernumber(4*k))
for k in range(10,13):
print(k,eulernumber(4*k))
输出结果为
( 1, Decimal('2.718'))
( 2, Decimal('2.7182818'))
( 3, Decimal('2.71828182846'))
( 4, Decimal('2.718281828459045'))
( 5, Decimal('2.7182818284590452354'))
(10, Decimal('2.718281828459045235360287471352662497757'))
(11, Decimal('2.7182818284590452353602874713526624977572471'))
(12, Decimal('2.71828182845904523536028747135266249775724709370'))
想要更专业的实现这个想法,可以查看 (unix/posix) 的 bc 数学库,它也包含了对数和三角函数的实现。指数函数的代码在手册页中也有作为示例提供。
如果你需要更精确的计算,可以试试使用 Fraction
:
from fractions import Fraction # use rational numbers, they are more precise than floats
e = Fraction(0)
f = Fraction(1)
n = Fraction(1)
while True:
d = Fraction(1) / f # this ...
if d < Fraction(1, 10**40): # don't continue if the advancement is too small
break
e += d # ... and this are the formula you wrote for "e"
f *= n # calculate factorial incrementally, faster than calling "factorial()" all the time
n += Fraction(1) # we will use this for calculating the next factorial
print(float(e))
或者使用 Decimal
:
from decimal import Decimal, getcontext
getcontext().prec = 40 # set the precision to 40 places
e = Decimal(0)
f = Decimal(1)
n = Decimal(1)
while True:
olde = e
e += Decimal(1) / f
if e == olde: # if there was no change in the 40 places, stop.
break
f *= n
n += Decimal(1)
print(float(e))
下面是 e
的1000位数字:
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702761838606261331384583000752044933826560297606737113200709328709127443747047230696977209310141692836819025515108657463772111252389784425056953696770785449969967946864454905987931636889230098793127736178215424999229576351482208269895193668033182528869398496465105820939239829488793320362509443117301238197068416140397019837679320683282376464804295311802328782509819455815301756717361332069811250996181881593041690351598888519345807273866738589422879228499892086805825749279610484198444363463244968487560233624827041978623209002160990235304369941849146314093431738143640546253152096183690888707016768396424378140592714563549061303107208510383750510115747704171898610687396965521267154688957035044
为了更清楚地展示它的功能,这里有一个简化版本:
e = f = 1.0
for i in range(2, 16):
e += 1.0 / f
f *= i
print(e)
显而易见的解决方案是
import math
def e(n=10):
return sum(1 / float(math.factorial(i)) for i in range(n))
但是在n=20的时候,它的精确度会下降(和math.e相比,误差大约在10^-16左右)
要达到40位数字的精确度可能会很困难,可能需要用到任意精度的数学运算
我其实不太明白为什么需要这么精确的“e”值,因为你根本无法用这么高的精度进行任何计算(如果你对它做任何操作,你都会失去这种精度,除非你用一些任意精度的数学运算来处理)。