【XSY2730】Ball 多项式exp 多项式ln 多项式开根 常系数线性递推 DP

题目大意

  一行有(n)个球,现在将这些球分成(k) 组,每组可以有一个球或相邻两个球。一个球只能在至多一个组中(可以不在任何组中)。求对于(1leq kleq m)的所有(k)分别有多少种分组方法。

  答案对(998244353)取模。

  (nleq {10}^9,m<2^{19})

题解

  因为(k>n)的项都是(0),所以我们钦定(mleq n)

  考虑DP。

  记(f_{i,j})为前(i)个球分为(j)组的方案数。

[f_{i,j}=f_{i-1,j}+f_{i-1,j-1}+f_{i-2,j-1} ]

  直接做是(O(nm))的。

  如果把(f_i)看成一个多项式,即

[F_i(x)=sum_{jgeq 0}f_{i,j}x^j ]

  那么转移就变成了

[F_i(x)=(1+x)F_{i-1}(x)+xF_{i-2}(x) ]

  这是一个常系数齐次线性递推,用FFT优化可以做到(O(mlog mlog n))

  考虑怎么求一个常系数齐次线性递推关系的通项公式。

  先求出这个转移矩阵的特征多项式:

[lambda^2-(1+x)lambda-x ]

  特征值为

[egin{align} lambda_1&=frac{1+x+sqrt{x^2+6x+1}}{2}\ lambda_2&=frac{1+x-sqrt{x^2+6x+1}}{2} end{align} ]

  我们钦定(F_{-1}(x)=0),设

[F_i(x)=c_1{lambda_1}^{i+1}+c_2{lambda_2}^{i+1} ]

  带入(F_{-1}(x),F_0(x))

[egin{cases} c_1&+c2&=0\ c_2lambda_1&+c_2lambda_2&=1 end{cases} ]

  解得

[egin{cases} c_1&=frac{1}{lambda_1-lambda_2}\ c_2&=frac{1}{lambda_2-lambda_1} end{cases} ]

  于是

[F_i(x)=frac{{lambda_1}^{i+1}-{lambda_2}^{i+1}}{lambda_1-lambda_2} ]

  直接用多项式开根求出(lambda_1,lambda_2),然后用多项式(ln exp)求出(F_n(x))

  时间复杂度:(O(mlog m))

  小优化:因为([x^0]lambda_2=0),所以(lambda_2)的前(n+1)项都是(0)。因为(mleq n),所以

[{lambda_2}^{n+1}equiv 0 pmod {x^{m+1}} ]

  所以我们不用计算({lambda_2}^{n+1})了。

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int maxn=300000;
const ll p=998244353;
const ll g=3;
const ll inv2=(p+1)/2;
ll fp(ll a,ll b)
{
	ll s=1;
	for(;b;b>>=1,a=a*a%p)
		if(b&1)
			s=s*a%p;
	return s;
}
int rev[maxn];
ll w1[maxn];
ll w2[maxn];
void ntt(ll *a,int n,int t)
{
	int i,j,k;
	ll u,v,w,wn;
	for(i=1;i<n;i++)
	{
		rev[i]=(rev[i>>1]>>1)|(i&1?n>>1:0);
		if(rev[i]>i)
			swap(a[i],a[rev[i]]);
	}
	for(i=2;i<=n;i<<=1)
	{
		wn=(t==1?w1[i]:w2[i]);
		for(j=0;j<n;j+=i)
		{
			w=1;
			for(k=j;k<j+i/2;k++)
			{
				u=a[k];
				v=a[k+i/2]*w%p;
				a[k]=(u+v)%p;
				a[k+i/2]=(u-v)%p;
				w=w*wn%p;
			}
		}
	}
	if(t==-1)
	{
		ll inv=fp(n,p-2);
		for(i=0;i<n;i++)
			a[i]=a[i]*inv%p;
	}
}
void getinv(ll *a,ll *b,int n)
{
	if(n==1)
	{
		b[0]=fp(a[0],p-2);
		return;
	}
	getinv(a,b,n>>1);
	static ll a1[maxn],a2[maxn];
	int i;
	for(i=0;i<n;i++)
		a1[i]=a[i];
	for(;i<n<<1;i++)
		a1[i]=0;
	for(i=0;i<n>>1;i++)
		a2[i]=b[i];
	for(;i<n<<1;i++)
		a2[i]=0;
	ntt(a1,n<<1,1);
	ntt(a2,n<<1,1);
	for(i=0;i<n<<1;i++)
		a1[i]=a2[i]*((2-a1[i]*a2[i])%p)%p;
	ntt(a1,n<<1,-1);
	for(i=0;i<n;i++)
		b[i]=a1[i];
}
void getsqrt(ll *a,ll *b,int n)
{
	if(n==1)
	{
		b[0]=1;
		return;
	}
	getsqrt(a,b,n>>1);
	static ll a1[maxn],a2[maxn],a3[maxn];
	int i;
	for(i=0;i<n>>1;i++)
		a1[i]=b[i];
	for(;i<n<<1;i++)
		a1[i]=0;
	for(i=0;i<n;i++)
		a2[i]=a[i];
	for(;i<n<<1;i++)
		a2[i]=0;
	getinv(a1,a3,n);
	ntt(a2,n<<1,1);
	ntt(a3,n<<1,1);
	for(i=0;i<n<<1;i++)
		a2[i]=a2[i]*a3[i]%p;
	ntt(a2,n<<1,-1);
	for(i=0;i<n;i++)
		b[i]=(a1[i]+a2[i])%p*inv2%p;
}
void getln(ll *a,ll *b,int n)
{
	static ll inv[maxn],a1[maxn],a2[maxn];
	int i;
	inv[0]=inv[1]=1;
	for(i=2;i<n;i++)
		inv[i]=-p/i*inv[p%i]%p;
	for(i=1;i<n;i++)
		a1[i-1]=a[i]*i%p;
	a1[n-1]=0;
	getinv(a,a2,n);
	for(i=n;i<n<<1;i++)
		a1[i]=a2[i]=0;
	ntt(a1,n<<1,1);
	ntt(a2,n<<1,1);
	for(i=0;i<n<<1;i++)
		a1[i]=a1[i]*a2[i]%p;
	ntt(a1,n<<1,-1);
	for(i=1;i<n;i++)
		b[i]=a1[i-1]*inv[i]%p;
	b[0]=0;
}
void getexp(ll *a,ll *b,int n)
{
	if(n==1)
	{
		b[0]=1;
		return;
	}
	getexp(a,b,n>>1);
	int i;
	for(i=n>>1;i<n;i++)
		b[i]=0;
	static ll a1[maxn],a2[maxn],a3[maxn];
	getln(b,a1,n);
	for(i=0;i<n>>1;i++)
	{
		a2[i]=b[i];
		a3[i]=a[i+(n>>1)]-a1[i+(n>>1)];
	}
	for(;i<n;i++)
		a2[i]=a3[i]=0;
	ntt(a2,n,1);
	ntt(a3,n,1);
	for(i=0;i<n;i++)
		a2[i]=a2[i]*a3[i]%p;
	ntt(a2,n,-1);
	for(i=0;i<n>>1;i++)
		b[i+(n>>1)]=a2[i];
}
ll a[maxn];
ll b[maxn];
ll c[maxn];
ll d[maxn];
ll e[maxn];
int n,m;
int k;
void solve()
{
	k=1;
	while(k<=m)
		k<<=1;
	int i;
	for(i=2;i<=k<<1;i++)
	{
		w1[i]=fp(g,(p-1)/i);
		w2[i]=fp(w1[i],p-2);
	}
	d[0]=1;
	d[1]=6;
	d[2]=1;
	getsqrt(d,c,k);
	for(i=0;i<k;i++)
		a[i]=c[i];
	a[0]++;
	a[1]++;
	for(i=0;i<k;i++)
		a[i]=a[i]*inv2%p;
	getinv(c,b,k);
	getln(a,e,k);
	for(i=0;i<k;i++)
		e[i]=e[i]*(n+1)%p;
	getexp(e,a,k);
	for(i=k;i<k<<1;i++)
		a[i]=b[i]=0;
	ntt(a,k<<1,1);
	ntt(b,k<<1,1);
	for(i=0;i<k<<1;i++)
		a[i]=a[i]*b[i]%p;
	ntt(a,k<<1,-1);
	for(i=1;i<=m;i++)
	{
		a[i]=(a[i]+p)%p;
		printf("%lld ",a[i]);
	}
}
int main()
{
#ifndef ONLINE_JUDGE
	freopen("b.in","r",stdin);
	freopen("b.out","w",stdout);
#endif
	int x=0;
	scanf("%d%d",&n,&m);
	if(m>n)
	{
		x=m-n;
		m=n;
	}
	solve();
	while(x--)
		printf("0 ");
	return 0;
}
原文地址:https://www.cnblogs.com/ywwyww/p/8513484.html