计算 1^X + 2^X + ... + N^X mod 1000000007
Calculating 1^X + 2^X + ... + N^X mod 1000000007
有什么算法可以计算(1^x + 2^x + 3^x + ... + n^x) mod 1000000007吗？
注：a^b是a的b次方。
限制条件为 1 <= n <= 10^16, 1 <= x <= 1000。所以N的值很大
如果m = 1000000007，我只能解决O(m log m)。很慢，因为时间限制是2秒。
你有什么高效的算法吗？
有人评论说它可能与 this question 重复，但绝对不同。
有几种方法可以加快模幂运算。从这里开始，我将用**表示"exponentiate"，用%表示"modulus"。
首先是一些观察。 (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，我们只是有多个重复的和，"tail" 部分重复。
我认为 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's formula 的有用 link。
你可以总结这个系列
1**X + 2**X + ... + N**X
在 Faulhaber's formula 的帮助下，您将得到一个具有 X + 1 幂的多项式来计算任意 N。
如果您不想计算伯努利数，您可以通过求解X + 2线性方程来求多项式]（对于 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）。
如果你想要一些易于实施和快速的东西，试试这个：
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)!）：
- if
j < i, cj = 0
- 否则，
cj = (j)_(i - 1)
第i个方程的独立值为(n)_(i - 1)。解释为什么有点乱，但你可以检查证明 here.
上述算法相当于求解这个线性方程组。
在 https://github.com/fcard/PolySum. The main drawback of this algorithm is that it consumes a lot of memory, even my most memory efficient version uses almost 1gb for n=3000. But it's faster than both SymPy and Mathematica, so I assume it's okay. Compare to Schultz's method 中可以找到大量的实现和进一步的解释，它使用了一组备用方程。
例子
对于小 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))，因此它使用前一行来计算当前行。
有什么算法可以计算(1^x + 2^x + 3^x + ... + n^x) mod 1000000007吗？
注：a^b是a的b次方。
限制条件为 1 <= n <= 10^16, 1 <= x <= 1000。所以N的值很大
如果m = 1000000007，我只能解决O(m log m)。很慢，因为时间限制是2秒。
你有什么高效的算法吗？
有人评论说它可能与 this question 重复，但绝对不同。
有几种方法可以加快模幂运算。从这里开始，我将用**表示"exponentiate"，用%表示"modulus"。
首先是一些观察。 (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，我们只是有多个重复的和，"tail" 部分重复。
我认为 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's formula 的有用 link。
你可以总结这个系列
1**X + 2**X + ... + N**X
在 Faulhaber's formula 的帮助下，您将得到一个具有 X + 1 幂的多项式来计算任意 N。
如果您不想计算伯努利数，您可以通过求解X + 2线性方程来求多项式]（对于 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）。
如果你想要一些易于实施和快速的东西，试试这个：
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)!）：
- if
j < i,cj=0 - 否则，
cj=(j)_(i - 1)
第i个方程的独立值为(n)_(i - 1)。解释为什么有点乱，但你可以检查证明 here.
上述算法相当于求解这个线性方程组。
在 https://github.com/fcard/PolySum. The main drawback of this algorithm is that it consumes a lot of memory, even my most memory efficient version uses almost 1gb for n=3000. But it's faster than both SymPy and Mathematica, so I assume it's okay. Compare to Schultz's method 中可以找到大量的实现和进一步的解释，它使用了一组备用方程。
例子
对于小 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))，因此它使用前一行来计算当前行。