有什么算法可以计算(1^x + 2^x + 3^x + ... + n^x) mod 1000000007吗？
注：a^b是 a 的 b 次方。
约束是1 <= n <= 10^16, 1 <= x <= 1000。所以N的值很大。
我只能解决O(m log m)if m = 1000000007。它非常慢，因为时间限制是 2 秒。
你有什么有效的算法吗？
有评论说它可能与这个问题重复，但肯定是不同的。
问问题
2808 次
4 回答
19
你可以总结一下这个系列
1**X + 2**X + ... + N**X
借助Faulhaber 公式，您将得到一个具有X + 1计算任意N.
如果您不想计算Bernoulli Numbers，您可以通过求解X + 2 线性方程(for N = 1, N = 2, N = 3, ..., N = X + 2) 找到多项式，这是一种较慢但更易于实现的方法。
让我们举个例子X = 2。在这种情况下，我们有一个X + 1 = 3阶多项式：
A*N**3 + B*N**2 + C*N + D
线性方程是
A + B + C + D = 1 = 1
A*8 + B*4 + C*2 + D = 1 + 4 = 5
A*27 + B*9 + C*3 + D = 1 + 4 + 9 = 14
A*64 + B*16 + C*4 + D = 1 + 4 + 9 + 16 = 30
解决了我们将得到的方程
A = 1/3
B = 1/2
C = 1/6
D = 0
最后的公式是
1**2 + 2**2 + ... + N**2 == N**3 / 3 + N**2 / 2 + N / 6
现在，您所要做的就是在公式中输入任意大。 N到目前为止，该算法具有O(X**2)复杂性（因为它不依赖于N）。
于 2017-01-17T11:41:20.907 回答
3
有几种方法可以加快模幂运算。从这里开始，我将使用**“取幂”和%“模数”来表示。
首先是一些观察。(a * b) % m情况总是如此((a % m) * (b % m)) % m。情况也总是与a ** n相同(a ** floor(n / 2)) * (a ** (n - floor(n/2))。这意味着对于 <= 1000 的指数，我们总是可以在最多 20 次乘法（和 21 个模数）中完成幂运算。
我们也可以跳过很多计算，因为(a ** b) % m与 相同 ((a % m) ** b) % m，如果 m 明显小于 n，我们只是有多个重复的和，带有部分重复的“尾巴”。
于 2017-01-17T11:19:04.087 回答
1
我认为 Vatine 的答案可能是要走的路，但我已经输入了这个，它可能对这个或其他人的类似问题有用。
今天早上我没有时间详细回答，但请考虑一下。
1^2 + 2^2 + 3^2 + ... + n^2将采取 O( n ) 步骤直接计算。但是，它等价于(n(n+1)(2n+1)/6)，可以在 O(1) 时间内计算。对于任何更高的积分幂
x都存在类似的等价性。
此类问题可能有一个通用的解决方案；我不知道，而且 Wolfram Alpha 似乎也不知道。但是，一般来说，等价表达式是x+1次的，可以通过计算一些样本值并求解一组系数的线性方程来计算。
但是，对于较大的x（例如您的问题中的 1000 ），这将很困难，并且可能无法在 2 秒内完成。
也许比我了解更多数学的人可以把它变成一个可行的解决方案？
编辑：哎呀，我看到 Fabian Pijcke 在我发布之前已经发布了一个关于 Faulhaber公式的有用链接。
于 2017-01-17T11:27:57.753 回答
0
如果您想要一些易于实现且快速的东西，请尝试以下操作：
Function Sum(x: Number, n: Integer) -> Number
P := PolySum(:x, n)
return P(x)
End
Function PolySum(x: Variable, n: Integer) -> Polynomial
C := Sum-Coefficients(n)
P := 0
For i from 1 to n + 1
P += C[i] * x^i
End
return P
End
Function Sum-Coefficients(n: Integer) -> Vector of Rationals
A := Create-Matrix(n)
R := Reduced-Row-Echelon-Form(A)
return last column of R
End
Function Create-Matrix(n: Integer) -> Matrix of Integers
A := New (n + 1) x (n + 2) Matrix of Integers
Fill A with 0s
Fill first row of A with 1s
For i from 2 to n + 1
For j from i to n + 1
A[i, j] := A[i-1, j] * (j - i + 2)
End
A[i, n+2] := A[i, n]
End
A[n+1, n+2] := A[n, n+2]
return A
End
解释
我们的目标是获得一个多项式Q，使得Q(x) = sum i^n for i from 1 to x。知道Q(x) = Q(x - 1) + x^n=> Q(x) - Q(x - 1) = x^n，我们就可以像这样建立一个方程组：
d^0/dx^0( Q(x) - Q(x - 1) ) = d^0/dx^0( x^n )
d^1/dx^1( Q(x) - Q(x - 1) ) = d^1/dx^1( x^n )
d^2/dx^2( Q(x) - Q(x - 1) ) = d^2/dx^2( x^n )
... .
d^n/dx^n( Q(x) - Q(x - 1) ) = d^n/dx^n( x^n )
假设Q(x) = a_1*x + a_2*x^2 + ... + a_(n+1)*x^(n+1)，我们将n+1得到具有未知数 的线性方程a1, ..., a_(n+1)，结果证明方程中与未知数cj相乘的aj系数i遵循以下模式（其中(k)_p= (k!)/(k - p)!）：
- 如果
j < i,cj=0 - 否则，
cj=(j)_(i - 1)
i并且第 th 个方程的独立值为(n)_(i - 1)。解释为什么会有点混乱，但你可以在这里查看证明。
上述算法等价于求解这个线性方程组。大量的实现和进一步的解释可以在https://github.com/fcard/PolySum中找到。该算法的主要缺点是它消耗大量内存，即使是我最节省内存的版本也几乎1gb使用n=3000. 但它比 SymPy 和 Mathematica 都快，所以我认为没关系。与舒尔茨方法相比，后者使用一组备用方程。
例子
对于小型 . 手动应用此方法很容易n。的矩阵n=1是
| (1)_0 (2)_0 (1)_0 | | 1 1 1 |
| 0 (2)_1 (1)_1 | = | 0 2 1 |
应用 Gauss-Jordan 消元，我们得到
| 1 0 1/2 |
| 0 1 1/2 |
结果 = {a1 = 1/2, a2 = 1/2}=>Q(x) = x/2 + (x^2)/2
注意矩阵总是已经是行梯形的，我们只需要减少它。
对于n=2：
| (1)_0 (2)_0 (3)_0 (2)_0 | | 1 1 1 1 |
| 0 (2)_1 (3)_1 (2)_1 | = | 0 2 3 2 |
| 0 0 (3)_2 (2)_2 | | 0 0 6 2 |
应用 Gauss-Jordan 消元，我们得到
| 1 1 0 2/3 | | 1 0 0 1/6 |
| 0 2 0 1 | => | 0 1 0 1/2 |
| 0 0 1 1/3 | | 0 0 1 1/3 |
结果 = {a1 = 1/6, a2 = 1/2, a3 = 1/3}=>Q(x) = x/6 + (x^2)/2 + (x^3)/3
算法速度的关键在于它不会为矩阵的每个元素计算阶乘，而是知道(k)_p= (k)_(p-1) * (k - (p - 1))，因此A[i,j]= (j)_(i-1)= (j)_(i-2) * (j - (i - 2))= A[i-1, j] * (j - (i - 2))，因此它使用前一行来计算当前行。
于 2019-03-26T07:22:15.150 回答