在Python中计算e

5 投票
4 回答
25211 浏览
提问于 2025-04-18 05:13

怎么表示下面这个方程式最简单呢?

我想说明一下,我的问题是想要一些代码来计算这个方程的答案。

这里有两个问题:

  1. 这个求和会导致一个无限循环,根本无法得到答案。

  2. 我希望能得到一个长一点、详细一点的答案(大概40位数字左右)。

4 个回答

0

为了得到和我之前关于近似 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 函数让你可以选择:

  1. 一个步骤数(“我希望它花这么长时间”);或者
  2. 一个期望的精度(“我希望它这么准确”)。

这里面有一些相对高级的 Python 内容(比如记忆化装饰器函数和生成器函数来生成序列),但你可以只关注主要的函数,在这里 next(series) 会给你求和的下一个项。

这样我就得到了输出:

math.e:     2.71828182845904509080.
Stepwise:   2.71828182845904553488.
Tolerated:  2.71828182844675936281.
0

最有效的方法是利用指数函数的特性。

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 数学库,它也包含了对数和三角函数的实现。指数函数的代码在手册页中也有作为示例提供。

8

如果你需要更精确的计算,可以试试使用 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)
6

显而易见的解决方案是

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”值,因为你根本无法用这么高的精度进行任何计算(如果你对它做任何操作,你都会失去这种精度,除非你用一些任意精度的数学运算来处理)。

撰写回答