python - 计算具有巨大中间值的总和

Question

我想计算

对于n高达1000000尽可能准确的值。这是一些示例代码。

from __future__ import division
from scipy.misc import comb
def M(n):
    return sum(comb(n,k,exact=True)*(1/n)*(1-k/n)**(2*n-k)*(k/n)**(k-1) for k in xrange(1,n+1))
for i in xrange(1,1000000,100):
    print i,M(i)

第一个问题是OverflowError: long int too large to convert当n = 1101. 这是因为comb(n,k,exact=True)太大而无法转换为浮点数。然而，最终结果始终是一个数字0.159。

我在如何计算具有大中间值的总和中提出了一个相关问题，但是由于三个主要原因，这个问题有所不同。

1. 我要计算的公式不同，这会导致不同的问题。
2. 从我给出的示例中可以看出，之前提出的使用 exact=True 的解决方案在这里没有帮助。编写我自己的梳子实现也行不通，因为我仍然需要执行浮点除法。
3. 我需要计算比以前更大的值的答案，这会导致新的问题。我怀疑如果不以某种巧妙的方式对总和进行编码，就无法做到这一点。

一个不会崩溃的解决方案是使用

from fractions import Fraction
def M2(n):
    return sum(comb(n,k,exact=True)*Fraction(1,n)*(1-Fraction(k,n))**(2*n-k)*Fraction(k,n)**(k-1) for k in xrange(1,n+1))
for i in xrange(1,1000000,100):
    print i, M2(i)*1.0

n=1101不幸的是，它现在太慢了，以至于我在合理的时间内 没有得到答案。

所以第二个问题是如何让它足够快地完成大型n.

Answer 1

您可以使用对数变换计算每个和，该对数变换分别用加法、减法和乘法代替乘法、除法和取幂。

def summand(n,k):
    lk=log(k)
    ln=log(n)
    a=(lk-ln)*(k-1)
    b=(log(n-k)-ln)*(2*n-k)
    c=-ln
    d=sum(log(x) for x in xrange(n-k+1,n+1))-sum(log(x) for x in xrange(1,k+1))
    return exp(a+b+c+d)

def M(n):
    return sum(summand(n,k) for k in xrange(1,n))

请注意，当 k=n 时，加数将为零，因此我不计算它，因为对数将是未定义的。

Answer 2

您可以使用gmpy2。它具有具有大指数界限的任意精度浮点运算。

from __future__ import division
from gmpy2 import comb,mpfr,fsum

def M(n):
    return fsum(comb(n,k)*(mpfr(1)/n)*(mpfr(1)-mpfr(k)/n)**(mpfr(2)*n-k)*(mpfr(k)/n)**(k-1) for k in xrange(1,n+1))
for i in xrange(1,1000000,100):
    print i,M(i)

这是输出的摘录：

2001 0.15857490038127975
2101 0.15857582611615381
2201 0.15857666768820194
2301 0.15857743607577454
2401 0.15857814042739268
2501 0.15857878842787806
2601 0.15857938657957615

免责声明：我维护 gmpy2。

Answer 3

一种相当粗暴的方法是计算所有因子，然后以这样的方式相乘以使结果保持在 1.0 左右（Python 3.x）：

def M(n):
    return sum(summand(n, k) for k in range(1, n + 1))

def f1(n, k):
    for i in range(k - 1):
        yield k
    for i in range(k):
        yield n - i

def f2(n, k):
    for i in range(k - 1):
        yield 1 / n
    for i in range(2 * n - k):
        yield 1 - k / n
    yield 1 / n
    for i in range(2, k + 1):
        yield 1 / i

def summand(n, k):
    result = 1.0
    factors1 = f1(n, k)
    factors2 = f2(n, k)
    while True:
        empty1 = False
        for factor in factors1:
            result *= factor
            if result > 1:
                break
        else:
            empty1 = True
        for factor in factors2:
            result *= factor
            if result < 1:
                break
        else:
            if empty1:
                break
    return result

因为M(1101)我得到0.15855899364641846了，但这需要几秒钟。M(2000)大约需要 14 秒并产生0.15857489065619598。

（我相信它可以优化。）