【XSY3344】连续段 DP 牛顿迭代 NTT

题目大意

　　对于一个长度为 (n) 的排列 (p)，我们称一个区间 ([l,r]) 是连续的当且仅当 ((max_{lleq ileq r}a_i)-(min_{lleq ileq r}a_i)=r-l)。

　　对于两个排列 (p_1,p_2)，我们称这两个排列是等价的，当且仅当他们的长度相同且连续区间的集合相同。

　　给你 (n)，对于 (1leq ileq n)，所有求 (i) 阶排列形成的等价类个数对 (p) 取模的值。

　　(nleq 100000,p=k imes2^{18}+1,kin mathbb {N},p) 是质数。

题解

　　对于一个区间 ([l,r])，如果这个区间的所有子区间都是连续的，就称这个区间为强连续区间。

　　每次我们找一个最大的强连续区间，把这个区间内所有点缩到一起。

　　如果找不到，就找一个最小的连续区间，把这个区间内所有点缩到一起。

　　这个过程就形成了一个树结构。

• 总共有 (n) 个叶子。
• 对于一个代表强连续区间的点，这个点的儿子个数 (geq 2)。
• 对于一个代表连续区间的点，这个点的儿子个数 (geq 4)。（可以发现，长度为 (leq 3) 的序列是一定有强连续区间的。）

　　每一个排列 (p) 对应了一棵树。

　　对于一棵树，可以构造出一个排列 (p)。

　　那么就只用对这棵树计数了。

　　显然可以 (O(n^2)) DP。

　　记 (f_i) 为长度为 (i) 的排列的方案数，(F(x)=sum_{igeq 0}f_ix^i)，那么就有：

[F(x)=frac{(F(x)^2+F(x)^4)}{1-F(x)}+x ]

　　用牛顿迭代法解这个方程即可。

[(F(x)^2+F(x)^4)frac{1}{1-F(x)}-F(x)+x=0\ G(y)=(y^2+y^4)frac{1}{1-y}-y+x\ G'(y)=(2y+4y^3)frac{1}{1-y}+(y^2+y^4)(-frac{1}{(1-y)^2})(-1)-1\ G(F(x))equiv 0pmod {x^n}\ G(F(x))equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x)) pmod {x^n}\ F(x)equiv F_0(x)-frac{G(F_0(x))}{G'(F_0(x))}pmod {x^n}\ ]

　　时间复杂度：(O(nlog n))

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<functional>
#include<cmath>
#include<vector>
#include<assert.h>
//using namespace std;
using std::min;
using std::max;
using std::swap;
using std::sort;
using std::reverse;
using std::random_shuffle;
using std::lower_bound;
using std::upper_bound;
using std::unique;
using std::vector;
typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef std::pair<int,int> pii;
typedef std::pair<ll,ll> pll;
void open(const char *s){
#ifndef ONLINE_JUDGE
	char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
void open2(const char *s){
#ifdef DEBUG
	char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
int rd(){int s=0,c,b=0;while(((c=getchar())<'0'||c>'9')&&c!='-');if(c=='-'){c=getchar();b=1;}do{s=s*10+c-'0';}while((c=getchar())>='0'&&c<='9');return b?-s:s;}
void put(int x){if(!x){putchar('0');return;}static int c[20];int t=0;while(x){c[++t]=x%10;x/=10;}while(t)putchar(c[t--]+'0');}
int upmin(int &a,int b){if(b<a){a=b;return 1;}return 0;}
int upmax(int &a,int b){if(b>a){a=b;return 1;}return 0;}
const int N=270000;
ll p,g;
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;
}
namespace ntt
{
	const int W=262144;
	int rev[N];
	int *w[20];
	void init()
	{
		ll s=fp(g,(p-1)/W);
		w[18]=new int[1<<17];
		w[18][0]=1;
		for(int i=1;i<W/2;i++)
			w[18][i]=w[18][i-1]*s%p;
		for(int i=17;i>=1;i--)
		{
			w[i]=new int[1<<(i-1)];
			for(int j=0;j<1<<(i-1);j++)
				w[i][j]=w[i+1][j<<1];
		}
	}
	void ntt(ll *a,int n,int t)
	{
		for(int 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(int i=2,l=1;i<=n;i<<=1,l++)
			for(int j=0;j<n;j+=i)
				for(int k=0;k<i/2;k++)
				{
					ll u=a[j+k];
					ll v=a[j+k+i/2]*w[l][k];
					a[j+k]=(u+v)%p;
					a[j+k+i/2]=(u-v)%p;
				}
		if(t==-1)
		{
			reverse(a+1,a+n);
			ll inv=fp(n,p-2);
			for(int i=0;i<n;i++)
				a[i]=a[i]*inv%p;
		}
	}
	void add(ll *a,ll *b,ll *c,int n,int m,int l)
	{
		static ll a1[N];
		int k=max(n,m);
		for(int i=0;i<=k;i++)
			a1[i]=0;
		for(int i=0;i<=n;i++)
			a1[i]=(a1[i]+a[i])%p;
		for(int i=0;i<=m;i++)
			a1[i]=(a1[i]+b[i])%p;
		for(int i=0;i<=l;i++)
			c[i]=a1[i];
	}
	void mul(ll *a,ll *b,ll *c,int n,int m,int l)
	{
		static ll a1[N],a2[N];
		int k=1;
		while(k<=n+m)
			k<<=1;
		for(int i=0;i<k;i++)
			a1[i]=a2[i]=0;
		for(int i=0;i<=n;i++)
			a1[i]=a[i];
		for(int i=0;i<=m;i++)
			a2[i]=b[i];
		ntt(a1,k,1);
		ntt(a2,k,1);
		for(int i=0;i<k;i++)
			a1[i]=a1[i]*a2[i]%p;
		ntt(a1,k,-1);
		for(int i=0;i<=l;i++)
			c[i]=(a1[i]+p)%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[N],a2[N];
		for(int i=0;i<n<<1;i++)
			a1[i]=a2[i]=0;
		for(int i=0;i<n;i++)
			a1[i]=a[i];
		for(int i=0;i<n>>1;i++)
			a2[i]=b[i];
		ntt(a1,n<<1,1);
		ntt(a2,n<<1,1);
		for(int i=0;i<n<<1;i++)
			a1[i]=a2[i]*(2-a1[i]*a2[i]%p)%p;
		ntt(a1,n<<1,-1);
		for(int i=0;i<n;i++)
			b[i]=(a1[i]+p)%p;
	}
}
int e[50],t;
int check(int x)
{
	for(int i=1;i<=t;i++)
		if(fp(x,(p-1)/e[i])==1)
			return 0;
	return 1;
}
void getg()
{
	int _=p-1;
	for(int i=2;i*i<=_;i++)
		if(_%i==0)
		{
			e[++t]=i;
			while(_%i==0)
				_/=i;
		}
	if(_!=1)
		e[++t]=_;
	for(int i=1;;i++)
		if(check(i))
		{
			g=i;
			return;
		}
}
void solve(ll *a,int n)
{
	if(n==1)
		return;
	solve(a,n>>1);
	static ll f1[N],f2[N],f3[N],f4[N],f5[N],a1[N],a2[N],a3[N],a4[N];
	for(int i=0;i<n;i++)
		a1[i]=(p-a[i])%p;
	a1[0]++;
	for(int i=0;i<n;i++)
		f1[i]=a[i];
	ntt::ntt(f1,n<<1,1);
	for(int i=0;i<n<<1;i++)
		f2[i]=f1[i]*f1[i]%p;
	ntt::ntt(f2,n<<1,-1);
	for(int i=0;i<n;i++)
		f4[i]=f2[i];
	ntt::ntt(f4,n<<1,1);
	for(int i=0;i<n<<1;i++)
		f3[i]=f4[i]*f1[i]%p;
	ntt::ntt(f3,n<<1,-1);
	for(int i=0;i<n<<1;i++)
		f4[i]=f4[i]*f4[i]%p;
	ntt::ntt(f4,n<<1,-1);
	for(int i=0;i<n;i++)
		f5[i]=f4[i];
	ntt::ntt(f5,n<<1,1);
	for(int i=0;i<n<<1;i++)
		f5[i]=f5[i]*f1[i]%p;
	ntt::ntt(f5,n<<1,-1);
	for(int i=0;i<n;i++)
		f1[i]=a[i];
	
	for(int i=0;i<n;i++)
		a1[i]=(-f1[i]+3*f2[i]-2*f3[i]+f4[i]-f5[i])%p;
	for(int i=1;i<n;i++)
		a1[i]=(a1[i]-2*f1[i-1]+f2[i-1])%p;
	a1[1]++;
	
	for(int i=0;i<n;i++)
		a2[i]=(4*f1[i]-2*f2[i]+4*f3[i]-3*f4[i])%p;
	a2[0]--;

	ntt::getinv(a2,a3,n);
	ntt::mul(a1,a3,a4,n-1,n-1,n-1);
		
	for(int i=0;i<n;i++)
		a[i]=(a[i]-a4[i]+p)%p;
}
int n;
ll s[N];
int main()
{
	open("b");
	scanf("%d%lld",&n,&p);
	getg();
	ntt::init();
	int k=1;
	while(k<=n)
		k<<=1;
	solve(s,k);
	for(int i=1;i<=n;i++)
	{
		s[i]=(s[i]+p)%p;
		printf("%lld
",s[i]);
	}
	return 0;
}