Description
给定 \(n\) 个排成一排的一模一样的球，要从中选出一些球并对它们进行分组，每个组中球的个数只能为 1,2，求取出 \(1,2,...,k \le 2^{16}\) 组的方案数分别为多少。
Solution
考虑 dp，设 \(f_{i,j}\) 表示考虑了前 \(i\) 个球并且已经取出了 \(j\) 组的方案数。
第一种转移方式是枚举最后一个球的去向，显然有三种：丢弃，单独成组，和上一个一起成组，则转移方程为
\[f_{i,j}=f_{i-1,j}+f_{i-1,j-1}+f_{i-2,j-1}
\]
第二种转移方式是，考虑从中间断开，即将 \(n\) 拆分为 \(a+b\)，那么如果第 \(a,a+1\) 个球不在同一组，则可以分成互不相干的两组来处理；如果是连着的，则对 \(1 \sim a-1, a+2 \sim a+b\) 分别处理，组数上再加上一组即可，转移方程为
\[f_{a+b,i} = \sum_{j=1}^{i-1} f_{a,j}f_{b,i-j}+\sum_{j=1}^{i-2}f_{a-1,j}f_{b-1,i-j-1}
\]
考虑转化为多项式形式，设 \(F_n(x)=\sum_{i=0}^k f_{n,i} x^i\)，则上述方程可以转化为
\[F_n(x)=F_{n-1}(x)+xF_{n-1}(x)+xF_{n-2}(x)
\]
\[F_{a+b}(x)=F_a(x)F_b(x)+xF_{a-1}(x)F_{b-1}(x)
\]
考虑令 \(a=n,b=n\)，则第二种转移方案的递推式改写为
\[F_{2n}(x)=F_n^2(x)+xF_{n-1}^2(x)
\]
同理有
\[F_{2n-1}(x)=F_n(x)F_{n-1}(x)+xF_{n-1}(x)F_{n-2}(x)
\]
\[F_{2n-2}(x)=F_{n-1}^2(x)+xF_{n-2}^2(x)
\]
整理一下，我们现在拥有的多项式形式的递推式有
\[F_n(x)=F_{n-1}(x)+xF_{n-1}(x)+xF_{n-2}(x)
\]
\[F_{2n}(x)=F_n^2(x)+xF_{n-1}^2(x)
\]
\[F_{2n-1}(x)=F_n(x)F_{n-1}(x)+xF_{n-1}(x)F_{n-2}(x)
\]
\[F_{2n-2}(x)=F_{n-1}^2(x)+xF_{n-2}^2(x)
\]
于是显然我们需要构建一个三元组 \(G_n = (F_n,F_{n-1},F_{n-2})\) 作为一个最小的单位。
根据第一个式子，我们可以做到 \(G_n \to G_{n+1}\)，根据后三个式子可以做到 \(G_n \to G_{2n}\)。
于是我们只需要对 \(n\) 分解一下即可。如果当前 \(n\) 是奇数那么就 \(-1\) 并且添加一个 \(G_n \to G_{n+1}\) 操作；如果当前 \(n\) 是偶数那么就 \(\div 2\) 并且添加一个 \(G_n \to G_{2n}\) 操作。循环执行直到 \(n=2\)，最后构造初始状态并倒序复原所有操作即可。
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = 262150;
const int mod = 998244353;
int qpow(int p,int q) {return (q&1?p:1)*(q?qpow(p*p%mod,q/2):1)%mod;}
namespace cipolla {
inline int le(int x) {return qpow(x,(mod-1)/2);}
int w;
struct comp {
int x,y;
comp(int a=0,int b=0) {x=a;y=b;}
};
comp operator + (comp a,comp b) {return comp((a.x+b.x)%mod,(a.y+b.y)%mod);}
comp operator - (comp a,comp b) {return comp((a.x-b.x+mod)%mod,(a.y-b.y+mod)%mod);}
comp operator * (comp a,comp b) {return comp((a.x*b.x+a.y*b.y%mod*w)%mod,(a.x*b.y+a.y*b.x)%mod);}
comp operator ^ (comp a,int b) {comp o(1,0); for(;b;a=a*a,b>>=1) if(b&1) o=o*a; return o;}
int calc(int x) {
x%=mod;
int a;
while(true) {
a=rand();
w=(a*a-x+mod)%mod;
if(le(w)==mod-1) break;
}
comp s=comp(a,1)^((mod+1)/2);
return min(s.x,mod-s.x);
}
}
namespace po {
int rev[N],inv[N],w[N],sz;
void presolve(int l) {
int len=1;
sz=0;
while(len<l) len<<=1, ++sz;
for(int i=1;i<len;i++) {
inv[i]=(i==1?1:inv[mod%i]*(mod-mod/i)%mod);
rev[i]=(rev[i>>1]>>1)|((i&1)<<(sz-1));
}
int wn=qpow(3,(mod-1)/len);
w[len/2]=1;
for(int i=len/2+1;i<len;i++) w[i]=w[i-1]*wn%mod;
for(int i=len/2-1;i;i--) w[i]=w[i<<1];
}
int pre(int l) {int g; for(g=1;g<l;g<<=1); return g;}
void ntt(int *a,int o,int n) {
static unsigned long long s[N];
int t=sz-__builtin_ctz(n),x;
for(int i=0;i<n;i++) s[rev[i]>>t]=a[i];
for(int l=1;l<n;l<<=1) for(int i=0;i<n;i+=l<<1) for(int j=0;j<l;j++) {
x=s[i+j+l]*w[j+l]%mod;
s[i+j+l]=s[i+j]+mod-x;
s[i+j]+=x;
}
for(int i=0;i<n;i++) a[i]=s[i]%mod;
if(o) {
x=qpow(n,mod-2);
for(int i=0;i<n;i++) a[i]=a[i]*x%mod;
reverse(a+1,a+n);
}
}
void mult(int n,int *x,int *y,int *z) {
static int a[N],b[N];
int l=pre(n<<1);
for(int i=0;i<l;i++) {
a[i]=(i<n?x[i]:0);
b[i]=(i<n?y[i]:0);
}
ntt(a,0,l); ntt(b,0,l);
for(int i=0;i<l;i++) z[i]=a[i]*b[i]%mod;
ntt(z,1,l);
for(int i=n;i<l;i++) z[i]=0;
}
void inve(int len,int *a,int *b) {
if(len==1) *b=qpow(*a,mod-2);
else {
inve((len+1)/2,a,b);
static int c[N];
int n=pre(len<<1);
for(int i=0;i<n;i++) i<len?c[i]=a[i]:b[i]=c[i]=0;
ntt(b,0,n);
ntt(c,0,n);
for(int i=0;i<n;i++) b[i]=((b[i]+b[i]-b[i]*b[i]%mod*c[i])%mod+mod)%mod;
ntt(b,1,n);
for(int i=len;i<n;i++) b[i]=0;
}
}
void sqrt(int n,int *a,int *b) {
if(n==1) *b=cipolla::calc(*a);
else {
sqrt((n+1)/2,a,b);
static int c[N];
inve(n,b,c);
mult(n,a,c,c);
for(int i=0;i<n;i++) b[i]=(b[i]+c[i])*inv[2]%mod;
}
}
void deri(int n,int *a,int *b) {
for(int i=0;i<n-1;i++) b[i]=a[i+1]*(i+1)%mod;
b[n-1]=0;
}
void inte(int n,int *a,int *b) {
for(int i=n-1;i>0;--i) b[i]=a[i-1]*inv[i]%mod;
b[0]=0;
}
void loge(int n,int *a,int *b) {
static int c[N];
inve(n,a,b);
deri(n,a,c);
mult(n,b,c,b);
inte(n,b,b);
}
void expr(int n,int *a,int *b) {
if(n==1) *b=1;
else {
expr((n+1)/2,a,b);
static int c[N];
loge(n,b,c);
for(int i=0;i<n;i++) c[i]=(a[i]-c[i]+mod)%mod;
c[0]=(c[0]+1)%mod;
mult(n,b,c,b);
}
}
}
int n,k,a[N],b[N],c[N];
struct poly
{
vector<int> x;
poly operator + (const poly &o)
{
poly res;
for(int i=0;i<k && i<x.size() && i<o.x.size();i++)
{
res.x.push_back((x[i]+o.x[i])%mod);
}
return res;
}
poly operator * (const poly &o)
{
poly res;
for(int i=0;i<k;i++) a[i]=b[i]=c[i]=0;
for(int i=0;i<x.size();i++) a[i]=x[i];
for(int i=0;i<o.x.size();i++) b[i]=o.x[i];
po::mult(k,a,b,c);
for(int i=0;i<k;i++) res.x.push_back(c[i]);
return res;
}
poly shift()
{
poly res;
res.x.push_back(0);
for(int i=1;i<k;i++) res.x.push_back(x[i-1]);
return res;
}
};
poly nullpoly()
{
poly res;
for(int i=0;i<k;i++) res.x.push_back(0);
return res;
}
poly unitpoly()
{
poly res;
res.x.push_back(1);
for(int i=1;i<k;i++) res.x.push_back(0);
return res;
}
void print(const vector<int> &vec)
{
for(int i=1;i<min(1ull*k,1ull*vec.size());i++) cout<<vec[i]<<" ";
cout<<endl;
}
void print(const poly &p)
{
print(p.x);
}
signed main() {
ios::sync_with_stdio(false);
cin>>n>>k;
++k;
po::presolve((k+1)<<2);
poly f2,f1,f0;
f0=unitpoly();
f1=f0+f0.shift();
f2=f1+f1.shift()+f0.shift();
if(n<=2)
{
if(n==1) print(f1);
if(n==2) print(f2);
}
else
{
swap(f2,f0);
int tmp=n;
vector<int> vop; // Vector of Operations, 0-Double, 1-Shift
while(tmp>2)
{
if(tmp&1)
{
vop.push_back(1);
tmp--;
}
else
{
vop.push_back(0);
tmp>>=1;
}
}
reverse(vop.begin(),vop.end());
for(int op:vop)
{
if(op==0)
{
poly g2,g1,g0;
g0=f0*f0+(f1*f1).shift();
g1=f0*f1+(f1*f2).shift();
g2=f1*f1+(f2*f2).shift();
f0=g0;
f1=g1;
f2=g2;
}
else
{
poly tmp;
tmp=f0+f0.shift()+f1.shift();
f2=f1;
f1=f0;
f0=tmp;
}
}
print(f0);
}
return 0;
}