题解 I. Ice Drinking "蔚来杯"2022牛客暑期多校训练营3
【分析】
先推一波公式：
答案 \(res\) 显然有公式：（其中 \(D_n\) 表示 \(n\) 个元素全部错排的方案数）
\(\begin{aligned}res&={1\over n!}\sum_{x=0}^n\dbinom n x x^kD_{n-x}\\&=\sum_{x=0}^n{x^k\over x!}\cdot {D_{n-x}\over (n-x)!}\end{aligned}\)
由于错排问题有公式：\(\displaystyle D_n=\lfloor{n!\over e}\rfloor=n!\cdot \sum_{i=0}^n{(-1)^i\over i!}\)
因此答案为：\(\displaystyle res=\sum_{x=0}^n {x^k\over x!}\cdot \sum_{i=0}^{n-x} {(-1)^i\over i!}\)
对 \(x^k\) 用第二类斯特林数展开为下降幂，得到：
\(\begin{aligned}res&=\sum_{x=0}^n {x^k\over x!}\cdot \sum_{i=0}^{n-x} {(-1)^i\over i!}\\&=\sum_{t=0}^k\begin{Bmatrix}k\\t\end{Bmatrix}\sum_{x=0}^n{x^{\underline t}\over x!}\cdot \sum_{i=0}^{n-x}{(-1)^i\over i!}\\&=\sum_{t=0}^{\min(n, k)}\begin{Bmatrix}k\\t\end{Bmatrix}\sum_{x=0}^{n-t}{1\over x!}\cdot \sum_{i=0}^{(n-t)-x}{(-1)^i\over i!}\\&=\sum_{t=0}^{\min(n, k)}\begin{Bmatrix}k\\t\end{Bmatrix}\sum_{T=0}^{n-t}\sum_{i+j=T}{1\over j!}\cdot {(-1)^i\over i!}\\&=\sum_{t=0}^{\min(n, k)}\begin{Bmatrix}k\\t\end{Bmatrix}[1+\sum_{T=1}^{n-t}{1\over T!}\sum_{i+j=T}\dbinom T {i, j}(-1)^i]\\&=\sum_{t=0}^{\min(n, k)}\begin{Bmatrix}k\\t\end{Bmatrix}[1+\sum_{T=1}^{n-t}{1\over T!}\cdot (1-1)^T]\\&=\sum_{t=0}^{\min(n, k)}\begin{Bmatrix}k\\t\end{Bmatrix}\\&=Bell_k-\sum_{t=n+1}^k\begin{Bmatrix}k\\t\end{Bmatrix}\end{aligned}\)
现在的问题化为，如何求解出 \(Bell_k\) 和 \(\begin{Bmatrix}k\\t\end{Bmatrix}\)
关于前者，根据题解提到的 Touchard 同余，可以得到贝尔数在模质数 \(p\) 意义下，有 \(Bell_n\equiv Bell_{n-p}+Bell_{n-(p-1)}\pmod {p}\)
这显然是一个 \(p\) 阶线性递推的形式，可以在 \(O(p^2\log k)\) 的时间内求解出 \(Bell_k\bmod p\) 的答案
而关于第二类斯特林数，有公式：
\(\displaystyle \begin{Bmatrix}x\\x-n\end{Bmatrix}=\sum_{k=0}^{n-1}\left<\left<\begin{matrix}n\\k\end{matrix}\right>\right>\dbinom {x+n-k-1} {2n}\)
其中，\(\displaystyle \left<\left<\begin{matrix}n\\k\end{matrix}\right>\right>\) 为第二类欧拉数，或二阶欧拉数。它表示多重集 \(\{1, 1, 2, 2, \cdots, n, n\}\) 构成的序列中，任意两个相同的 \(i\) 之间，所有数字都大于 \(i\) 的序列方案数。如 \(123321\) 是满足要求的，但 \(123123\) 不是，因为两个 \(2\) 之间含有小于 \(2\) 的 \(1\) ；两个 \(3\) 之间含有小于 \(3\) 的 \(1,2\)
第二类欧拉数有递推公式： \(\displaystyle \left<\left<\begin{matrix}n\\k\end{matrix}\right>\right>=(2n-k-1)\left<\left<\begin{matrix}n-1\\k-1\end{matrix}\right>\right>+(k+1)\left<\left<\begin{matrix}n-1\\k\end{matrix}\right>\right>\)
因此，模数为小质数 \(p\) 时，我们可以采用 \(O(p^2\log k)\) 的时间求解 \(Bell_k\bmod p\) ；\(O(p^2)\) 递推 \(\dbinom n m (0\leq n,m<p)\) 的数值，然后再在 \(O((k-n)^2)\) 的时间内地推出所有 \(O(k-n)\) 的第二类欧拉数，最后对于每个 \(\begin{Bmatrix}k\\t\end{Bmatrix}\) 利用第二类欧拉数和 Lucas 定理进行 \(O((k-t)\log k)\) 求解
总复杂度为 \(\displaystyle O(p^2\log k)+O(p^2)+O((k-n)^2)+\sum_{t=n+1}^k O((k-t)\log k)=O(p^2\log k+(k-n)^2\log k)\)
然而，我们的模数 \(862118861=857\times 997\times 1009\) 是个合数，不适用上面的算法。这当然很好解决，我们单独跑每个质数的结果，最后用中国剩余定理合并即可
【代码】
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef vector<int> vi;
#define de(x) cout << #x <<" = "<<(x)<<endl
#define dd(x) cout << #x <<" = "<<(x)<<" "
#define all(a) a.begin(), a.end()
#define rep(i, a, b) for(int i=(a); i<(b); ++i)
#define per(i, a, b) for(int i=(b)-1; i>=(a); --i)
const int P=862118861, p[]={0, 857, 997, 1009};
ll n, k, S[1024][1024], Bell[1024];
inline ll kpow(ll a,ll x, ll p) { ll ans=1; for(;x;x>>=1, a=a*a%p) if(x&1) ans=ans*a%p; return ans; }
inline ll inv(ll a, ll p) { return kpow(a, p-2, p); }
int linear_recurrence(ll n, int m, vi a, vi c, int P) {
if (n<m) return (a[n]+P)%P;
vector<ll> v(m, 0), u(m<<1, 0);
v[0] = 1;
for(ll x = 0, W = n ? 1ll<<(63 - __builtin_clzll(n)) : 0; W; W >>= 1, x <<= 1) {
fill(all(u), 0);
int b = !!(n & W); if(b) x++;
if(x < m) u[x] = 1;
else {
rep(i, 0, m) rep(j, 0, m) (u[i + b + j] += v[i] * v[j]) %= P;
per(i, m, 2*m) rep(j, 0, m) (u[i - m + j] += c[j] * u[i]) %= P;
}
copy(u.begin(), u.begin() + m, v.begin());
}
ll ans = 0;
rep(i, 0, m) (ans += v[i] * a[i]) %= P;
return (ans+P)%P;
}
vi linear_a, linear_c;
inline int work_bell(ll k, int p) {
linear_a.clear();
linear_c.clear();
S[0][0]=1;
Bell[0]=1;
linear_a.push_back(Bell[0]);
for(int i=1; i<=p; ++i) {
for(int j=1; j<=i; ++j)
S[i][j]=(S[i-1][j-1]+(ll)S[i-1][j]*j)%p;
Bell[i]=0;
for(int j=1; j<=i; ++j)
Bell[i]=(Bell[i]+S[i][j])%p;
linear_a.push_back(Bell[i]);
}
linear_c.resize(p);
linear_c[0]=linear_c[1]=1;
return linear_recurrence(k, p, linear_a, linear_c, p);
}
int Euler2[5010][5010], C[1024][1024];
inline void pre_Euler2(int p) {
Euler2[0][0]=1;
for(int n=1; n<=5000; ++n) {
Euler2[n][0]=1;
for(int k=1; k<=n; ++k)
Euler2[n][k]=((2ll*n-k-1)*Euler2[n-1][k-1]+(k+1ll)*Euler2[n-1][k])%p;
}
}
inline void pre_C(int p) {
C[0][0]=1;
for(int i=1; i<p; ++i) {
C[i][0]=C[i][i]=1;
for(int j=1; j<i; ++j)
C[i][j]=(C[i-1][j]+C[i-1][j-1])%p;
}
}
inline int Lucas(ll n, ll m, int p) {
int res=1;
for(; n||m; n/=p, m/=p)
res=(ll)res*C[n%p][m%p]%p;
return res;
}
inline int S2(ll x, int n, int p) {
int res=0;
for(int k=0; k<=n; ++k)
res=(res+(ll)Euler2[n][k]*Lucas(x+n-k-1, n+n, p))%p;
return res;
}
inline int work(int p) {
int res=work_bell(k, p);
pre_Euler2(p);
pre_C(p);
for(ll t=n+1; t<=k; ++t)
res=(res-S2(k, k-t, p))%p;
return res<0?res+p:res;
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0); cout.tie(0);
cin>>n>>k;
int x1=work(p[1]);
int x2=work(p[2]);
int x3=work(p[3]);
ll m1=P/p[1], m2=P/p[2], m3=P/p[3];
ll res=x1*m1*inv(m1, p[1])+x2*m2*inv(m2, p[2])+x3*m3*inv(m3, p[3]);
cout<<res%P;
cout.flush();
return 0;
}
