Python中的高效求和

Kel*_*ndy 53

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

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

我是如何到达那里的：

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

如果在步骤 1 和 2 之后得出它的另一种方法。您知道它是 5 次多项式：

用简单的实现计算六个值。 从具有六个未知数（多项式系数）的六个方程计算多项式。我的做法与此类似，但与此相比，我的矩阵A是左右镜像的，我称其为 y-vector b。

代码：

from fractions import Fraction
import math
from functools import reduce

def naive(n):
    return sum(x**2 * sum(range(x+1)) for x in range(n+1))

def lcm(ints):
    return reduce(lambda r, i: r * i // math.gcd(r, i), ints)

def polynomial(xys):
    xs, ys = zip(*xys)
    n = len(xs)
    A = [[Fraction(x**i) for i in range(n)] for x in xs]
    b = list(ys)
    for _ in range(2):
        for i0 in range(n):
            for i in range(i0 + 1, n):
                f = A[i][i0] / A[i0][i0]
                for j in range(i0, n):
                    A[i][j] -= f * A[i0][j]
                b[i] -= f * b[i0]
        A = [row[::-1] for row in A[::-1]]
        b.reverse()
    coeffs = [b[i] / A[i][i] for i in range(n)]
    denominator = lcm(c.denominator for c in coeffs)
    coeffs = [int(c * denominator) for c in coeffs]
    horner = str(coeffs[-1])
    for c in coeffs[-2::-1]:
        horner += ' * n'
        if c:
            horner = f"({horner} {'+' if c > 0 else '-'} {abs(c)})"
    return f'{horner} // {denominator}'

print(polynomial((x, naive(x)) for x in range(6)))

输出（在线尝试！）：

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

@Adam是的，我认为在这种情况下，您已经简化了实际问题，在问题中解释给定公式只是一个示例非常重要，但您的目标实际上是弄清楚如何快速计算总和，而不是得到该特定公式的实际答案。否则，您将面临获得类似解决方案的风险，这无疑是解决您提出的问题的最佳方法，但对您真正遇到的问题完全没有帮助。 (15认同) @Adam 我想这解释了为什么您在非 NumPy 解决方案中使用了所有这些函数，这看起来确实很奇怪。也许更一般的情况仍然允许类似的优化，但这取决于更一般的程度。也许你真的有通用公式问另一个问题？就像，使用 `f(x)` 和 `g(y)` 而不是 `x^2` 和 `y` 左右，其中 `f` 和 `g` 是未知函数（尽管可能某些属性是已知的并且可以被利用）。 (2认同)

dan*_*444 18

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

在快速 NumPy 方法中，您需要指定dtype=np.object以便 NumPy 不会将 Python 转换int为其自己的 dtypes（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.aumulate(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)

我也会尝试等效的非 NumPy 版本，您可能会发现它比 NumPy 版本*快*。例如`result1 = sum(x*x * ysum for x, ysum in enumerate(itertools.aumulate(range(n+1))))`或`ysum = 0; result1 = sum(x*x * (ysum := ysum + x) for x in range(n+1))` (6认同) @diggusbickus 因为`np.int64`只有64位来存储整数，Python`int`可以和你的RAM一样大。通过使用“通用”`np.object`，您可以确保 numpy 不会将 `int` 转换为 `np.int64`。 (4认同) 那为什么 np.object 而不是 np.int64 呢？ (3认同) 我不认为它适合我的（因为我只是在谈论我的方法并希望保持这种方式），但它会让你的方法更好。 (3认同) 纯 python 版本在比较中有点作弊，使用 `x*x` 而不是 `x**2` 像其他方法一样。将 numpy 解决方案更改为 `x*x` 是较大 n 的最快方法（在我的计算机上）。 (3认同) @Forss 嗯，对。我出于习惯使用了`x*x`，尽管部分原因是速度。既然你指出了，我很失望 NumPy 并没有一次发现它是一次乘法，而是将顿悟应用于整个数组。但后来我记得我们正在使用“对象”，我猜 NumPy 并没有在那里做出假设/类型分析。如果没有 `dtype`，它会同样快地执行 `x*x` 和 `x**2`。对于使用 `x*x` 的两种解决方案，NumPy 在 n=1000 时对我来说快一点，在 n=10000 时差不多快，在 n=100000 时慢一点。 (2认同)

Peq*_*que 7

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

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

from sympy import summation
from sympy import symbols

n, x, y = symbols("n,x,y")
eq = summation(x ** 2 * summation(y, (y, 0, x)), (x, 0, n))
eq.evalf(subs={"n": 1000})

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

@Kelly Bundy 的回答很棒，但是如果您像我一样使用计算器进行计算2 + 2，那么您会喜欢 SymPy 吗？如您所见，只需 3 行代码即可获得相同的结果，并且该解决方案也适用于其他更复杂的情况。