问题
给定 $n, k$,求出:
1^k+2^k+\cdots+n^k=\sum_{i=1}^{n}{i^k}\bmod 1000000007
$$
$n\leq 10^9$,$k\leq 10^6$.
解决这个问题需要用到一些简单的多项式技巧。
拉格朗日插值
如果给定一个 $k$ 次多项式上的 $k+1$ 个点值,就可以通过拉格朗日差值求出这个多项式的表达式:
f(x)=\sum_{i=1}^{k+1}y_i\prod_{j\neq i}{\frac{x-x_j}{x_i-x_j}}
$$
对于每个点,我们构造一个 $k$ 次多项式 $g_i(x)$,满足:
g_i(x_t)=
\left\{
\begin{align}
&0 &t=i\\
&y_t &t\neq i
\end{align}
\right.
$$
显然这样的多项式可以为:
g_i(x) = y_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}
$$
我们把这 $k+1$ 个 $g$ 加起来,每个点值只会被加上一次,所以得到:
f(x)=\sum_{i=1}^{k+1}g_i(x)=\sum_{i=1}^{k+1}y_i\prod_{j\neq i}{\frac{x-x_j}{x_i-x_j}}
$$
$k$ 次幂和是 $k+1$ 次多项式
$f_0(n) =n$.
$f_1(n) =\dfrac{n(n+1)}{2}$.
$f_2(n)=\dfrac{n(n+1)(2n+1)}{6}$.
$f_3(n)=\dfrac{n^2(n+1)^2}{4}$.
众所周知,$f_k(n)=\sum_{i=1}^{n}{i^k}$,是一个 $k+1$ 次多项式,可以用数学归纳法证明。
分析
回到题目,设 $S_k(n)$ 为题目所求。由上面的性质,$S_k(n)$ 是 $k+1$ 次多项式,需要 $k+2$ 个点值进行构造。
如果直接选取 $k+1$ 个点,利用拉格朗日插值计算,复杂度为 $\mathcal{O}(k^2)$,无法通过。
优化
考虑选连续的一段点,进行优化。
我们选择 $k+2$ 个自变量,为 $[1, k+2]$ 中的整数。也就是说选择的点集为:$\{(i, S_k(i))\mid i\in \Bbb{N}^{+}, i\leq k+2\}$.
将这些点代入拉格朗日插值公式,得到:
S_k(n)=\sum_{i=1}^{k+2}{S_k(i)\prod_{j\neq i}{\frac{n-j}{i-j}}}
$$
简单变换后得到:
S_k(n)=\sum_{i=1}^{k+2}{S_k(i){\frac{\prod\limits_{j\neq i}(n-j)}{\prod\limits_{j\neq i}(i-j)}}}
$$
之后开始找规律)
对于累乘项的分子,暴力拆开:
\begin{align}
\prod_{j\neq i}(n-j)&= n(n-1)\dots (n-(i-1))\times(n-(i+1))\dots(n-(k+2))\\
&= (\prod_{j=1}^{i-1}(n-j))(\prod_{j=i+1}^{k+2}(n-j))
\end{align}
$$
所以我们维护 $(n-j)$ 的前缀积和后缀积,就能直接得到分子的值了。
对于分母,假设 $i$ 已经固定,那么有:
\begin{aligned}
&\dfrac{1}{\prod\limits_{i\not ={j}}(i-j)}\\ =
&\dfrac{1}{i(i-1)(i-2) \dots 1 \times(-1) \dots (k+2-i-1)(k+2-i)}\\=&(-1)^{k+2-i}\dfrac{1}{i! (k+2-i)!}
\end{aligned}
$$
所以我们预处理阶乘,再求逆元,就可以快速得到分母的值。
最后,把变换后的式子代入 $S_k(n)$ 的表达式:
S_k(n)=\sum_{i=1}^{k+2}S_k(i)(-1)^{k+2-i}\dfrac{(\prod\limits_{j=1}^{i-1}(n-j))(\prod\limits_{j=i+1}^{k+2}(n-j))}{i!(k+2-i)!}
$$
时间复杂度是 $\mathcal{O}(k\log k)$.
代码
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int MAXK = 1e6 + 10, mod = 1e9 + 7;
int n, k;
int pre[MAXK], suf[MAXK], fac[MAXK];
int power(int x, int k) {
int ret = 1;
while(k) {
if(k & 1) ret = ret * x % mod;
k >>= 1;
x = x * x % mod;
}
return ret;
}
inline int inv(int x) {
return power(x, mod - 2);
}
signed main() {
scanf("%lld%lld", &n, &k);
fac[0] = pre[0] = 1;
for(int i = 1; i <= k + 2; i++) {
pre[i] = pre[i - 1] * (n - i) % mod;
fac[i] = fac[i - 1] * i % mod;
}
suf[k + 3] = 1;
for(int i = k + 2; i >= 1; i--) suf[i] = suf[i + 1] * (n - i) % mod;
int tmp = 0, ans = 0;
for(int i = 1; i <= k + 2; i++) {
tmp = (tmp + power(i, k)) % mod;
int a = pre[i - 1] * suf[i + 1] % mod;
int b = fac[i - 1] * fac[k + 2 - i] % mod;
if((k - i) % 2 == 0) ans = (ans + tmp * a % mod * inv(b) % mod) % mod;
else ans = (ans - tmp * a % mod * inv(b) % mod) % mod;
}
cout << (ans + mod) % mod << endl;
return 0;
}
进一步优化
上面做法有两个部分的复杂度带 $\log k$,阶乘的逆元可以线性预处理,问题就只剩:如何快速预处理
1^k, 2^k, 3^k, \cdots, k^k
$$
发现 $f(t)=t^k$ 是完全积性函数,于是可以只求所有质数的 $k$ 次幂,然后线性筛。
这样整体复杂度大约是 $\mathcal{O}(k)$ 的。
参考资料