「Codeforces 755G」PolandBall and Many Other Balls

Codeforces 755G. PolandBall and Many Other Balls

题意

有 $n$ 个球编号为 $1,2,\dotsc,n$，一组可以是一个球 ${i}$ 或者是两个相邻的球 ${i,i+1}$

对于 $i=1,2,\dotsc,k$ 求 $n$ 个球划分成 $i$ 组的方案数，每个球至多在一个组内，并且可以不在任何一个组内

对 $998244353$ 取模

$n\le 10^9,k\le 2^{15}$

做法

显然有 $\mathcal O(nk)$ 的 dp 做法

令 $f_{i,j}$ 表示前 $i$ 个球分成 $k$ 组的方案数

转移是

$$
f_{i,j}=f_{i-1,j}+f_{i-1,j-1}+f{i-2,j-1}
$$

用生成函数表示第二维，有

$$
f_i(x)=f_{i-1}(x)+x\cdot f_{i-1}(x)+x\cdot f_{i-2}(x)
$$

矩阵

这个东西可以矩阵转移

$$
\begin{bmatrix}
x+1 & x \
1 & 0
\end{bmatrix}
\begin{pmatrix}
f_i(x) \
f_{i-1}(x)
\end{pmatrix}

\begin{pmatrix}
f_{i+1}(x) \
f_i(x)
\end{pmatrix}
$$

快速幂一下就好了

复杂度 $\mathcal O(k\log k\log n)$

倍增

上述算法写得不好就被卡常了，直接倍增常数小得多

大概要讨论一下最中间的几个球的情况，没写

复杂度还是 $\mathcal O(k\log k\log n)$

通项

写出递推关于 $z$ 特征方程

$$
z^2=(x+1)z+x
$$

其中系数是多项式

解得

$$
\begin{align}
z_0&=\frac{x+1-\sqrt{x^2+6x+1}}{2} \
z_1&=\frac{x+1+\sqrt{x^2+6x+1}}{2}
\end{align}
$$

因此 $f_n(x)$ 可以被表示为 $A z_0^n+B z_1^n$

根据

$$
\begin{cases}
f_0(x) = A+B =1 \
f_1(x) = A z_0 + B z_1 = x+1
\end{cases}
$$

解得

$$
\begin{cases}
A=\frac{-x-1+\sqrt{x^2+6x+1}}{2\sqrt{x^2+6x+1}} \
B=\frac{x+1+\sqrt{x^2+6x+1}}{2\sqrt{x^2+6x+1}}
\end{cases}
$$

事实上由于 $z_0$ 和 $A$ 的常数项为 $0$，$A z_0^n$ 对答案的前 $n$ 项没有影响，而答案除了前 $n$ 项都为 $0$，可以特判后忽略这部分

于是

$$
\begin{align}
f_n(x) &\equiv B z_1^n &\pmod{x^{n+1}} \
&= \frac{z_1^{n+1}}{\sqrt{x^2+6x+1}} &
\end{align}
$$

直接做复杂度也是 $\mathcal O(k\log k\log n)$，常数想必更小

由于这里 $z_1$ 的常数项为 $1$，可以方便地用 $\ln$ 和 $\exp$ 计算

总复杂度 $\mathcal O(k\log k)$

代码

矩乘

#include<cstdio>
#include<algorithm>
#include<cctype>
#include<string.h>
#include<cmath>
#include<vector>

using namespace std;
#define ll long long

const int N = 1<<15, M = N<<1, P = 998244353;
int n, k, w[M];
vector<int> ans;
inline int Pow(ll x, int y=P-2){
	int ans=1;
	for(; y; y>>=1, x=x*x%P) if(y&1) ans=ans*x%P;
	return ans;
}
inline void DFT(vector<int> &f, int n){
	static unsigned ll F[M];
	for(int i=0, j=0; i<n; ++i){
		F[i]=f[j];
		for(int k=n>>1; (j^=k)<k; k>>=1);
	}
	for(int i=1; i<n; i<<=1) for(int j=0; j<n; j+=i<<1) for(int k=j; k<j+i; ++k){
		int t=F[k+i]*w[i+k-j]%P;
		F[k+i]=F[k]+P-t, F[k]+=t;
	}
	for(int i=0; i<n; ++i) f[i]=F[i]%P;
}
inline void IDFT(vector<int> &f, int n){
	reverse(f.begin()+1, f.end()), DFT(f, n);
	for(int i=0, I=Pow(n); i<n; ++i) f[i]=(ll)f[i]*I%P;
}
inline int Calc(int x){ int ans=1; while(ans<=x) ans<<=1; return ans;}
inline vector<int> operator *(const vector<int> &x, const vector<int> &y){
	if(!x.size() || !y.size()){
		vector<int> a;
		return a.push_back(0), a;
	}
	vector<int> a=x, b=y;
	int n=Calc(x.size()+y.size()-2);
	a.resize(n), b.resize(n), DFT(a, n), DFT(b, n);
	for(int i=0; i<n; ++i) a[i]=(ll)a[i]*b[i]%P;
	IDFT(a, n);
	return a.resize(min((int)x.size()+(int)y.size()-1, k+1)), a;
}
inline void operator *=(vector<int> &x, const vector<int> &y){ x=x*y;}
inline vector<int> operator +(const vector<int> &x, const vector<int> &y){
	vector<int> ans=x;
	if(y.size()>x.size()) ans.resize(y.size());
	for(unsigned i=0; i<y.size(); ++i) (ans[i]+=y[i])%=P;//
	return ans;
}
inline void operator +=(vector<int> &x, const vector<int> &y){ x=x+y;}
struct matrix{
	vector<int> a[2][2];
	inline matrix operator *(const matrix &rhs)const{
		matrix ans;
		for(int i=0; i<2; ++i) for(int k=0; k<2; ++k) for(int j=0; j<2; ++j)
			ans.a[i][j]+=a[i][k]*rhs.a[k][j];
		return ans;
	}
} A;
inline matrix Pow(matrix x, int y){
	matrix ans=x;
	--y;
	for(; y; y>>=1, x=x*x) if(y&1) ans=ans*x;
	return ans;
}
int main() {
	for(int i=1; i<M; i<<=1){
		w[i]=1, w[i+1]=Pow(3, (P-1)/i/2);
		for(int j=2; j<i; ++j) w[i+j]=(ll)w[i+j-1]*w[i+1]%P;
	}
	scanf("%d%d", &n, &k);
	A.a[0][0]={1, 1}, A.a[0][1]={0, 1}, A.a[1][0]={1};
	A=Pow(A, n), ans=A.a[1][0]*vector<int>{1,1}+A.a[1][1];
	for(int i=1; i<=k; ++i) printf("%d ", i>n?0:ans[i]);
	return 0;
}

通项

纯手写

刚刚卡了下常数还是没最快= =

数组迭代实现好快= =

#include<cstdio>
#include<algorithm>
#include<cctype>
#include<string.h>
#include<cmath>
#include<vector>

using namespace std;
#define ll long long

const int OUT_LEN = 1000000;
char obuf[OUT_LEN], *ooh=obuf;
inline void print(char c) {
	if (ooh==obuf+OUT_LEN) fwrite(obuf, 1, OUT_LEN, stdout), ooh=obuf;
	*ooh++=c;
}
template<class T>
inline void print(T x) {
	static int buf[30], cnt;
	if (x==0) print('0');
	else {
		if (x<0) print('-'), x=-x;
		for (cnt=0; x; x/=10) buf[++cnt]=x%10+48;
		while(cnt) print((char)buf[cnt--]);
	}
}
inline void flush() { fwrite(obuf, 1, ooh - obuf, stdout); }

const int N = 1<<15, M = N<<1, P = 998244353;
int n, k, w[M], inv[M];
inline int Pow(ll x, int y=P-2){
	int ans=1;
	for(; y; y>>=1, x=x*x%P) if(y&1) ans=ans*x%P;
	return ans;
}
inline void DFT(vector<int> &f, int n){
	static unsigned ll F[M];
	for(int i=0, j=0; i<n; ++i){
		F[i]=f[j];
		for(int k=n>>1; (j^=k)<k; k>>=1);
	}
	for(int i=1; i<n; i<<=1) for(int j=0; j<n; j+=i<<1){
		int *W=w+i;
		unsigned ll *F0=F+j, *F1=F+j+i;
		for(int k=j; k<j+i; ++k, ++W, ++F0, ++F1){
			int t=*F1**W%P;
			*F1=*F0+P-t, *F0+=t;
		}
	}
	for(int i=0; i<n; ++i) f[i]=F[i]%P;
}
inline void IDFT(vector<int> &f, int n){
	reverse(f.begin()+1, f.end()), DFT(f, n);
	for(int i=0, I=Pow(n); i<n; ++i) f[i]=(ll)f[i]*I%P;
}
inline int Calc(int x){ int ans=1; while(ans<=x) ans<<=1; return ans;}
inline vector<int> operator *(const vector<int> &x, const vector<int> &y){
	if(!x.size() || !y.size()) return {0};
	if((unsigned ll)x.size()*y.size()<=1<<8){
		vector<int> ans(x.size()+y.size()-1);
		for(unsigned i=0; i<x.size(); ++i) for(unsigned j=0; j<y.size(); ++j)
			ans[i+j]=(ans[i+j]+(ll)x[i]*y[j])%P;
		return ans;
	}
	vector<int> a=x, b=y;
	int n=Calc(x.size()+y.size()-2);
	a.resize(n), b.resize(n), DFT(a, n), DFT(b, n);
	for(int i=0; i<n; ++i) a[i]=(ll)a[i]*b[i]%P;
	IDFT(a, n);
	return a.resize(min((int)x.size()+(int)y.size()-1, k+1)), a;
}
inline void operator *=(vector<int> &x, const vector<int> &y){ x=x*y;}
inline vector<int> operator +(const vector<int> &x, const vector<int> &y){
	vector<int> ans=x;
	if(y.size()>x.size()) ans.resize(y.size());
	for(unsigned i=0; i<y.size(); ++i) (ans[i]+=y[i])%=P;
	return ans;
}
inline void operator +=(vector<int> &x, const vector<int> &y){ x=x+y;}
inline vector<int> operator -(const vector<int> &x, const vector<int> &y){
	vector<int> ans=x;
	if(y.size()>x.size()) ans.resize(y.size());
	for(unsigned i=0; i<y.size(); ++i) (ans[i]+=P-y[i])%=P;
	return ans;
}
inline vector<int> PolyDiv2(const vector<int> &x){
	vector<int> ans(x.size());
	for(unsigned i=0; i<x.size(); ++i) ans[i]=(x[i]&1?x[i]+P:x[i])>>1;
	return ans;
}
inline vector<int> Ext(const vector<int> &a, int n){
	if(n<=(int)a.size()) return vector<int>(a.begin(), a.begin()+n);
	vector<int> ans=a;
	return ans.resize(n), ans;
}
vector<int> PolyInv(const vector<int> &a, int n=-1){
	if(n==-1) n=a.size();
	if(n==1) return {Pow(a[0])};
	vector<int> ans=PolyInv(a, (n+1)/2), tmp=Ext(a, n);
	int m=Calc(n*2-1);
	ans.resize(m), tmp.resize(m), DFT(ans, m), DFT(tmp, m);
	for(int i=0; i<m; ++i) ans[i]=(2+(ll)(P-tmp[i])*ans[i])%P*ans[i]%P;
	IDFT(ans, m);
	return ans.resize(n), ans;
}
vector<int> PolySqrt(const vector<int> &a, int n=-1){
	if(n==-1) n=a.size();
	if(n==1) return {1};
	vector<int> ans=PolySqrt(a, (n+1)/2);
	return PolyDiv2(Ext(ans+Ext(a, n)*PolyInv(Ext(ans, n)), n));
}
inline vector<int> D(const vector<int> &a){
	vector<int> ans(a.size()-1);
	for(unsigned i=1; i<a.size(); ++i) ans[i-1]=(ll)a[i]*i%P;
	return ans;
}
inline vector<int> Int(const vector<int> &a){
	vector<int> ans(a.size()+1);
	for(unsigned i=0; i<a.size(); ++i) ans[i+1]=(ll)a[i]*inv[i+1]%P;
	return ans;
}
inline vector<int> PolyLn(const vector<int> &a){
	return Int(Ext(D(a)*PolyInv(a), a.size()-1));
}
vector<int> PolyExp(const vector<int> &a, int n=-1){
	if(n==-1) n=a.size();
	if(n==1) return {1};
	vector<int> ans=PolyExp(a, (n+1)/2);
	return Ext(ans*(Ext(a, n)-PolyLn(Ext(ans, n))+vector<int>{1}), n);
}
int main() {
	for(int i=1; i<M; i<<=1){
		w[i]=1, w[i+1]=Pow(3, (P-1)/i/2);
		for(int j=2; j<i; ++j) w[i+j]=(ll)w[i+j-1]*w[i+1]%P;
	}
	inv[1]=1;
	for(int i=2; i<M; ++i) inv[i]=(ll)(P-P/i)*inv[P%i]%P;
	scanf("%d%d", &n, &k);
	vector<int> s={1, 6, 1}, c;
	s.resize(k+1), s=PolySqrt(s);
	c=PolyLn(PolyDiv2(vector<int>{1, 1}+s));
	for(int &i:c) i=(ll)i*(n+1)%P;
	c=PolyExp(c)*PolyInv(s);
	for(int i=1; i<=k; ++i) (i>n?print('0'):print(c[i])), print(' ');
	return flush(), 0;
}