本题是本人目前出的题中难度最高的题。
加强版紫色是肯定可以的。
题目链接：无所畏惧的求和 - 洛谷
尽情享受吧！
这道题其实做法有很多：
待定系数法 + 矩阵求解 推代数公式（不够优秀）
组合数学 + 待定系数法 推组合公式（受限于 \(k\) 的范围）
在模意义下的特殊做法（受限于模数范围）
拉格朗日插值法（正解）
第一类斯特林数（时间复杂度可能有点问题，为 \(O(k^2)\)
伯努利数（奇奇怪怪）
……
那么这里讲解前四种方式
自然数幂方求和公式？在高等教育出版社出版的《数学手册》中有这么一些公式：
采用万能的数学归纳法可以一一证明上述公式。此处不提。
但是观察上述公式，可以发现一个特征：自然数 \(k\) 次幂求和公式是关于 \(n\) 的 \(k + 1\) 次有理多项式。
也就是
知道上述结果之后，可采用待定系数法，也就是写出 \(n = 1, 2, 3, \dots, k+1\) 的 \(k + 1\) 个代数式，利用矩阵求解即可。
举个例子，对于 \(k = 3\) 的情况：
那么借此求出每一项的系数即可对于每一个询问在 \(O(k)\) 的复杂度内完成计算。
在获取系数的部分：
vector<int> mat[K];
void get_coefficient(int k, vector<int> &cctk) {
int h = k + 1, w = k + 2;
cctk.assign(h + 1, 0);
for (int i(1); i <= h; ++i) mat[i].assign(w + 1, 0);
int sum(0);
for (int n(1); n <= h; ++n) {
(sum += qpow(n, k)) %= mod;
for (int x(n), j(1); j < w; ++j, (x *= n) %= mod) {
mat[n][j] = x;
} mat[n][w] = sum;
}
// Guass 消元
for (int i(1); i <= h; ++i) {
int inv = qpow(mat[i][i], mod - 2);
for (int j(i); j <= w; ++j)
(mat[i][j] *= inv) %= mod;
for (int j = 1; j <= h; ++j) {
if (j ^ i) {
inv = mat[j][i];
for (int c(i); c <= w; ++c) {
(mat[j][c] -= mat[i][c] * inv) %= mod;
if (mat[j][c] < 0) mat[j][c] += mod;
}
}
}
}
for (int i = 1; i <= h; ++i) {
cctk[i] = mat[h + 1 - i][w];
}
}
在最终求答案的部分：
void work() {
int T;
cin >> T;
while (T--) {
int n, k;
cin >> n >> k;
vector<int> &cctk = ccts[k];
if (cctk.empty())
get_coefficient(k, cctk);
int ans(0);
n %= mod;
for (int i = 1; i <= k + 1; ++i) {
(ans = ans * n % mod + cctk[i]) %= mod;
}
cout << ans * n % mod << '\n';
}
}
总的复杂度为 \(O(k^4 + T \cdot k)\)。还不够优秀，普通版只能过前10个点，期望得分30，加强版无法得分。
但是简单版更快……
这个方法相对优秀一点。标程就是用的此写法。
其实不难发现，对于 \(x^k\)，我们可以改写为：
那么依据某些公式推导：
所以，类似的，我们也可以枚举 \(n = 1, 2, 3, \dots, k\) 来寻找其系数：
以 \(k = 3\) 为例
同时我们规定 \(\binom nr = 0\ (n \lt r)\)。所以上式也可以写为
这就是一个下三角矩阵，每一次扫一遍即可。
代码也非常简单，常数比第一种方法小很多：
void get_coefficient(int k, int * ccts) {
int sum = 0;
for (int n = 1; n <= k; ++n) {
sum += qpow(n, k, MOD);
int rest = sum;
// 由于我们已经知道了前 (n-1) 个系数，直接通过总和一一减去即可。
for (int i = 1; i < n; ++i) {
// 注意 k < MOD，所以此处不需要Lucas
(rest -= ccts[i] * C(n + 1, i + 1) % MOD) %= MOD;
}
if (rest < 0) rest += MOD;
// 明显可知，C(n, n) = 1，所以剩下的即是系数
ccts[n] = rest;
}
}
核心部分也非常简单，只是模数较小，需要用到 \(Lucas\) 定理。
int ccts[K][K]; // 用于保存系数
void work() {
int T, n, k;
cin >> T;
while (T--) {
cin >> n >> k;
int * cctk = ccts[k];
if (!cctk[1]) // 其实不难发现，a1 一定为 1，所以借此判断即可
get_coefficient(k, cctk);
int ans = 0;
for (int i = 1; i <= k; ++i)
(ans += cctk[i] * Lucas(n + 1, i + 1)) %= MOD;
cout << ans << '\n';
}
}
总的复杂度为 \(O(k^3 + T \cdot \log_pn)\)，由于 \(p = 111121\)，所以我们可以认为复杂度为 \(O(2k^3 + 2T)\)。普通版期望得分100分，加强版受限于 \(k\) 的大小，期望得分30分。
由于数据原因需要取模……所以有了这么一个做法。
做法来源于 fanghaoyu
在模意义下有如下恒等式：
意味着
那么，很明显，我们可以预处理出 \(\lfloor \frac np \rfloor \sum_{x=1}^p x^k\)
然后对于每一次求解求后一项即可，每一次查询时间为 \(O(p \log k)\)。
当然，也可以采用空间换时间的方法，保存后一项，使得查询时间复杂度为 \(O(1)\)。
而预处理需要 \(O(p)\)，所以总的复杂度大概为 \(O(p k\log k + T p\log k)\) 或者 \(O(p k\log k + T)\)
而空间复杂度……
其实这个方法特别好卡，当初我为了提示使用 Lucas 定理，专门设置的小质数，也就成了这个做法的来源，所以只需要小小的修改一下模数，改大一点，这就熄火了 QwQ。
为了尊重这个方法的想出者，所以贴出未修改的源代码。
#include <bits/stdc++.h>
using namespace std;
const int N = 1e5 + 5,MOD = 8887;
typedef long long ll;
ll n,k,a[1001][MOD + 5];
inline ll ksm(ll base,ll pts)
{
ll ret = 1;
for(;pts > 0;pts >>= 1,base = base * base % MOD)
if(pts & 1)
ret = ret * base % MOD;
return ret;
}
int main()
{
ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
int T;
cin>>T;
for(int t = 1;t <= 1000;t++)
for(int i = 1;i <= MOD;i++)
a[t][i] = a[t][i - 1] + ksm(i,t),a[t][i] %= MOD;
while(T--)
{
cin>>n>>k;
ll ans = 0;
ans = ans + a[k][MOD] * (n / MOD) % MOD;
ans = ans + a[k][n % MOD];
cout<<ans<<endl;
}
return 0;
}
在初中学过，两点确定一线（一次多项式），三点可以确定一个抛物线（二次多项式）。扩展出来，\(n + 1\) 个点可以确定 \(n\) 次多项式。
在上文中提及，自然数 \(k\) 次幂求和公式是关于 \(n\) 的 \(k + 1\) 次有理多项式。
那么我们只需要确定 \(k + 2\) 个点，就可以求出这个多项式。加之我们只需要求在 \(n\) 处的值，所以想到使用拉格朗日插值定理求解。
首先，对于一个 \(k\) 次多项式，我们有了 \(k + 1\) 个点 \((x_i, y_i)\)。
定义：
容易构造出
因此最终我们得到：
那么考虑为什么？
由于我们对于每一个 \(y_if_i(x)\)，可以保证在 \(x_i\) 点处，取值为 \(y_i\)，其余的点上取值为 \(0\)。
那么 \(f(x) = \sum_{i=1}^{k+1} y_i fi(x)\) 可以一一穿过这三个点。
也就是说，我们可以借此求得 \(f(x)\) 了。
更严谨的证明此处略过……
很明显，对于每一个分母，是两个阶乘的积，所以预处理阶乘。
对于分母，预处理前缀积和后缀积即可。
核心代码如下：
typedef long long lint;
lint fac[K], rfac[K];
inline lint qpow(lint a, int x) {
lint r = 1;
for (; x; x >>= 1, a = a * a % mod)
if (x & 1) r = r * a % mod;
return r;
}
inline void prepare() {
fac[0] = rfac[0] = 1;
// 预处理下方阶乘
for (int i = 1; i < K; ++i) {
fac[i] = fac[i - 1] * i % mod;
// rfac[i] = rfac[i - 1] * (-i) % mod;
rfac[i] = rfac[i - 1] * (mod - i) % mod;
}
}
// prefix, suffix
lint pfx[K], sfx[K];
int calc(int n, int k) {
pfx[0] = sfx[k + 3] = 1;
for (int i = 1; i <= k + 2; ++i)
pfx[i] = pfx[i - 1] * (n - i) % mod;
for (int i = k + 2; i; --i)
sfx[i] = sfx[i + 1] * (n - i) % mod;
lint y = 0, ans = 0, frup, frdown;
for (int i = 1; i <= k + 2; ++i) {
y = (y + qpow(i, k)) % mod;
frup = pfx[i - 1] * sfx[i + 1] % mod;
frdown = fac[i - 1] * rfac[k + 2 - i] % mod;
ans = (ans + y * frup % mod * qpow(frdown, mod - 2) % mod) % mod;
}
return (ans + mod) % mod;
}
时间复杂度：\(O(\sum k)\)，期望得分 100。