[LOJ 6055]「from CommonAnts」一道数学题 加强版（神仙拉格朗日插值） | 错题本

题目

「from CommonAnts」一道数学题 加强版

分析

记 $F_k(n)=f(k,n)$，统计每个 $x^k$ 的贡献可以得到 $$F_k(n)=\sum _{i=0}^{n-1}\left(2^{n-i-1}\cdot i^k\right)+n^k$$ 具体来说可以画个二叉树理解：树上深度为 $i$ 的结点的权值为 $(n-i+1{)}^k$，根以及他的右子树权值和为 $F_k(n)$，根的左儿子和根的左儿子的右子树权值和 $F_k(n-1)$，以此类推。

显然不能直接插，因为 $n$ 在指数上。化简一下：
$$\begin{array}{rl}{F}_k(n)& =\sum _{i=0}^{n-1}\left(2^{n-i-1}\cdot i^k\right)+n^k\\ & =2^{n-1}\sum _{i=0}^{n-1}\frac{i^k}{2^i}+n^k\end{array}$$ 令 $a=\frac{1}{2}$，则 $$F_k(n)=2^{n-1}\sum _{i=0}^{n-1}{a}^i{i}^k+n^k$$ 指数转下降幂可以方便错位相减找递推式，进而找到指数函数的多项式形式：$$\begin{array}{rl}{F}_k(n)& =2^{n-1}\sum _{i=0}^{n-1}\left(a^i\sum _{j=0}^k\left\{\genfrac{}{}{0ex}{}{k}{j}\right\}{i}^{\underline{j}}\right)+n^k\\ & =2^{n-1}\sum _{j=0}^k\left\{\genfrac{}{}{0ex}{}{k}{j}\right\}\sum _{i=0}^{n-1}{a}^i{i}^{\underline{j}}+n^k\end{array}$$ 令 $S_k(n)={\sum }_{i=0}^{n-1}{a}^i{i}^{\underline{k}}$，则 $$F_k(n)=2^{n-1}\sum _{j=0}^k\left\{\genfrac{}{}{0ex}{}{k}{j}\right\}{S}_j(n)+n^k$$

指数转下降幂：
$\left\{\genfrac{}{}{0ex}{}{n}{m}\right\}$ 表示第二类斯特林数，即“$n$ 个不同的球，放入 $m$ 个相同的盒子，盒子不可以为空的方案数”。有公式 $$n^k=\sum _{i=0}^k\left\{\genfrac{}{}{0ex}{}{k}{i}\right\}\cdot i!\cdot C_n^i$$ 等式左边是“$k$ 个不同的球放入 $n$ 个不同的盒子，盒子可以为空的方案数”，右边即为枚举不为空的盒子数，分别计算，因此两边相等。
把组合数拆开可以得到：$$n^k=\sum _{i=0}^k\left\{\genfrac{}{}{0ex}{}{k}{i}\right\}\cdot n^{\underline{i}}$$

考虑求 $S_k(n)$ 的递推式，用错位相减法：
$$\begin{array}{rl}{S}_k(n)& =\sum _{i=0}^{n-1}{a}^i{i}^{\underline{k}}\\ a{S}_k(n)& =\sum _{i=1}^n{a}^i(i-1{)}^{\underline{k}}\\ \Rightarrow (1-a)S_k(n)& =\sum _{i=1}^{n-1}{a}^i\left[i^{\underline{k}}-(i-1{)}^{\underline{k}}\right]+a^0{0}^{\underline{k}}-a^n(n-1{)}^{\underline{k}}\\ & =\sum _{i=1}^n{a}^i\left[i^{\underline{k}}-(i-1{)}^{\underline{k}}\right]-a^n{n}^{\underline{k}}+a^n(n-1{)}^k+a^0{0}^{\underline{k}}-a^n(n-1{)}^k\\ & =\sum _{i=1}^n{a}^i\left[i^{\underline{k}}-(i-1{)}^{\underline{k}}\right]-a^n{n}^{\underline{k}}+a^0{0}^{\underline{k}}\\ & =\sum _{i=1}^n{a}^i(i-1{)}^{\underline{k-1}}[i-(i-k)]-a^n{n}^{\underline{k}}+a^0{0}^{\underline{k}}\\ & =ak\sum _{i=0}^{n-1}{a}^i{i}^{\underline{k-1}}-a^n{n}^{\underline{k}}+a^0{0}^{\underline{k}}\\ & =ak{S}_{k-1}(n)-a^n{n}^{\underline{k}}+a^0{0}^{\underline{k}}\\ \Rightarrow S_k(n)& =\frac{ak}{1-a}{S}_{k-1}(n)-\frac{a^n{n}^{\underline{k}}+a^0{0}^{\underline{k}}}{1-a}\end{array}$$ 特别地，$S_0(n)={\sum }_{i=0}^{n-1}{a}^i=\frac{1-a^n}{1-a}$。将 $a$ 替换为 $\frac{1}{2}$ 得到 $$S_k(n)=\{\begin{array}{ll}2\left(1-\frac{1}{2^n}\right)& k=1\\ k{S}_{k-1}(n)-\frac{n^{\underline{k}}}{2^{n-1}}& k>1\end{array}$$ 这时候手刨几项：

• $S_0(n)=2(1-\frac{1}{2^n})$
• $S_1(n)=1\cdot 2(1-\frac{1}{2^n})-\frac{n}{2^{n-1}}$
• $S_2(n)=2\left[1\cdot 2(1-\frac{1}{2^n})-\frac{n}{2^{n-1}}\right]-\frac{n(n-1)}{2^n}$
• $S_3(n)=3\left\{2\left[1\cdot 2(1-\frac{1}{2^n})-\frac{n}{2^{n-1}}\right]-\frac{n(n-1)}{2^n}\right\}-\frac{n(n-1)(n-2)}{2^n}$

接下来就是神仙操作了，还记得常系数齐次线性递推么？ 将常数项分成两部分，将含 $n$ 的部分提出一个 $\frac{1}{2^n}$：$$\begin{array}{rl}{S}_3(n)& =3\left\{2\left[1\cdot 2(1-\frac{1}{2^n})-\frac{n}{2^{n-1}}\right]-\frac{n(n-1)}{2^n}\right\}-\frac{n(n-1)(n-2)}{2^n}\\ & =12-\frac{1}{2^n}[12+6n+3n(n-1)+n(n-1)(n-2)]\end{array}$$ 两个 $12$ 是一样的！其实原因就在于 $S_0(n)=2(1-\frac{1}{2^n})$，两个 $1$ 可以分给两边，使得新的一个 $S_k(n)$ 可以写成 $$S_k(n)=G_k(0)-\frac{G_k(n)}{2^n}$$ 其中 $G_k(n)$ 显然是一个 $k$ 次多项式，而 $G_k(0)$ 正式它的常数项！

可能的疑问是，为什么不能直接提出一个 $\frac{1}{2^n}$，变成 $$S_k(n)=\frac{T_k(n)}{2^n}$$ 原因是，$S_0(n)$ 事实上是两部分，如果不分开的话，$T_k(n)$ 中会有 $2^n$ 这样的东西（$S_0(n)$ 中的 “$1$”造成的），所以不可以直接提 $\frac{1}{2^n}$。分开过后，$G_k(n)$ 就不含指数函数了！

回代进 $F_k(n)$：$$\begin{array}{rl}{F}_k(n)& =2^{n-1}\sum _{j=0}^k\left\{\genfrac{}{}{0ex}{}{k}{j}\right\}{S}_j(n)+n^k\\ & =2^{n-1}\sum _{j=0}^k\left\{\genfrac{}{}{0ex}{}{k}{j}\right\}\left(G_j(0)-\frac{G_j(n)}{2^n}\right)+n^k\\ & =2^{n-1}\left(\sum _{i=0}^k\left\{\genfrac{}{}{0ex}{}{k}{i}\right\}{G}_i(0)-\frac{{\sum }_{i=0}^k\left\{\genfrac{}{}{0ex}{}{k}{i}\right\}{G}_i(n)}{2^n}\right)+n^k\end{array}$$ 令 $$H_k(n)=\sum _{i=0}^k\left\{\genfrac{}{}{0ex}{}{k}{i}\right\}{G}_i(n)$$ 显然 $H_i(n)$ 为 $k$ 次多项式，且 $$F_k(n)=2^{n-1}\left(H_k(0)-\frac{H_k(n)}{2^n}\right)+n^k$$ 接下来求出 $R_k(n)=H_k(0)-\frac{1}{2^n}{H}_k(n)$ 即可。

• $H_k(0)$ 的部分：
  这个函数我们一项都算不出来，回顾一下开头，发现 $R_k(n)$ 可以算：$$R_k(n)=\sum _{i=0}^{n-1}\frac{i^k}{2^i}$$ 把式子弄一下 $$H_k(n)=2^n{H}_k(0)-2^n{R}_k(n)$$ 那么用一个二元组 { x , y } \{x, y\} {x,y} 表示 $H_k(n)$ 即 $H_k(n)=x{H}_k(0)-y$，这个二元组显然可以快速地进行加、减和乘一个常的运算。有了这些我们就可以用 $H_k(0),H_k(1),\cdots ,H_k(k)$ 的二元表示插出 $H_k(k+1)$ 的二元表示记为 { a 1 , b 1 } \{a_1,b_1\} {a1​,b1​}。 $H_k(n)=2^n{H}_k(0)-2^n{R}_k(n)$ 也能得到一个 $H_k(k+1)$ 的二元表示记为 { a 2 , b 2 } \{a_2, b_2\} {a2​,b2​}。由 $H_k(k+1)=H_k(k+1)$ 得到 $$a_1{H}_k(0)-b_1=a_2{H}_k(0)-b_2$$ 于是 $$H_k(0)=\frac{b_1-b_2}{a_1-a_2}$$

  注意这两个二元表示是不等价的！可以理解为只要给出次数，拉格朗日插值就可以形式化地创造一个表示方法，这个新的值的计算方式与你原来计算点值的方式没有任何关系，因此这个方程有解。

• $H_k(n)$ 的部分：
  得到 $H_k(0)$，就知道了 $H_k(1),H_k(2),\cdots ,H_k(k)$，直接拉格朗日插值即可。

---

还有一个细节是算 $2^n$ 什么的，$n$ 有一百万位显然不能快速幂……要用欧拉定理：$$a^{\phi (p)}\equiv 1(\text{mod}p)$$ 这里 $p=10^9+7$ 是个质数，所以 $\phi (p)=p-1$，所以 $2^n\equiv 2^{n\text{mod}(p-1)}(\text{mod}p)$。

---

做完了。这题太夸张了。据说 高阶差分 可以简单一点点，或者最开始就往 常系数线性递推 化。太菜了看不怎么懂……

代码

变量名有出入，但思路是完全一致的。

#include <algorithm>
#include <cstdio>
#include <cstring>

typedef std::pair<int, int> PII;

const int MAXN = 1000000;
const int MOD = 1000000007;
const int PHI = MOD - 1;

inline int Add(int x, const int &y) {
	x += y; if (x >= MOD) x -= MOD; return x;
}

inline int Sub(int x, const int &y) {
	x -= y; if (x < 0) x += MOD; return x;
}

inline int Mul(const int &x, const int &y) {
	return (long long)x * y % MOD;
}

int Pow(int x, int y) {
	int ret = 1;
	while (y) {
		if (y & 1)
			ret = Mul(ret, x);
		y >>= 1;
		x = Mul(x, x);
	}
	return ret;
}

const int INV2 = Pow(2, MOD - 2);

PII operator + (const PII &x, const PII &y) {
	return std::make_pair(Add(x.first, y.first), Add(x.second, y.second));
}

PII operator - (const PII &x, const PII &y) {
	return std::make_pair(Sub(x.first, y.first), Sub(x.second, y.second));
}

PII operator * (const PII &x, const int &y) {
	return std::make_pair(Mul(x.first, y), Mul(x.second, y));
}

int N, P, K;
int Fac[MAXN + 5], InvFac[MAXN + 5];

int H[MAXN + 5];
int F[MAXN + 5];
PII G[MAXN + 5];

int Pre[MAXN + 5], Suf[MAXN + 5];

void Init(int n) {
	Pre[0] = n;
	for (int i = 1; i <= K; i++)
		Pre[i] = Mul(Pre[i - 1], Sub(n, i));
	Suf[K + 1] = 1;
	for (int i = K; i >= 0; i--)
		Suf[i] = Mul(Suf[i + 1], Sub(n, i));
}

PII Lagrange1(int n) {
	Init(n);
	PII res(0, 0);
	for (int i = 0; i <= K; i++) {
		int A = Mul(i ? Pre[i - 1] : 1, Suf[i + 1]);
		int B = Mul(InvFac[i], InvFac[K - i]);
		PII C = G[i] * Mul(A, B);
		res = ((K - i) & 1) ? (res - C) : (res + C);
//		printf("%d %d (%d, %d)\n", A, B, res.first, res.second);
	}
	return res;
}

int Lagrange2(int n) {
	Init(n);
	int res = 0;
	for (int i = 0; i <= K; i++) {
		int A = Mul(i ? Pre[i - 1] : 1, Suf[i + 1]);
		int B = Mul(InvFac[i], InvFac[K - i]);
		int C = Mul(H[i], Mul(A, B));
		res = ((K - i) & 1) ? Sub(res, C) : Add(res, C);
	}
	return res;
}

int main() {
	char c = getchar();
	while (c < '0' || c > '9')
		c = getchar();
	while (c >= '0' && c <= '9') {
		P = ((long long)P * 10 + c - '0') % PHI;
		N = ((long long)N * 10 + c - '0') % MOD;
		c = getchar();
	}
	scanf("%d", &K);

	Fac[0] = 1;
	for (int i = 1; i <= K; i++)
		Fac[i] = Mul(Fac[i - 1], i);
	InvFac[K] = Pow(Fac[K], MOD - 2);
	for (int i = K - 1; i >= 0; i--)
		InvFac[i] = Mul(InvFac[i + 1], i + 1);

	for (int i = 1; i <= K + 1; i++)
		F[i] = Add(F[i - 1], Mul(Pow(INV2, i - 1), Pow(i - 1, K)));
	for (int i = 0; i <= K + 1; i++)
		G[i] = std::make_pair(Pow(2, i), Mul(Pow(2, i), F[i]));

	PII x = G[K + 1], y = Lagrange1(K + 1);
	int A = Sub(y.second, x.second);
	int B = Sub(x.first, y.first);
	H[0] = Mul(A, Pow(B, MOD - 2));
	for (int i = 1; i <= K; i++)
		H[i] = Add(Mul(G[i].first, H[0]), G[i].second);

	int Ans = Sub(Mul(Lagrange2(N), Pow(INV2, P)), H[0]);
	Ans = Add(Mul(Pow(2, P ? (P - 1) : (MOD - 2)), Ans), Pow(N, K));
	printf("%d", Ans);
	return 0;
}