Newcoder 146 B.Filling pools（CDQ分治+生成函数+FFT）

Description

对于一个$n\times n$的池塘，初始每个格子都没有水，每行和每列只能选取一个位置填满水，每个时刻对于一个没有水的位置，只要其上下左右四个位置中至少有两个有水，那么这个位置也会填满水，问有多少种选取方案可以使得整个水池最终填满水

Input

一个整数$n(1\le n<262144)$

Output

输出方案数，结果模$988244353$

Sample Input

3

Sample Output

6

Solution1

对于$n\times n$的区域，假设$n>1$，那么一个合法的方案可以按以下两种方式被分成两个合法方案的并：

1.左上方$x\times x$是合法方案和右下方$(n-x)\times (n-x)$是合法方案

2.右上方$x\times x$是合法方案和左下方$(n-x)\times (n-x)$是合法方案

设$f(n)$为$n\times n$区域的合法方案数，$g(n)$表示把$n\times n$区域划分成一个$x\times x$合法区域加上一个$(n-x)\times (n-x)$合法区域的方案数，由于第一种划分方式和第二种划分方式一一对应，故有$g(n)=\frac{f(n)+[n=1]}{2}$

由于一个$n\times n$区域的合法方案有多重划分方案，为避免记重，我们取$x$最小的划分方案，此时这个$x\times x$区域的合法划分方案的划分方式要与$n\times n$区域的划分方式不同（即一个是按主对角线划分另一个就要按斜对角线划分），进而该$x\times x$区域的划分方案数即为$g(x)$，而$(n-x)\times (n-x)$区域的划分方案即为$f(n-x)$，进而有转移$g(n)=\sum _{i=1}^{n-1}g(i)f(n-i)$

$CDQ$分治+$FFT$，为避免除法令$g(n)=2g(n)$，那么有$g(n)=f(n)+[n=1]$，假设当前考虑区间为$[l,r]$，目的是求$f(n),g(n),n\in [l,r]$，将区间二分为$[l,mid]\cup [mid+1,r]$，假设已经求出$f(x),g(x),x\in [l,mid]$，考虑这些值对$f[mid+1],...,f(r)$的贡献

令$t\in [mid+1,r]$，则$f(t)=\sum _{i=1}^{t-1}g(i)f(t-i)$，贡献分成两部分

1.$f(l),...,f(mid)$对$f(t)$的贡献，$\sum _{i=l}^{mid}f(i)g(t-i)$，此时需要$g(1),...,g(r-l)$

2.$g(l),...,g(mid)$对$f(t)$的贡献，$\sum _{i=l}^{mid}g(i)f(t-i)$，此时需要$f(1),...,f(r-l)$

若$r-l\ge l$，那么两部分贡献可以合并为$f(l),...,f(mid),g(l),...,g(mid)$对$f(t)$的贡献，否则两部分贡献要单独计算，贡献显然为卷积形式，用$FFT$加速即可，时间复杂度$O(nlo{g}^{2}n)$

Code1

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;
#define maxfft 524288+5
#define mod 998244353
int mul(int x,int y)
{
	ll z=1ll*x*y;
	return z-z/mod*mod;
}
int add(int x,int y)
{
	x+=y;
	if(x>=mod)x-=mod;
	return x;
}
const double pi=acos(-1.0);
struct cp
{
	double a,b;
	cp operator +(const cp &o)const {return (cp){a+o.a,b+o.b};}
	cp operator -(const cp &o)const {return (cp){a-o.a,b-o.b};}
	cp operator *(const cp &o)const {return (cp){a*o.a-b*o.b,b*o.a+a*o.b};}
	cp operator *(const double &o)const {return (cp){a*o,b*o};}
	cp operator !() const{return (cp){a,-b};}
}w[maxfft];
int pos[maxfft];
void fft_init(int len)
{
	int j=0;
	while((1<<j)<len)j++;
	j--;
	for(int i=0;i<len;i++)
		pos[i]=pos[i>>1]>>1|((i&1)<<j);
}
void fft(cp *x,int len,int sta)
{
	for(int i=0;i<len;i++)
		if(i<pos[i])swap(x[i],x[pos[i]]);
	w[0]=(cp){1,0};
	for(int i=2;i<=len;i<<=1)
	{
		cp g=(cp){cos(2*pi/i),sin(2*pi/i)*sta};
		for(int j=i>>1;j>=0;j-=2)w[j]=w[j>>1];
		for(int j=1;j<i>>1;j+=2)w[j]=w[j-1]*g;
		for(int j=0;j<len;j+=i)
		{
			cp *a=x+j,*b=a+(i>>1);
			for(int l=0;l<i>>1;l++)
			{
				cp o=b[l]*w[l];
				b[l]=a[l]-o;
				a[l]=a[l]+o;
			}
		}
	}
	if(sta==-1)for(int i=0;i<len;i++)x[i].a/=len,x[i].b/=len;
}
cp x[maxfft],y[maxfft],z[maxfft];
int temp[maxfft];
void FFT(int *a,int *b,int n,int m,int *c)
{
	if(n<=100&&m<=100||min(n,m)<=5)
	{
		for(int i=0;i<n+m-1;i++)temp[i]=0;
		for(int i=0;i<n;i++)
			for(int j=0;j<m;j++)
				temp[i+j]=add(temp[i+j],mul(a[i],b[j]));
		for(int i=0;i<n+m-1;i++)c[i]=temp[i];
		return ;
	}
	int len=1;
	while(len<n+m)len<<=1;
	fft_init(len);
	for(int i=0;i<len;i++)
	{
		int aa=i<n?a[i]:0,bb=i<m?b[i]:0;
		x[i]=(cp){(double)(aa>>15),(double)(aa&32767)},y[i]=(cp){(double)(bb>>15),(double)(bb&32767)};
	}
	fft(x,len,1),fft(y,len,1);
	for(int i=0;i<len;i++)
	{
		int j=len-1&len-i;
		z[i]=((x[i]+!x[j])*(y[i]-!y[j])+(x[i]-!x[j])*(y[i]+!y[j]))*(cp){0,-0.25};
	}
	fft(z,len,-1);
	for(int i=0;i<n+m-1;i++)
	{
		ll ta=(ll)(z[i].a+0.5)%mod;
		ta=(ta<<15)%mod;
		c[i]=ta;
	}
	for(int i=0;i<len;i++)
	{
		int j=len-1&len-i;
		z[i]=(x[i]-!x[j])*(y[i]-!y[j])*(cp){-0.25,0}+(x[i]+!x[j])*(y[i]+!y[j])*(cp){0,0.25};
	}
	fft(z,len,-1);
	for(int i=0;i<n+m-1;i++)
	{
		ll ta=(ll)(z[i].a+0.5)%mod,tb=(ll)(z[i].b+0.5)%mod;
		ta=(ta+(tb<<30))%mod;
		c[i]=add(c[i],ta);
	}
}
int f[maxfft],g[maxfft],h[maxfft];
void update(int l1,int r1,int l2,int r2,int l3,int r3)
{
	FFT(f+l1,g+l2,r1-l1+1,r2-l2+1,h);
	for(int i=0;i<=r1-l1+r2-l2&&l1+l2+i<=r3;i++)
		if(l1+l2+i>=l3)f[l1+l2+i]=add(f[l1+l2+i],h[i]);
}
void Solve(int l,int r)
{
	if(l==r)
	{
		if(l>1)g[l]=f[l];
		return ;
	}
	int mid=(l+r)/2;
	Solve(l,mid);
	if(r-l>=l)update(l,mid,l,mid,mid+1,r);
	else update(1,r-l,l,mid,mid+1,r),update(l,mid,1,r-l,mid+1,r);
	Solve(mid+1,r);
}
int main()
{
	int n;
	scanf("%d",&n);
	f[1]=1,g[1]=2;
	Solve(1,n);
	printf("%d\n",f[n]);
}

Solution2

生成函数+多项式开方，令$F(x)=\sum _{n=1}^{\infty }f(n)x^n,G(x)=\sum _{n=1}^{\infty }g(n)x^n$，那么有$G(x)=F(x)+\frac{x}{2}$，而上面的转移$g(n)=\sum _{i=1}^{n-1}g(i)f(n-i)$即为$[x^n]G(x)=[x^n](G(x)F(x)),n>1$，而$[x]G(x)=[x](G(x)F(x))+1$，故有$G(x)=G(x)F(x)+x$，将$G(x)=F(x)+\frac{x}{2}$代入得$F^2(x)+(x-1)F(x)+x=0$，解该一元二次方程得$F(x)=\frac{1-x\pm \sqrt{x^2-6x+1}}{2}$，注意到多项式开方的常数项为$1$，而$F(x)$没有常数项，故舍去另一个解，得到$F(x)=\frac{1-x-\sqrt{x^2-6x+1}}{2}$，时间复杂度$O(nlogn)$

Code2

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;
#define maxfft 524288+5
#define mod 998244353
#define inv2 499122177
const double pi=acos(-1.0);
struct cp
{
	double a,b;
	cp operator +(const cp &o)const {return (cp){a+o.a,b+o.b};}
	cp operator -(const cp &o)const {return (cp){a-o.a,b-o.b};}
	cp operator *(const cp &o)const {return (cp){a*o.a-b*o.b,b*o.a+a*o.b};}
	cp operator *(const double &o)const {return (cp){a*o,b*o};}
	cp operator !() const{return (cp){a,-b};}
}w[maxfft];
int pos[maxfft];
void fft_init(int len)
{
	int j=0;
	while((1<<j)<len)j++;
	j--;
	for(int i=0;i<len;i++)
		pos[i]=pos[i>>1]>>1|((i&1)<<j);
}
void fft(cp *x,int len,int sta)
{
	for(int i=0;i<len;i++)
		if(i<pos[i])swap(x[i],x[pos[i]]);
	w[0]=(cp){1,0};
	for(int i=2;i<=len;i<<=1)
	{
		cp g=(cp){cos(2*pi/i),sin(2*pi/i)*sta};
		for(int j=i>>1;j>=0;j-=2)w[j]=w[j>>1];
		for(int j=1;j<i>>1;j+=2)w[j]=w[j-1]*g;
		for(int j=0;j<len;j+=i)
		{
			cp *a=x+j,*b=a+(i>>1);
			for(int l=0;l<i>>1;l++)
			{
				cp o=b[l]*w[l];
				b[l]=a[l]-o;
				a[l]=a[l]+o;
			}
		}
	}
	if(sta==-1)for(int i=0;i<len;i++)x[i].a/=len,x[i].b/=len;
}
cp x[maxfft],y[maxfft],z[maxfft];
void FFT(int *a,int *b,int n,int m,int *c)
{
	int len=1;
	while(len<n+m)len<<=1;
	fft_init(len);
	for(int i=0;i<len;i++)
	{
		int aa=i<n?a[i]:0,bb=i<m?b[i]:0;
		x[i]=(cp){(double)(aa>>15),(double)(aa&32767)},y[i]=(cp){(double)(bb>>15),(double)(bb&32767)};
	}
	fft(x,len,1),fft(y,len,1);
	for(int i=0;i<len;i++)
	{
		int j=len-1&len-i;
		z[i]=((x[i]+!x[j])*(y[i]-!y[j])+(x[i]-!x[j])*(y[i]+!y[j]))*(cp){0,-0.25};
	}
	fft(z,len,-1);
	for(int i=0;i<n+m-1;i++)
	{
		ll ta=(ll)(z[i].a+0.5)%mod;
		ta=(ta<<15)%mod;
		c[i]=ta;
	}
	for(int i=0;i<len;i++)
	{
		int j=len-1&len-i;
		z[i]=(x[i]-!x[j])*(y[i]-!y[j])*(cp){-0.25,0}+(x[i]+!x[j])*(y[i]+!y[j])*(cp){0,0.25};
	}
	fft(z,len,-1);
	for(int i=0;i<n+m-1;i++)
	{
		ll ta=(ll)(z[i].a+0.5)%mod,tb=(ll)(z[i].b+0.5)%mod;
		ta=(ta+(tb<<30))%mod;
		c[i]=(c[i]+ta)%mod;
	}
}
int temp1[maxfft],temp2[maxfft];
void Poly_Inv(int *poly,int n,int *ans)
{
	ans[0]=1;
	for(int i=2;i<=n;i<<=1)
	{
		FFT(poly,ans,i,i/2,temp1);
		FFT(ans,temp1+i/2,i/2,i/2,temp1);
		for(int j=0;j<i/2;j++)ans[j+i/2]=temp1[j]==0?0:mod-temp1[j];
	}
}
void Poly_Root(int *poly,int n,int *ans)
{
	ans[0]=1;
	for(int i=2;i<=n;i<<=1)
	{
		Poly_Inv(ans,i,temp2);
		FFT(ans,ans,i/2,i/2,ans);
		for(int j=0;j<i;j++)ans[j]=(ll)(ans[j]+poly[j])*inv2%mod;
		FFT(ans,temp2,i,i,ans);
		for(int j=i;j<2*i;j++)ans[j]=0;
	}
}
int mul(int x,int y)
{
	ll z=1ll*x*y;
	return z-z/mod*mod;
}
int add(int x,int y)
{
	x+=y;
	if(x>=mod)x-=mod;
	return x;
}
int f[maxfft],g[maxfft];
int main()
{
	int n;
	scanf("%d",&n);
	if(n==1)
	{
		printf("1\n");
		return 0;
	}
	int len=1;
	while(len<=n)len<<=1;
	f[0]=1,f[1]=mod-6,f[2]=1;
	Poly_Root(f,len,g);
	printf("%d\n",add(0,mod-mul(g[n],inv2)));
	return 0;
}