BZOJ3456城市规划

BZOJ3456

http://www.lydsy.com/JudgeOnline/problem.php?id=3456

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
using namespace std;
const int mod=1004535809,maxn=1<<18|1;
int p[maxn],gn[233];
int n;
int fac[maxn],inv[maxn],tp[maxn];
inline int fp(int a,int b){
	int res=1;
	while(b){
		if(b&1)res=1ll*res*a%mod;
		a=1ll*a*a%mod;b>>=1;
	}
	return res;
}
inline void init(){
	for(register int i=0;i<=23;++i)gn[i]=fp(3,(mod-1)/(1<<(i+1)));
	fac[0]=1;for(register int i=1;i<=n+1;++i)fac[i]=1ll*fac[i-1]*i%mod;
	inv[1]=1;for(register int i=2;i<=n+1;++i)inv[i]=mod-1ll*mod/i*inv[mod%i]%mod;
	inv[0]=1;for(register int i=1;i<=n+1;++i)inv[i]=1ll*inv[i]*inv[i-1]%mod;
	for(register int i=0;i<=n+1;++i)tp[i]=fp(2,(1ll*i*(i-1)/2)%(mod-1));
}
struct poly{
	int a[maxn],n;
	poly(){memset(a,0,sizeof(a));}
	inline int &operator[](const int x){return a[x];}
	inline void ntt(int d,int f){
		for(register int i=0;i<d;++i)if(i<p[i])swap(a[i],a[p[i]]);
		for(register int i=1,t=0,w,v;i<d;i<<=1,++t){
			for(register int j=0;j<d;j+=(i<<1)){
				w=1;
				for(register int k=j;k<i+j;++k,w=1ll*w*gn[t]%mod)
					v=1ll*w*a[i+k]%mod,a[i+k]=(a[k]-v+mod)%mod,a[k]=(a[k]+v)%mod;
			}
		}
		if(f==1)return;
		reverse(a+1,a+d);
		register int ny=fp(d,mod-2);
		for(register int i=0;i<d;++i)a[i]=1ll*a[i]*ny%mod;
	}
	friend inline poly operator*(poly &A,poly &B){
		register poly res;res.n=A.n+B.n;int d,lg2;
		for(d=1,lg2=0;d<=res.n;d<<=1)++lg2;
		for(register int i=0;i<d;++i)p[i]=(p[i>>1]>>1)^((i&1)<<(lg2-1));
		A.ntt(d,1);B.ntt(d,1);
		for(register int i=0;i<d;++i)res[i]=1ll*A[i]*B[i]%mod;
		res.ntt(d,-1);A.ntt(d,-1);B.ntt(d,-1);
		return res;
	}
}A,B,tmp,invA,Ans;
inline void poly_inv(int deg){
	if(deg==1){
		invA[0]=1;invA[0]=fp(A[0],mod-2);
		return ;
	}
	poly_inv((deg+1)>>1);
	register int d=1,lg2=0;for(d=1;d<=(deg<<1);d<<=1)++lg2;
	for(register int i=0;i<d;++i)p[i]=(p[i>>1]>>1)^((i&1)<<(lg2-1));
	for(register int i=0;i<d;++i)tmp[i]=i<deg?A[i]:0;
	for(register int i=(deg+1)>>1;i<d;++i)invA[i]=0;
	invA.ntt(d,1);tmp.ntt(d,1);
	for(register int i=0;i<d;++i)invA[i]=(mod+2ll*invA[i]%mod-1ll*invA[i]*invA[i]%mod*tmp[i]%mod)%mod;
	invA.ntt(d,-1);invA.n=deg-1;
}
int main(){
	scanf("%d",&n);init();A.n=B.n=n;
	for(register int i=0;i<=n;++i)A[i]=1ll*tp[i]*inv[i]%mod;
	for(register int i=1;i<=n;++i)B[i]=1ll*tp[i]*inv[i-1]%mod;
	poly_inv(n+1);
	printf("%d\n",1ll*(invA*B)[n]*fac[n-1]%mod);
	return 0;
}

原文地址：https://www.cnblogs.com/Stump/p/8454369.html