题意
传送门 Codeforces 622F The Sum of the k-th Powers
题解
∑
i
=
1
n
i
k
\sum\limits_{i=1}^{n}i^k
i=1∑nik 是一个关于
n
n
n 的
k
+
1
k + 1
k+1 次多项式。用
1
,
2
,
3
⋯
k
+
2
1,2,3\cdots k + 2
1,2,3⋯k+2 共
k
+
2
k+2
k+2 个值,进行拉格朗日插值。令
f
(
n
)
=
∑
i
=
1
n
i
k
f(n) = \sum\limits_{i=1}^{n}i^k
f(n)=i=1∑nik。得到
f
(
x
)
=
∑
i
=
1
k
+
2
f
(
i
)
∏
j
=
1
,
j
≠
i
k
+
2
x
−
j
i
−
j
f(x) = \sum\limits_{i=1}^{k+2}f(i)\prod\limits_{j=1,j\neq i}^{k+2}\frac{x-j}{i-j}
f(x)=i=1∑k+2f(i)j=1,j=i∏k+2i−jx−j 分子可以表示为
∏
j
=
1
k
+
2
x
−
j
x
−
i
\frac{\prod\limits_{j=1}^{k+2}x - j}{x - i}
x−ij=1∏k+2x−j 可以通过预处理前缀与后缀积
O
(
1
)
O(1)
O(1) 得到。分母是
(
i
−
1
)
!
(
k
+
2
−
i
)
!
(
−
1
)
k
+
2
−
i
(i-1)!(k+2-i)!(-1)^{k+2-i}
(i−1)!(k+2−i)!(−1)k+2−i 可以通过预处理阶乘的逆元
O
(
1
)
O(1)
O(1) 得到。对于
f
(
i
)
f(i)
f(i),满足
f
(
i
)
=
f
(
p
)
f
(
i
/
p
)
,
p
∣
i
f(i) = f(p)f(i/p),p\vert i
f(i)=f(p)f(i/p),p∣i,可以通过线性筛递推。由于
[
1
,
n
]
[1,n]
[1,n] 内质数规模为
O
(
n
/
ln
n
)
O(n/\ln n)
O(n/lnn),且仅当筛到质数时进行快速幂,故求解
f
(
i
)
f(i)
f(i) 总时间复杂度接近于线性。
总时间复杂度 O ( k ) O(k) O(k)。
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
constexpr int MOD = (int)1e9 + 7;
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
int n, k;
cin >> n >> k;
int mx = k + 2;
vector<ll> fac(mx + 1), inv(mx + 1), invf(mx + 1);
fac[0] = fac[1] = 1;
inv[1] = 1;
invf[0] = invf[1] = 1;
for (int i = 2; i <= mx; ++i) {
fac[i] = fac[i - 1] * i % MOD;
inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD;
invf[i] = invf[i - 1] * inv[i] % MOD;
}
vector<ll> pre(mx + 2), suf(mx + 2);
pre[0] = suf[mx + 1] = 1;
for (int i = 1; i <= mx; ++i) {
pre[i] = pre[i - 1] * (n - i) % MOD;
}
for (int i = mx; i > 0; --i) {
suf[i] = suf[i + 1] * (n - i) % MOD;
}
auto power = [&](ll x, int n) {
ll res = 1;
while (n > 0) {
if (n & 1) (res *= x) %= MOD;
(x *= x) %= MOD, n >>= 1;
}
return res;
};
vector<ll> f(mx + 1);
vector<int> prime, minp(mx + 1);
f[1] = 1;
for (int i = 2; i <= mx; ++i) {
if (minp[i] == 0) {
minp[i] = i;
f[i] = power(i, k);
}
for (int j = 0, x; j < (int)prime.size() && (x = prime[j] * i) <= mx; ++j) {
minp[x] = prime[j];
f[x] = f[i] * f[prime[j]] % MOD;
if (minp[i] == prime[j]) break;
}
}
for (int i = 2; i <= mx; ++i) (f[i] += f[i - 1]) %= MOD;
ll res = 0;
for (int i = 1; i <= mx; ++i) {
ll up = pre[i - 1] * suf[i + 1] % MOD;
ll down = invf[i - 1] * invf[k + 2 - i] % MOD;
if ((k + 2 - i) & 1) down = MOD - down;
(res += up * down % MOD * f[i] % MOD) %= MOD;
}
res = (res + MOD) % MOD;
cout << res << '\n';
return 0;
}