Python中的高效求和

Question

我正在尝试在Python中有效地计算求和的总和：

WolframAlpha能够计算出过高的 n 值：sum 的总和。

我有两种方法：for循环方法和np.sum方法。我认为 np.sum 方法会更快。然而，它们在 n 很大之前都是相同的，之后 np.sum 会出现溢出错误并给出错误的结果。

我正在尝试找到计算这个总和的最快方法。

import numpy as np
import time

def summation(start,end,func):
    sum=0
    for i in range(start,end+1):
        sum+=func(i)
    return sum

def x(y):
    return y

def x2(y):
    return y**2

def mysum(y):
    return x2(y)*summation(0, y, x)

n=100

# method #1
start=time.time()
summation(0,n,mysum)
print('Slow method:',time.time()-start)

# method #2
start=time.time()
w=np.arange(0,n+1)
(w**2*np.cumsum(w)).sum()
print('Fast method:',time.time()-start)

Answer 1

这是一个非常快速的方法：

\n

result = ((((12 * n + 45) * n + 50) * n + 15) * n - 2) * n // 120\n

\n

我是如何到达那里的：

\n

\n
1. 将内部总和重写为众所周知的x*(x+1)//2。于是整个事情就变成了sum(x**2 * x*(x+1)//2 for x in range(n+1))。
\n
2. 重写为sum(x**4 + x**3 for x in range(n+1)) // 2.
\n
3. 查找和的公式。sum(x**4)sum(x**3)
\n
4. 将由此产生的混乱简化(12*n**5 + 45*n**4 + 50*n**3 + 15*n**2 - 2*n) // 120为.
\n
5. 霍纳吧。
\n

\n

如果在步骤 1. 和 2. 之后您知道它是 5 次多项式，则可以采用另一种推导方法：

\n

\n
3. 通过简单的实现计算六个值。
\n
4. 根据具有六个未知数（多项式系数）的六个方程计算多项式。我的做法与此类似，但A与此相比，我的矩阵是左右镜像的，我称之为 y-vector b。
\n

\n

代码：

\n

from fractions import Fraction\nimport math\nfrom functools import reduce\n\ndef naive(n):\n    return sum(x**2 * sum(range(x+1)) for x in range(n+1))\n\ndef lcm(ints):\n    return reduce(lambda r, i: r * i // math.gcd(r, i), ints)\n\ndef polynomial(xys):\n    xs, ys = zip(*xys)\n    n = len(xs)\n    A = [[Fraction(x**i) for i in range(n)] for x in xs]\n    b = list(ys)\n    for _ in range(2):\n        for i0 in range(n):\n            for i in range(i0 + 1, n):\n                f = A[i][i0] / A[i0][i0]\n                for j in range(i0, n):\n                    A[i][j] -= f * A[i0][j]\n                b[i] -= f * b[i0]\n        A = [row[::-1] for row in A[::-1]]\n        b.reverse()\n    coeffs = [b[i] / A[i][i] for i in range(n)]\n    denominator = lcm(c.denominator for c in coeffs)\n    coeffs = [int(c * denominator) for c in coeffs]\n    horner = str(coeffs[-1])\n    for c in coeffs[-2::-1]:\n        horner += \' * n\'\n        if c:\n            horner = f"({horner} {\'+\' if c > 0 else \'-\'} {abs(c)})"\n    return f\'{horner} // {denominator}\'\n\nprint(polynomial((x, naive(x)) for x in range(6)))\n

\n

输出（在线尝试！）：

\n

((((12 * n + 45) * n + 50) * n + 15) * n - 2) * n // 120\n

\n

Answer 2

（最快的方法 3 和 4 在最后）

在快速 NumPy 方法中，您需要指定，dtype=np.object以便 NumPy 不会将 Python 转换int为其自己的数据类型（np.int64或其他数据类型）。现在它会给你正确的结果（检查到 N=100000）。

# method #2
start=time.time()
w=np.arange(0, n+1, dtype=np.object)
result2 = (w**2*np.cumsum(w)).sum()
print('Fast method:', time.time()-start)

您的快速解决方案比慢速解决方案要快得多。是的，对于较大的 N，但在 N=100 时，速度快了 8 倍：

start=time.time()
for i in range(100):
    result1 = summation(0, n, mysum)
print('Slow method:', time.time()-start)

# method #2
start=time.time()
for i in range(100):
    w=np.arange(0, n+1, dtype=np.object)
    result2 = (w**2*np.cumsum(w)).sum()
print('Fast method:', time.time()-start)

Slow method: 0.06906533241271973
Fast method: 0.008007287979125977

编辑：更快的方法（由KellyBundy，南瓜）是使用纯Python。事实证明，NumPy 在这里没有优势，因为它没有np.objects.

# method #3
import itertools
start=time.time()
for i in range(100):
    result3 = sum(x*x * ysum for x, ysum in enumerate(itertools.accumulate(range(n+1))))
print('Faster, pure python:', (time.time()-start))

Faster, pure python: 0.0009944438934326172

EDIT2：Forss 注意到 numpy 快速方法可以通过使用x*x而不是 来优化x**2。因为N > 200它比纯Python方法更快。因为N < 200它比纯Python方法慢（边界的确切值可能取决于机器，在我的机器上是200，最好自己检查一下）：

# method #4
start=time.time()
for i in range(100):
    w = np.arange(0, n+1, dtype=np.object)
    result2 = (w*w*np.cumsum(w)).sum()
print('Fast method x*x:', time.time()-start)

Answer 3

像这样将 Python 与 WolframAlpha 进行比较是不公平的，因为 Wolfram 在计算之前会简化方程。

\n

幸运的是，Python 生态系统没有限制，因此您可以使用SymPy：

\n

from sympy import summation\nfrom sympy import symbols\n\nn, x, y = symbols("n,x,y")\neq = summation(x ** 2 * summation(y, (y, 0, x)), (x, 0, n))\neq.evalf(subs={"n": 1000})\n

\n

它几乎会立即计算出预期结果：100375416791650。这是因为 SymPy 为您简化了方程，就像 Wolfram 一样。查看 的值eq：

\n

\n

@Kelly Bundy\'s 的答案很棒，但如果你像我一样使用计算器来计算2 + 2，那么你会喜欢 SymPy \xe2\x9d\xa4。正如您所看到的，只需 3 行代码即可获得相同的结果，并且该解决方案也适用于其他更复杂的情况。

\n