有没有算法来计算(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)，我只能解决m = 1000000007。它非常慢，因为时间限制是2秒。
你有任何有效的算法吗？
有评论说它可能是this question的副本，但它肯定是不同的。
你可以总结一下这个系列
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）。
有几种方法可以加速模幂运算。从这里开始，我将使用**来表示“指数”，并使用%来表示“模数”。
先是几点意见。总是这样的情况，(a * b) % m是((a % m) * (b % m)) % m。 a ** n与(a ** floor(n / 2)) * (a ** (n - floor(n/2))相同也是如此。这意味着对于指数<= 1000，我们总是可以在最多20次乘法（和21次mod）中完成取幂。
我们也可以跳过相当多的计算，因为(a ** b) % m与((a % m) ** b) % m相同，如果m明显低于n，我们只需要多个重复和，具有部分重复的“尾部”。
我认为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的有用链接。
如果你想要一些易于实现和快速的东西，试试这个：
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 = 0cj = (j)_(i - 1)并且ith方程的独立值是(n)_(i - 1)。解释为什么有点凌乱，但你可以检查证明here。
上述算法相当于求解该系统的线性方程组。在https://github.com/fcard/PolySum中可以找到大量的实现和进一步的解释。这个算法的主要缺点是它消耗了大量的内存，即使我的大多数内存效率版本使用几乎1gb为n=3000。但它比SymPy和Mathematica都要快，所以我认为没关系。与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))，因此它使用前一行来计算当前行。