题解 | #Ice Drinking#

【分析】

先推一波公式：

答案 $r e s$ 显然有公式：（其中 $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}$

现在的问题化为，如何求解出 $B e l l_{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;
}