2022牛客暑期多校第三场 I. Ice Drinking
作者：互联网
2022牛客暑期多校第三场 I. Ice Drinking
题意
按随机顺序摆放 \(1,2,3,...,n\)，设随机变量 \(x\) 为数字与位置（第几个）相等的个数，给定非负整数 \(k\)，求 \(x^k\) 的期望。
分析
设错排方案为 \(P(n)\)，根据组合意义，枚举正确的个数，剩下的全部错排的方案数之和就是全排列
\[n! = \sum_{i=0}^n\binom{n}{i}P(i) \]二项式反演得到错排的容斥公式
\[P(n)=\sum_{i=0}^n(-1)^{n-i}\binom{n}{i}i! \]那么恰好有 \(j\) 个位置正确的情况数为
\[\binom{n}{j}P(n-j) =\sum_{i=0}^{n-j}(-1)^{n-i-j}\binom{n}{j}\binom{n-j}{i}i! \]于是期望为
\[\frac{1}{n!}\sum_{j=0}^n\sum_{i=0}^{n-j}(-1)^{n-i-j}\binom{n}{j}\binom{n-j}{i}i!j^k \]稍微化简一下
\[\sum_{j=0}^n\sum_{i=0}^{n-j}(-1)^{n-i-j}\frac{1}{j!(n-j-i)!}j^k \]众所周知，\(j^k\) 可以利用第二类斯特林数的组合意义，用第二类斯特林数的和式代替（组合含义为 \(j\) 个不同盒子，\(k\) 个不同小球的放法等价于枚举非空盒个数 \(t\)，在 \(j\) 个不同盒子里选择 \(t\) 个盒子，\(k\) 个不同小球划分到 \(t\) 个相同盒子里（盒子不允许为空）的方案乘上非空盒的全排列）
\[j^k=\sum_{t=0}^j \binom{j}{t}t!S(k,t) \]我们不直接代入，继续反演
\[j!S(k,j)=\sum_{t=0}^j (-1)^{j-t}\binom{j}{t}t^k \]把 \(j!\) 移到右边去，得到 \(S(k,j)\) 的表达式
\[S(k,j)=\sum_{t=0}^j (-1)^{j-t}\frac{1}{t!(j-t)!}t^k \]再次观察期望的式子，发现形式很像，还差一点。
\[\sum_{j=0}^n\sum_{i=0}^{n-j}(-1)^{n-i-j}\frac{1}{j!(n-j-i)!}j^k \]把 \(i,j\) 计算次序交换一下
\[\sum_{i=0}^n\sum_{j=0}^{n-i}(-1)^{n-i-j}\frac{1}{j!(n-j-i)!}j^k \]发现了不得了的事情，里层求和符号的形式与 \(S(k,n-i)\) 一模一样。
于是期望就变成了
\[\sum_{i=0}^n S(k,n-i)=\sum_{i=0}^n S(k,i) \]仅仅推出这个式子还不够，\(n\) 范围非常大，难以直接计算。
由 \(n,k\) 范围可以知道 \(k\) 最多比 \(n\) 大 \(5000\)，也就是如果知道一行的斯特林数之和，最多只要减掉 \(5000\) 个数。
一行斯特林数的和就是贝尔数，所以 \(n > k\) 就是第 \(k\) 个贝尔数，而贝尔数有一个很好的同余性质，可以快速计算模素数意义下的贝尔数，复杂度是 \(O(p^2 \log_p n)\) 的，而模数 \(862118861 = 857 * 997 * 1009\)，我们可以求出三个质因子为模数的答案，然后中国剩余定理解决。具体做法，参照这个链接 Bell数。
问题就变成如何逆推斯特林数的一行的后面几个，只需要 \(O((k-n)^2)\) 的递推就够了。这在《具体数学》这本书里有详细的讲解。
由于我写这道题时还不知道《具体数学》书里有相应记载，所以自行采用了另外一种方法进行递推（与书上方法不同，推出的结论看起来也不太相同，但是殊途同归，想看书上方法的请自行查阅这本书的“特殊的数”这一章节），以下是方法。
众所周知，斯特林数有一个递推式（组合意义也十分明确）
\[S(n,m)=S(n-1,m-1)+mS(n-1,m)\\ S(0,0)=1\\ S(n,0)=0(n>0) \]我们尝试找一找规律，
有
还有
\[\begin{aligned} &S(n,n-1)\\ =&S(n-1,n-2)+(n-1)\\ =&S(n-2,n-3)+(n-2)+(n-1)\\ ...\\ =&S(1,0)+1+2+...+(n-1)\\ =&\binom{n}{2} \end{aligned} \]还有
\[\begin{aligned} &S(n,n-2)\\ =&S(n-1,n-3)+(n-2)\binom{n-1}{2}\\ =&S(n-2,n-4)+(n-3)\binom{n-2}{2}+(n-2)\binom{n-1}{2}\\ ...\\ =&S(2,0)+1\cdot\binom{2}{2}+2\cdot\binom{3}{2}+...+(n-2)\binom{n-1}{2}\\ =&\sum_{i=1}^{n-2}i\binom{i+1}{2}\\ =&\sum_{i=1}^{n-2}\binom{i+1}{2}+\sum_{i=1}^{n-2}(i-1)\binom{i+1}{2}\\ =&\sum_{i=1}^{n-2}\binom{i+1}{2}+\frac{3!}{2!}\sum_{i=1}^{n-2}\binom{i+1}{3}\\ =&\binom{n}{3}+3\binom{n}{4} \end{aligned} \]还有
\[\begin{aligned} &S(n,n-3)\\ =&S(n-1,n-4)+(n-3)\left(\binom{n-1}{3}+3\binom{n-1}{4}\right)\\ =&S(n-2,n-5)+(n-4)\left(\binom{n-2}{3}+3\binom{n-2}{4}\right)\\ ...\\ =&\sum_{i=1}^{n-3}i\left(\binom{i+2}{3}+3\binom{i+2}{4}\right)\\ =&\sum_{i=1}^{n-3}i\binom{i+2}{3}+3\sum_{i=1}^{n-3}i\binom{i+2}{4}\\ =&\sum_{i=1}^{n-3}(i-1)\binom{i+2}{3}+\sum_{i=1}^{n-3}\binom{i+2}{3}+3\sum_{i=1}^{n-3}(i-2)\binom{i+2}{4}+6\sum_{i=1}^{n-3}\binom{i+2}{4}\\ =&\frac{4!}{3!}\sum_{i=1}^{n-3}\binom{i+2}{4}+\sum_{i=1}^{n-3}\binom{i+2}{3}+3\cdot\frac{5!}{4!}\sum_{i=1}^{n-3}\binom{i+2}{5}+6\sum_{i=1}^{n-3}\binom{i+2}{4}\\ =&4\binom{n}{5}+\binom{n}{4}+15\binom{n}{6}+6\binom{n}{5}\\ =&\binom{n}{4}+10\binom{n}{5}+15\binom{n}{6} \end{aligned} \]推了 \(3\) 个，发现这些都是若干组合数的和（先不管系数具体多少），记系数为 \(p_{i,j}\) （角标是为了后面推这个系数的递推式用的），具体地归纳一下，有
\[S(n,n-k)=\sum_{i=1}^kp_{k,i}\binom{n}{k+i} \]怎么求系数 \(p_{i,j}\) 呢？利用数学归纳法，就能求递推式了。
假定如下式子对于任意 \(n\) 成立（\(k\) 固定）
\[S(n,n-k+1)=\sum_{i=1}^{k-1}p_{k-1,i}\binom{n}{k-1+i} \]我们开始由上面的式子推 \(S(n,n-k)\)，方法和上面找规律的推法很像。
\[\begin{aligned} &S(n,n-k)\\ =&S(n-1,n-k-1)+(n-k)S(n-1,n-k)\\ =&S(n-2,n-k-2)+(n-k-1)S(n-2,n-k-1)+(n-k)S(n-1,n-k)\\ ...\\ =&\sum_{i=1}^{n-k}iS(i+k-1,i)\\ =&\sum_{i=1}^{n-k}\sum_{j=1}^{k-1}i\cdot p_{k-1,j}\binom{i+k-1}{k-1+j}\\ =&\sum_{j=1}^{k-1}\sum_{i=1}^{n-k}i\cdot p_{k-1,j}\binom{i+k-1}{k-1+j}\\ =&\sum_{j=1}^{k-1}\sum_{i=1}^{n-k}(i-j)\cdot p_{k-1,j}\binom{i+k-1}{k-1+j}+\sum_{j=1}^{k-1}\sum_{i=1}^{n-k}j\cdot p_{k-1,j}\binom{i+k-1}{k-1+j}\\ =&\sum_{j=1}^{k-1}(k+j)\cdot p_{k-1,j}\sum_{i=1}^{n-k}\binom{i+k-1}{k+j}+\sum_{j=1}^{k-1}j\cdot p_{k-1,j}\sum_{i=1}^{n-k}\binom{i+k-1}{k-1+j}\\ =&\sum_{j=1}^{k-1}(k+j)\cdot p_{k-1,j}\binom{n}{k+j+1}+\sum_{j=1}^{k-1}j\cdot p_{k-1,j}\binom{n}{k+j}\\ =&\sum_{j=2}^{k}(k+j-1)\cdot p_{k-1,j-1}\binom{n}{k+j}+\sum_{j=1}^{k-1}j\cdot p_{k-1,j}\binom{n}{k+j}\\ \end{aligned} \]为了优雅（方便）起见，定义 \(p_{i,0}=p_{i,i+1}=0\)
\[\begin{aligned} &S(n,n-k)\\ =&\sum_{j=1}^{k}(k+j-1)\cdot p_{k-1,j-1}\binom{n}{k+j}+\sum_{j=1}^{k}j\cdot p_{k-1,j}\binom{n}{k+j}\\ =&\sum_{j=1}^{k}\left(j\cdot p_{k-1,j}+(k+j-1)\cdot p_{k-1,j-1}\right)\binom{n}{k+j}\\ \end{aligned} \]和我们猜想的结论对比一下
\[S(n,n-k)=\sum_{i=1}^kp_{k,i}\binom{n}{k+i} \]递推式就得到了（讲真，这个式子很优美）
\[p_{k,i}=i\cdot p_{k-1,i}+(k+i-1)\cdot p_{k-1,i-1} \]上面的数归也就非常完美的证明的这个猜想。
因为 \(S(n,n-1)=\binom{n}{2}\)，所以 \(p_{1,1}=1\)，加之定义的 \(p_{i,0}=p_{i,i+1}=0\)，所有的递推边界条件就齐了。
用卢卡斯定理+中国剩余定理计算组合数取模，这道题就结束了（非常毒瘤）。
代码
#include <iostream>
#include <vector>
using namespace std;
typedef long long Lint;
constexpr Lint fpow(Lint a, Lint b, Lint mod) {
a %= mod;
Lint res = 1;
for (; b; b >>= 1) {
if (b & 1)
res = res * a % mod;
a = a * a % mod;
}
return res;
}
constexpr Lint mod = 857 * 997 * 1009;
constexpr Lint p1 = 857;
constexpr Lint p2 = 997;
constexpr Lint p3 = 1009;
constexpr Lint mod_p1 = mod / p1;
constexpr Lint mod_p2 = mod / p2;
constexpr Lint mod_p3 = mod / p3;
constexpr Lint inv_mod_p1 = fpow(mod_p1, p1 - 2, p1);
constexpr Lint inv_mod_p2 = fpow(mod_p2, p2 - 2, p2);
constexpr Lint inv_mod_p3 = fpow(mod_p3, p3 - 2, p3);
constexpr Lint c1 = mod_p1 * inv_mod_p1 % mod;
constexpr Lint c2 = mod_p2 * inv_mod_p2 % mod;
constexpr Lint c3 = mod_p3 * inv_mod_p3 % mod;
Lint CRT(Lint a1, Lint a2, Lint a3) {
return (a1 * c1 % mod + a2 * c2 % mod + a3 * c3 % mod) % mod;
}
Lint fac_p1[p1], fac_p2[p2], fac_p3[p3];
Lint inv_fac_p1[p1], inv_fac_p2[p2], inv_fac_p3[p3];
Lint bell_p1[p1 + 1], bell_p2[p2 + 1], bell_p3[p3 + 1];
Lint temp[2][2000];
Lint C_p(Lint n, Lint k, Lint p) {
if (n < k || n < 0 || k < 0)
return 0;
if (p == p1) {
return fac_p1[n] * inv_fac_p1[k] % p * inv_fac_p1[n - k] % p;
} else if (p == p2) {
return fac_p2[n] * inv_fac_p2[k] % p * inv_fac_p2[n - k] % p;
} else {
return fac_p3[n] * inv_fac_p3[k] % p * inv_fac_p3[n - k] % p;
}
}
Lint Lucas(Lint n, Lint k, Lint p) {
if (n < k || n < 0 || k < 0)
return 0;
return n == 0 ? 1 : Lucas(n / p, k / p, p) * C_p(n % p, k % p, p) % p;
}
Lint C(Lint n, Lint k) {
if (n < k || n < 0 || k < 0)
return 0;
return CRT(Lucas(n, k, p1), Lucas(n, k, p2), Lucas(n, k, p3));
}
void init() {
fac_p1[0] = fac_p2[0] = fac_p3[0] = 1;
for (int i = 1; i < p1; i++) {
fac_p1[i] = fac_p1[i - 1] * i % p1;
}
for (int i = 1; i < p2; i++) {
fac_p2[i] = fac_p2[i - 1] * i % p2;
}
for (int i = 1; i < p3; i++) {
fac_p3[i] = fac_p3[i - 1] * i % p3;
}
inv_fac_p1[p1 - 1] = fpow(fac_p1[p1 - 1], p1 - 2, p1);
inv_fac_p2[p2 - 1] = fpow(fac_p2[p2 - 1], p2 - 2, p2);
inv_fac_p3[p3 - 1] = fpow(fac_p3[p3 - 1], p3 - 2, p3);
for (int i = p1 - 2; i >= 0; i--) {
inv_fac_p1[i] = inv_fac_p1[i + 1] * (i + 1) % mod;
}
for (int i = p2 - 2; i >= 0; i--) {
inv_fac_p2[i] = inv_fac_p2[i + 1] * (i + 1) % mod;
}
for (int i = p3 - 2; i >= 0; i--) {
inv_fac_p3[i] = inv_fac_p3[i + 1] * (i + 1) % mod;
}
bell_p1[0] = bell_p2[0] = bell_p3[0] = 1;
temp[0][0] = 1;
for (int i = 1; i <= p1; i++) {
int u = i & 1;
bell_p1[i] = temp[u][0] = temp[u ^ 1][i - 1];
for (int j = 1; j <= i; j++) {
temp[u][j] = (temp[u][j - 1] + temp[u ^ 1][j - 1]) % p1;
}
}
temp[0][0] = 1;
for (int i = 1; i <= p2; i++) {
int u = i & 1;
bell_p2[i] = temp[u][0] = temp[u ^ 1][i - 1];
for (int j = 1; j <= i; j++) {
temp[u][j] = (temp[u][j - 1] + temp[u ^ 1][j - 1]) % p2;
}
}
temp[0][0] = 1;
for (int i = 1; i <= p3; i++) {
int u = i & 1;
bell_p3[i] = temp[u][0] = temp[u ^ 1][i - 1];
for (int j = 1; j <= i; j++) {
temp[u][j] = (temp[u][j - 1] + temp[u ^ 1][j - 1]) % p3;
}
}
}
Lint get_Bell_p(Lint n, Lint p) {
if (p == p1)
return bell_p1[n];
else if (p == p2)
return bell_p2[n];
else
return bell_p3[n];
}
Lint Bell_p(Lint n, Lint p) {
vector<Lint> a;
for (; n; n /= p) {
a.push_back(n % p);
}
for (int i = 0; i <= p; i++) {
temp[0][i] = get_Bell_p(i, p);
}
for (int i = 1; i < a.size(); i++) {
for (int j = 1; j <= a[i]; j++) {
for (int k = 0; k < p; k++) {
temp[1][k] = (i * temp[0][k] % p + temp[0][k + 1]) % p;
}
temp[1][p] = (temp[1][0] + temp[1][1]) % p;
for (int k = 0; k <= p; k++)
temp[0][k] = temp[1][k];
}
}
return temp[0][a[0]];
}
Lint Bell(Lint n) {
return CRT(Bell_p(n, p1), Bell_p(n, p2), Bell_p(n, p3));
}
Lint p[5005][5005];
int main() {
init();
Lint n, k;
cin >> n >> k;
if (n >= k) {
cout << Bell(k) << '\n';
return 0;
}
Lint sum = Bell(k);
p[1][1] = 1;
for (int i = 2; i <= k - n; i++) {
for (int j = 1; j <= i; j++) {
p[i][j] = (i + j - 1) * p[i - 1][j - 1] % mod + j * p[i - 1][j] % mod;
if (p[i][j] >= mod)
p[i][j] -= mod;
}
}
Lint rm = 1;
for (int i = 1; i < k - n; i++) {
for (int j = 1; j <= i; j++) {
rm += p[i][j] * C(k, i + j) % mod;
if (rm >= mod)
rm -= mod;
}
}
sum = (sum - rm + mod) % mod;
cout << sum << '\n';
return 0;
}
标签：p1,sum,Lint,多校,Ice,fac,binom,第三场,mod 来源： https://www.cnblogs.com/Bamboo-Wind/p/16523152.html