The sum of kth Powers

一道利用 Lagrange 插值法解决的多项式问题。

题源:Codeforces Edu Round 7 F

题解

不妨记 $p(x)=\sum\limits_{i=1}^x i^k$ ,不难发现,这是一个关于 $x$ 的 $k+1$ 次多项式。

(自然数幂次和公式如下,但发现上述性质只需仔细观察题目。)
$$
S_p(n) = 1^p + 2^p + \cdots+n^p= \frac{1}{p+1} \sum_{k=0}^{p} \binom{p+1}{k} B_k n^{p+1-k}
$$


我们希望能够求出 $p(n)$ 。

考虑 Lagrange 插值法,我们只需要确定 $k+2$ 个点,就能完全确定这个 $k+1$ 次多项式。

我们可以选取 $1,2,\cdots,k+2$ 这连续的 $k+2$ 个点,求出这些点的坐标所需时间复杂度为 $O(k\log k)$。


对于连续点的 Lagrange 插值,可以做如下优化:
$$
p(x)=\sum\limits_{i=1}^{k+2}y_i\prod\limits_{j=1,j\ne i}^{k+2}\frac{x-x_j}{x_i-x_j}=\sum\limits_{i=1}^{k+2}y_i\prod\limits_{j=1,j\ne i}^{k+2}\frac{x-j}{i-j}
$$
连乘项中,分子维护住前缀积 $pre[i]=\prod\limits_{j=1}^{i}(x-j)$,后缀积 $suf[i]=\prod\limits_{j=i}^{k+2}(x-j)$。

分母维护住阶乘前缀积 $pref[i]=\prod\limits_{j=1}^ij$,这样,我们的连乘项可以化成:
$$
\prod\limits_{j=1,j\ne i}^{k+2}\frac{x-j}{i-j}=\frac{pre[i-1]\times suf[i+1]}{pref[i-1]\times (-1)^{k+2-i}\cdot pref[k+2-i]}
$$
这样,就能在 $O(k)$ 时间完成 $p(n)$ 的计算。


总复杂度 $O(k\log k)$,瓶颈在于获取前 $k+2$ 项的初值。代码实现:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#include<bits/stdc++.h>
#define int long long
using namespace std;

constexpr int mod = 1e9 + 7;

int qmi(int a, int k, int p){
int res = 1;
while(k){
if(k & 1) res = res * a % p;
k >>= 1;
a = a * a % p;
}
return res;
}

int inv(int x){
return qmi(x, mod - 2, mod);
}

struct modint{
int v;
modint(int x = 0) : v((x % mod + mod) % mod) {}

modint operator+ (const modint &o){
return modint(v + o.v);
}
modint operator- (const modint &o){
return modint(v - o.v);
}
modint operator* (const modint &o){
return modint(v * o.v);
}
modint operator/ (const modint &o){
return modint(v * inv(o.v));
}
};
using Z = modint;

void solve(){
int n, k; cin >> n >> k;
vector<Z>p(k + 3);
p[1] = 1;
for(int i = 2; i <= k + 2; i++){
p[i] = p[i - 1] + Z(qmi(i, k, mod));
}

if(n <= k + 2){
cout << p[n].v << '\n';
return;
}

vector<Z>pre(k + 3), suf(k + 4);
pre[0] = 1;
for(int i = 1; i <= k + 2; i++){
pre[i] = pre[i - 1] * Z(n - i);
}
suf[k + 3] = 1;
for(int i = k + 2; i >= 1; i--){
suf[i] = suf[i + 1] * Z(n - i);
}

vector<Z>pref(k + 3);
pref[0] = 1;
for(int i = 1; i <= k + 2; i++){
pref[i] = pref[i - 1] * Z(i);
}

Z res(0);
for(int i = 1; i <= k + 2; i++){
Z mo = pre[i - 1] * suf[i + 1];
Z de = pref[i - 1] * pref[k + 2 - i] * Z((k + 2 - i) % 2 ? -1 : 1);
res = res + p[i] * mo / de;
}
cout << res.v << '\n';
}

signed main(){
ios::sync_with_stdio(0);
cin.tie(0);
solve();
return 0;
}