求自然数幂的和 – CF622F

问题

给定 $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)$ 的。

参考资料

地靈殿 – 洛谷博客

HenryHuang 的博客 – 洛谷博客

拉格朗日插值 – OI Wiki

 

暂无评论

发送评论 编辑评论

|´・ω・)ノ
ヾ(≧∇≦*)ゝ
(☆ω☆)
(╯‵□′)╯︵┴─┴
 ̄﹃ ̄
(/ω\)
∠( ᐛ 」∠)_
(๑•̀ㅁ•́ฅ)
→_→
୧(๑•̀⌄•́๑)૭
٩(ˊᗜˋ*)و
(ノ°ο°)ノ
(´இ皿இ`)
⌇●﹏●⌇
(ฅ´ω`ฅ)
(╯°A°)╯︵○○○
φ( ̄∇ ̄o)
ヾ(´・ ・`。)ノ"
( ง ᵒ̌皿ᵒ̌)ง⁼³₌₃
(ó﹏ò。)
Σ(っ °Д °;)っ
( ,,´・ω・)ノ"(´っω・`。)
╮(╯▽╰)╭
o(*////▽////*)q
>﹏<
( ๑´•ω•) "(ㆆᴗㆆ)
😂
😀
😅
😊
🙂
🙃
😌
😍
😘
😜
😝
😏
😒
🙄
😳
😡
😔
😫
😱
😭
💩
👻
🙌
🖕
👍
👫
👬
👭
🌚
🌝
🙈
💊
😶
🙏
🍦
🍉
😣
Source: github.com/k4yt3x/flowerhd
颜文字
Emoji
小恐龙
花!
上一篇
下一篇
window.pjaxLoaded = function(){ } window.pjaxLoaded();