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]=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 ",1ll*(invA*B)[n]*fac[n-1]%mod); return 0; }