题解 C. Lucky Sequence 2020 银川 Regional
之前现场没来得及开这一题，现在挂 cf 上就补了
第二类斯特林数 + 贝尔数 + 多项式 exp + 线段树分治 FFT + cdq 分治 FFT
感觉出出来就没打算让选手在现场过了...
【大意】
询问有多少个长度为 \(n\) 的序列 \(a_1, a_2, \cdots, a_n\) 满足如下条件：
- \(0\leq {a_i\over i}<{\sqrt 5+1\over 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;
}
