题解 C. Lucky Sequence 2020 银川 Regional

传送门

之前现场没来得及开这一题，现在挂 cf 上就补了

第二类斯特林数 + 贝尔数 + 多项式 exp + 线段树分治 FFT + cdq 分治 FFT

感觉出出来就没打算让选手在现场过了...

---

【大意】

询问有多少个长度为 \(n\) 的序列 \(a_1, a_2, \cdots, a_n\) 满足如下条件：

1. \(0\leq {a_i\over i}<{\sqrt 5+1\over 2}\)
2. 任意非零的两个不同数 \(a_i, a_j\) ，满足 \(a_i\neq a_j\)

---

【分析】

方便起见，直接记 \(a_i\) 的上限为 \(\displaystyle \lfloor{{\sqrt 5+1\over 2}\cdot i}\rfloor\to b_i-1\)

考虑进行 dp ，设 \(f_{i, j}\) 表示前 \(i\) 个数中，有 \(j\) 个 \(0\) 的方案数

很显然，可以有 dp 转移：

当当前数为 \(0\) 时，方案数为 \(f_{i-1, j-1}\)

当当前数不为 \(0\) 时，则当前数不能与前面的 \((i-1)-j\) 个非 \(0\) 数相同，故方案数为 \(\lfloor{{\sqrt 5+1\over 2}\cdot i}\rfloor-((i-1)-j)=b_i-i+j\to c_i+j\)

因此，总方案数为 \(f_{i, j} = f_{i-1, j-1} + f_{i-1, j} \cdot (c_i+j)=(f_{i-1, j-1}+j\cdot f_{i-1, j}) + c_i\cdot f_{i-1, j}\)

这个转移方程的两部分等价于告诉我们，在转移过程中，\(f\) 的值有两种决策，一部分是 \(f_{i-1, j-1}+j\cdot f_{i-1, j}\) ，一部分是 \(c_i\cdot f_{i-1, j}\)

---

我们先考虑第一部分的决策：

其等价于 \(c_i=0\) ，则方程满足 \(f_{i, j} = f_{i-1, j-1}+j\cdot f_{i-1, j}\) ，即第二类斯特林数 \(\begin{Bmatrix} i \\ j \end{Bmatrix} = \begin{Bmatrix} i-1 \\ j-1 \end{Bmatrix} +j \begin{Bmatrix} i-1 \\ j \end{Bmatrix}\)

因此在题目给定 \(n\) 的情况下，我们需要对 \(0\) 的个数进行求和，即对应第二类斯特林数的行和，也就是贝尔数 \(Bell_n\)

因此，该部分的生成函数与贝尔数一致，为 \(\hat G(x)=\exp(e^x-1)\)

---

在第一部分的决策确定的情况下，我们考虑第二部分的决策：

根据上述的描述，我们可以变形式子：

\(\displaystyle \sum_j f_{i, j} =\sum_j (f_{i-1, j-1}+j\cdot f_{i-1, j}) + \sum_j c_i\cdot f_{i-1, j}\)

则 \(\displaystyle F_i=[{x^i\over i!}]\hat G(x)+c_i\cdot F_{i-1}\)

这个式子的组合意义表示，对于 \(i\) 个数的答案而言，它要不然走第二类斯特林数部分的贡献（答案为贝尔数），要不然走 \(c_i\) 的累乘部分（答案为 \(c_1, c_2, \cdots, c_i\) 的某个子集的累乘）

因此，原答案的生成函数为 \(\displaystyle [x^i]F(x)=\sum_{p+q=i}[{x^p\over p!}]\hat G(x)\cdot [x^q]\prod_{j=1}^i(1+c_jx)\)

我们考虑将 \(\hat G(x)\) 的 EGF 转 OGF ，再记 \(\displaystyle P_{[l, r]}(x)=\prod_{i=l}^r (1+c_ix)\)

则答案的生成函数优化为 \(\displaystyle [x^i]F(x)=\sum_{p+q=i}[x^p]G(x)\cdot [x^q]P_{[1, i]}(x)=[x^i](G\cdot P_{[l, r]})\)

---

当然，先 \(O(n\log n)\) 预处理 \(G(x)\) ，再对每次询问 \(n\) ，用启发式合并 FFT 或线段树分治 FFT 暴力计算 \([x^n](G\cdot P_{[1, n]})\) ，复杂度是 \(O(n\log^2 n)\) 的

所以最坏情况下，答案会劣化为 \(O(n^2\log^2 n)\) ，是不能通过本题的

---

我们考虑最后实际需要的答案为 \([x^n](G\cdot P_{[1, n]})\)

我们可以先通过线段树分治 FFT ，求出每个线段树端点 \([l, r]\) 区间内的答案 \(P_{[l, r]}(x)\)

现在我们考虑通过 cdq 分治 FFT 计算 \([1, n]\) 的答案：

对于左半部分的答案，它可以沿着线段树分治 FFT 的结点继续递归求解

对于右半部分的答案，它们的贡献中，一定有一部分来自于 \(G\cdot P_{[1, mid]}\)

例如对于右半部分的某个位置 \(n\) ，它的贡献 \([x^n](G\cdot P_{[1, n]})=[x^n]((G\cdot P_{[1, mid]})\cdot P_{[mid+1, n]})\)

故我们计算右半部分时，可以预算出 \(G\cdot P_{[1, mid]}\) ，再递归下去和右半部分的卷积进行求解

但是，这样的做法，复杂度仍然为 \(O(n^2\log^2 n)\)

为保证复杂度，考虑到右半部分的数值，卷积时只取答案的 \([x^{mid+1}]\) 到 \([x^n]\) 系数；因此我们可以计算出 \(G\cdot P_{[1, mid]}\) 后，直接删去次数较低的 \(mid+1\) 项

这样的操作下，能保证递归时，多项式长度减半

而对于左半部分，为保证复杂度，同理需要删除次数高于 \(mid\) 项的所有系数，使多项式长度减半

递归到叶子结点 \(n\) 后，则当前维护的多项式为：\([x^{n-1}, x^n](G\cdot P_{[1, n-1]})\) ，较低次项与较高次项已经全部删除了

因此，最后答案为：

\(\begin{aligned} &[x^n](G\cdot P_{[1, n]}) \\=&[x^{n-1}](G\cdot P_{[1, n-1]})\cdot [x^1](1+c_ix)+[x^n](G\cdot P_{[1, n-1]})\cdot [x^0](1+c_ix) \\=&[x^{n-1}](G\cdot P_{[1, n-1]})\cdot c_i+[x^n](G\cdot P_{[1, n-1]}) \end{aligned}\)

可以 \(O(1)\) 求出

因此，时间复杂度为 \(T(n)=2T({n\over 2})+O(n\log n)=O(n\log^2 n)\)

---

【代码】

为了方便计算偏移量，在 cdq 分治 FFT 时，数组采用了 0 开始的下标（也许可能是不需要的？）

#include <bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define sz(a) (int)a.size()
#define de(a) cout << #a <<" = "<<a<<endl
#define dd(a) cout << #a <<" = "<<a<<" "
#define all(a) a.begin(), a.end()
#define pw(x) (1ll<<(x))
#define lc(x) ((x)<<1)
#define rc(x) ((x)<<1|1)
#define rsz(a, x) (a.resize(x))
typedef unsigned uint;
typedef unsigned long long ull;
typedef long long ll;
typedef pair<int, int> pii;
typedef vector<int> vi;
typedef long double db;

const uint P=998244353;
inline uint kpow(uint a, ull x, uint p=P) { uint ans=1; for(; x; x>>=1, a=(ull)a*a%p) if(x&1) ans=(ull)ans*a%p; return ans; }
inline int exgcd(int a, int b, int &x, int &y) {
	static int g;
	return b?(exgcd(b, a%b, y, x), y-=a/b*x, g):(x=1, y=0, g=a);
}
inline int inv(int a, int p=P) {
	static int x, y;
	return (exgcd(a, p, x, y)==1)?(x<0?x+p:x):(-1);
}
const int LimBit=18;
const int M=1<<LimBit<<1;

namespace Poly {
	const uint G=3;
	inline uint norm(uint v, uint p=P) { return v>=p?v-p:v; }
	struct vir {
		uint v;
		inline vir(uint v_=0):v(norm(v_)) {}
		inline vir& operator += (const vir &x) { v=norm(v+x.v); return *this; }
		inline vir& operator -= (const vir &x) { v=norm(v+P-x.v); return *this; }
		inline vir& operator *= (const vir &x) { v=(ull)v*x.v%P; return *this; }
		inline vir operator + (const vir &x) const { vir y=*this; return y+=x; }
		inline vir operator - (const vir &x) const { vir y=*this; return y-=x; }
		inline vir operator * (const vir &x) const { vir y=*this; return y*=x; }
		
		inline vir operator - () const { return vir(P-v); }
		inline vir operator ! () const { return vir(inv(v)); }
		inline operator uint() const { return v; }
	};
	typedef vir Z;
	struct poly : public vector<Z> {
		inline friend ostream& operator << (ostream& out, const poly &p) {
			if(!p.empty()) out<<(uint)p[0];
			for(int i=1; i<sz(p); ++i) out<<' '<<(uint)p[i];
			return out;
		}
	};
	
	int N;
	Z Inv[M];
	vir invN, w[M>>1];
	inline void init() {
		Inv[1]=1;
		for(int i=2; i<M; ++i)
			Inv[i]=-Z(P/i)*Inv[P%i];
		
		w[0]=1; w[M>>2]=kpow(G, (P-1)/M);
		for(int i=M>>3; i; i>>=1) w[i]=w[i<<1]*w[i<<1];
		for(int i=1; i<(M>>1); ++i) w[i]=w[i&(i-1)]*w[i&-i];
	}
	
	inline void dft(vir *a) {
		for(int i=N>>1; i; i>>=1)
			for(vir *j=a, *o=w, x; j!=a+N; j+=i<<1, ++o)
				for(vir *k=j; k!=j+i; ++k)
					x=*o*k[i], k[i]=*k-x, *k+=x;
	}
	inline void idft(vir *a) {
		for(int i=1; i<N; i<<=1)
			for(vir *j=a, *o=w, x; j!=a+N; j+=i<<1, ++o)
				for(vir *k=j; k!=j+i; ++k)
					x=*k+k[i], k[i]=*o*(*k-k[i]), *k=x;
		reverse(a+1, a+N);
		for(int i=0; i<N; ++i) a[i]*=invN;
	}
	
	vir p1[M], p0[M];
	inline void get_mul(poly &a, poly &b, int na, int nb) {
		na=min(na, sz(a)); nb=min(nb, sz(b));
		int d=__lg(na+nb-2)+1; N=pw(d);
		invN=!vir(N);
		for(int i=0; i<na; ++i) p1[i]=(vir)a[i]; for(int i=na; i<N; ++i) p1[i]=vir();
		for(int i=0; i<nb; ++i) p0[i]=(vir)b[i]; for(int i=nb; i<N; ++i) p0[i]=vir();
		dft(p1); dft(p0);
		for(int i=0; i<N; ++i) p1[i]*=p0[i];
		idft(p1);
		rsz(a, na+nb-1); for(int i=0; i<sz(a); ++i) a[i]=p1[i];
	}
	
	poly a, b;
	inline void get_inv(poly &f, poly &g, int n) {
		b=f; rsz(b, n);
		rsz(g, 1); g[0]=!b[0];
		for(int l=2; l<n<<1; l<<=1) {
			get_mul(a=g, g, l>>1, l>>1); rsz(a, l);
			get_mul(a, b, l, l);
			rsz(g, l);
			for(int i=0; i<l; ++i) g[i]=g[i]+g[i]-a[i];
		}
		rsz(g, n);
	}
	
	inline void get_der(poly &f, poly &g) {
		g=f;
		for(int i=0; i<sz(g); ++i) g[i]*=Z(i);
		g.erase(g.begin());
	}
	inline void get_int(poly &f, poly &g, int C=0) {
		g=f; g.insert(g.begin(), Z(C));
		for(int i=1; i<sz(g); ++i) g[i]*=Inv[i];
	}
	
	poly c;
	inline void get_ln(poly &f, poly &g, int n, int ln0=0) {
		get_inv(f, c, n);
		get_der(f, a);
		get_mul(c, a, n, n);
		get_int(c, g, ln0);
		rsz(g, n); rsz(c, 0);
	}
	
	poly d, e;
	inline void get_exp(poly &f, poly &g, int n, int exp0=1) {
		d=f; rsz(d, n);
		rsz(g, 1); g[0]=exp0;
		for(int l=2; l<n<<1; l<<=1) {
			get_ln(g, e, l);
			e.erase(e.begin(), e.begin()+(l>>1));
			for(int i=0, j=(l>>1); j<l && j<n; ++i, ++j) e[i]-=d[j];
			get_mul(e, g, l>>1, l>>1);
			rsz(g, l); for(int i=l>>1, j=0; i<l; ++i, ++j) g[i]-=e[j];
		}
		rsz(d, 0); rsz(e, 0); rsz(g, n);
	}
	
	inline void gen_exp(poly &f, int n) {
		Z now=1;
		rsz(f, n);
		for(int i=0; i<n; ++i) {
			f[i]=now;
			now*=Inv[i+1];
		}
	}
	
	poly T[M], smt[M+M];
	inline poly& get_smt(int id, int l, int r) {
		poly &now = smt[id];
		if(l==r) {
			now=T[l];
			return now;
		}
		int mid=(l+r) / 2;
		poly &lcp = get_smt(lc(id), l, mid);
		poly &rcp = get_smt(rc(id), mid+1, r);
		get_mul(now=lcp, rcp, sz(lcp), sz(rcp));
		return now;
	}
	
	poly flr[M+M];
	void work_cdq(int id, poly &g, int l, int r) {
		if(l>r) return ;
		if(l==r) {
			b[l]=g[0] * T[l][1] + g[1];
			return ;
		}
		flr[id]=g;
		int mid = (l + r) / 2;
		g.erase(g.begin()+mid-l+2, g.end());
		work_cdq(lc(id), g, l, mid);
		g=flr[id];
		get_mul(g, smt[lc(id)], sz(g), sz(smt[lc(id)]));
		g.erase(g.begin(), g.begin()+mid-l+1);
		work_cdq(rc(id), g, mid+1, r);
		g=flr[id];
	}
	inline void get_cdq(poly &g, int n) {
		for(int i=0; i<n; ++i) b[i]=Z();
		work_cdq(1, g, 0, n-1);
		b.insert(b.begin(), Z());
		g=b;
	}
}
using Poly::poly;
using Poly::Z;

const int Lim=1e5, MAXN=Lim+10;
int a[MAXN];
poly f, g;
inline void work(int n=Lim) {
	Poly::gen_exp(f, n+1);
	f[0]-=Z(1);
	Poly::get_exp(f, f, n+1);
	
	Z fac=1;
	for(int i=0; i<=n; ++i) {
		f[i] *= fac;
		fac *= Z(i+1);
	}
	
	for(int i=0; i<n; ++i) {
		poly &now = Poly::T[i];
		rsz(now, 2);
		now[0]=1;
		now[1]=a[i+1];
	}
	Poly::get_smt(1, 0, n-1);
	
	Poly::get_cdq(f, n);
}
inline void init() {
	Poly::init();
	db phi=(1+sqrtl(5))/2;
	for(int i=1; i<=Lim; ++i) a[i]=ceill(phi*i)+1e-6;
	for(int i=1; i<=Lim; ++i) a[i]-=i;
}
int main() {
	ios::sync_with_stdio(0);
	cin.tie(0); cout.tie(0);
	init();
	work();
	int T, n; cin>>T;
	while(T--&&cin>>n)
		cout<<f[n]<<"\n";
	return 0;
}